In [6]:
import telluricpy, numpy as np, gc
import scipy
import VTKUtil as pvtkUtil
In [7]:
%matplotlib qt
In [8]:
def simpeg2vtk(mesh,modDict):
from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION
from vtk.util.numpy_support import numpy_to_vtk
# Deal with dimensionalities
if mesh.dim >= 1:
vX = mesh.vectorNx
xD = mesh.nNx
yD,zD = 1,1
vY, vZ = np.array([0,0])
if mesh.dim >= 2:
vY = mesh.vectorNy
yD = mesh.nNy
if mesh.dim == 3:
vZ = mesh.vectorNz
zD = mesh.nNz
# Use rectilinear VTK grid.
# Assign the spatial information.
vtkObj = rectGrid()
vtkObj.SetDimensions(xD,yD,zD)
vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
# Assign the model('s) to the object
if modDict is not None:
for item in modDict.iteritems():
# Convert numpy array
vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
# Set the active scalar
vtkObj.GetCellData().SetActiveScalars(modDict.keys()[0])
return vtkObj
In [20]:
def plot3dSetion(mesh,mod,figName,lutName = 'Con', camera = None):
# Convert the mesh and model
fullvtkObj = simpeg2vtk(mesh,mod)
import VTKUtil as pvtkUtil
# Make scalar bar and lookup table
lutRes, scalarBarRes, vecNameRes = pvtkUtil.makeLookupTable(lutName)
scalarBarRes.SetPosition(0.01,0.15)
# Set some sizes
renSize = [1280,800]
axesBounds = [556800.0, 557800.0, 7133100.0, 7134100.0, -500.0, 500.0]
screenSize = 14.0
xRange = [556.8, 557.8]
yRange = [7133.1, 7134.1]
zRange = [-0.5, 0.5]
# Read the model
boxImp = vtk.vtkBox()
boxImp.SetBounds(556800.0, 557800.0, 7133100.0, 7134100.0, -500.0, 500.0)
extractFilt = vtk.vtkExtractGeometry()
extractFilt.SetExtractInside(1)
extractFilt.SetExtractBoundaryCells(1)
extractFilt.SetInputData(fullvtkObj)
extractFilt.SetImplicitFunction(boxImp)
# extractFilt.Update()
# Remove air
thresFilt = vtk.vtkThreshold()
thresFilt.SetInputConnection(extractFilt.GetOutputPort())
thresFilt.ThresholdByUpper(1e-8)
thresFilt.AllScalarsOn()
thresFilt.SetInputArrayToProcess(0,0,0,vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS,vecNameRes)
# thresFilt.Update()
vtkObjIn = thresFilt.GetOutput()
# print vtkObjIn.GetNumberOfCells()
# Set the camera
if camera is None:
camera = vtk.vtkCamera()
camera.SetClippingRange(85.01306538179968, 85013.06538179968)
camera.SetFocalPoint(557170.1657806067, 7133465.5306818895, -8.100855224539828)
camera.SetPosition(555046.658444768, 7131550.331097811, 1643.4398007476757)
camera.SetViewUp(0.35030191124890564, 0.3579734646304746, 0.8655308022224387)
camera.SetParallelScale(1.0)
else:
camera = camera
# Make a renderer
ren = vtk.vtkRenderer()
ren.SetActiveCamera(camera)
# Make renderwindow.
renwin = vtk.vtkRenderWindow()
# Set to off screen rendering
renwin.AddRenderer(ren)
renwin.SetSize(renSize)
iren = vtk.vtkRenderWindowInteractor()
iren.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
iren.SetRenderWindow(renwin)
iren.Initialize()
# Add the axes
axes = pvtkUtil.addAxes(screenSize,ren,xRange,yRange,zRange,axesBounds)
## Organize the data
# Plane 1
global plane, actor1
plane = vtk.vtkPlane()
plane.SetOrigin(557100,7133600,0)
plane.SetNormal(0,-1,0)
vtkObjClip1 = vtk.vtkClipDataSet()
# vtkObjClip1.SetInputData(vtkObjIn)
vtkObjClip1.SetInputConnection(thresFilt.GetOutputPort())
vtkObjClip1.SetClipFunction(plane)
vtkObjClip1.InsideOutOn()
# vtkObjClip1.Update()
vtkObj1 = vtkObjClip1.GetOutput()
vtkObj1.GetCellData().SetActiveScalars(vecNameRes)
# Set the mapper's
# Clip 1
mapper1 = vtk.vtkDataSetMapper()
mapper1.SetInputData(vtkObj1)
mapper1.SetScalarVisibility(1)
mapper1.SetLookupTable(lutRes)
mapper1.UseLookupTableScalarRangeOn()
mapper1.SetInterpolateScalarsBeforeMapping(1)
actor1 = vtk.vtkLODActor()
actor1.SetMapper(mapper1)
actor1.VisibilityOff()
# actor1.GetProperty().SetEdgeColor(1,0.5,0)
actor1.GetProperty().SetEdgeVisibility(0)
# actor1.SetScale(1.01, 1.01, 1.01)
# actor1.GetProperty().SetRepresentationToSurface()
if False:
# Create the widget, its representation, and callback
def MovePlane(widget, event_string):
rep.GetPlane(plane)
rep = vtk.vtkImplicitPlaneRepresentation()
rep.SetPlaceFactor(1.0);
rep.PlaceWidget(vtkObjClip1.GetOutput().GetBounds())
rep.DrawPlaneOn()
rep.SetOrigin(557100,7133600,0)
rep.SetNormal(0,-1,0)
# rep.SetPlane(plane)
planeWidget = vtk.vtkImplicitPlaneWidget2()
planeWidget.SetInteractor(iren)
planeWidget.SetRepresentation(rep)
planeWidget.SetEnabled(1)
# planeWidget.PlaceWidget()
planeWidget.AddObserver("InteractionEvent",MovePlane)
else:
# Callback function
def movePlane(obj, event):
global plane, actor1
obj.GetPlane(plane)
actor1.VisibilityOn()
# Associate the line widget with the interactor
planeWidget = vtk.vtkImplicitPlaneWidget()
# planeWidget.SetInputConnection(thresFilt.GetOutputPort())
planeWidget.SetInputConnection(vtkObjClip1.GetOutputPort())
planeWidget.SetInteractor(iren)
planeWidget.SetPlaceFactor(1.05) # Increases the size of the widget bounds
# b1,b2,b3 = vtkObj1.GetBounds()[::2]
# planeWidget.SetOrigin(b1,b2,b3)
planeWidget.SetOutsideBounds(0) # Not allow the widget to move outside the input bounds
planeWidget.SetScaleEnabled(0) # Ability to scale with the mouse
planeWidget.SetEnabled(1) # Starts the widget
planeWidget.SetOutlineTranslation(0) # Abiltiy to move the widget with the mouse
planeWidget.GetPlaneProperty().SetOpacity(0.1)
planeWidget.PlaceWidget()
planeWidget.AddObserver("InteractionEvent",movePlane)
# Orientation widget
oriWid = pvtkUtil.addDirectionWidget(iren,ren,150,35)
# Set renderer options
ren.SetBackground(1.0,1.0,1.0)
ren.AddActor(actor1)
ren.AddActor2D(scalarBarRes)
ren.AddViewProp(axes)
# Fix the colorbar title
title = scalarBarRes.GetTitle()
scalarBarRes.SetTitle('')
# Add the title at the top of the figure manually.
titText = vtk.vtkTextActor()
titText.GetPositionCoordinate().SetCoordinateSystemToNormalizedViewport()
titText.GetTextProperty().SetFontSize(35)
titText.GetTextProperty().SetColor(0.0,0.0,0.0)
titText.SetPosition(0.02,0.87)
titText.SetInput(title)
ren.AddActor(titText)
# Start the render Window
renwin.Render()
iren.Start()
# Save the fig
planeWidget.SetEnabled(0)
w2i = vtk.vtkWindowToImageFilter()
w2i.SetMagnification(1)
w2i.SetInput(renwin)
w2i.Update()
writer = vtk.vtkTIFFWriter()
writer.SetCompressionToNoCompression()
writer.SetInputConnection(w2i.GetOutputPort())
writer.SetFileName(figName + '.tif')
writer.Write()
if True:
camera = ren.GetActiveCamera()
# For playing around with the locations of the figures
# For printing out view information.
print('camera.GetClippingRange' + str(camera.GetClippingRange()))
print('camera.GetFocalPoint' + str(camera.GetFocalPoint()))
print('camera.GetPosition' + str(camera.GetPosition()))
print('camera.GetViewUp' + str(camera.GetViewUp()))
print('camera.GetParallelScale(' + str(camera.GetParallelScale()) +')')
# Close the window when exited
iren.TerminateApp()
renwin.Finalize()
del iren, renwin
# Gargage collect
gc.collect()
In [ ]:
In [22]:
# Load the model and mesh
In [23]:
mesh, modDict = simpeg.Mesh.TensorMesh.readVTK('MTwork/inv3d_HPK1/run_thibaut4_off/recoveredMod_run_thibaut4_off_it10.vtr')
In [24]:
# Plot the section
plot3dSetion(mesh,modDict,'GeologicalModel')
In [26]:
# Plot the background model
mesh, modDict = simpeg.Mesh.TensorMesh.readVTK('MTwork/inv3d_HPK1/run_thibaut4_off/nsmesh_CoarseHKPK1_NoExtension.vtr')
plot3dSetion(mesh,modDict,'MTwork/OriginalModel')
In [91]:
import simpegViz
In [92]:
modView = simpegViz.vtkView(mesh,{'C':modDict})
In [93]:
modView.limits = 1e-5, 0.01
modView.range = 1e-5, 0.01
modView.extent = [8,28,8,33,8,55]
In [94]:
modView.Show()
In [ ]: