From 6a101483a914e03c6c6ef34019ecb0be220ce1c9 Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Mon, 25 Nov 2013 17:22:46 -0800 Subject: [PATCH 01/16] Drasl --- SimPEG/InteractiveView.py | 335 ++++++++++++++++++++++ SimPEG/visulize/__init__.py | 2 + SimPEG/visulize/vtk/__init__.py | 2 + SimPEG/visulize/vtk/vtkTools.py | 361 ++++++++++++++++++++++++ SimPEG/visulize/vtk/vtkView.py | 137 +++++++++ notebooks/3DRenderingWithvtkTools.ipynb | 59 ++++ 6 files changed, 896 insertions(+) create mode 100644 SimPEG/InteractiveView.py create mode 100644 SimPEG/visulize/__init__.py create mode 100644 SimPEG/visulize/vtk/__init__.py create mode 100644 SimPEG/visulize/vtk/vtkTools.py create mode 100644 SimPEG/visulize/vtk/vtkView.py create mode 100644 notebooks/3DRenderingWithvtkTools.ipynb diff --git a/SimPEG/InteractiveView.py b/SimPEG/InteractiveView.py new file mode 100644 index 00000000..ba4f1b63 --- /dev/null +++ b/SimPEG/InteractiveView.py @@ -0,0 +1,335 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +from mpl_toolkits.mplot3d import Axes3D +from utils import mkvc + + +class TensorView(object): + """ + Provides viewing functions for TensorMesh + + This class is inherited by TensorMesh + """ + def __init__(self): + pass + + def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): + """ + Mesh.plotImage(I) + + Plots scalar fields on the given mesh. + + Input: + + :param numpy.array I: scalar field + + Optional Input: + + :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz + :param int figNum: number of figure to plot to + :param matplotlib.axes.Axes ax: axis to plot to + :param str direction: slice dimensions, 3D only ('x', 'y', 'z') + :param bool numbering: show numbering of slices, 3D only + :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' + :param bool showIt: call plt.show() + + .. plot:: examples/mesh/plot_image_2D.py + :include-source: + + .. plot:: examples/mesh/plot_image_3D.py + :include-source: + """ + assert type(I) == np.ndarray, "I must be a numpy array" + assert type(numbering) == bool, "numbering must be a bool" + assert direction in ["x", "y","z"], "direction must be either x,y, or z" + + + if imageType == 'CC': + assert I.size == self.nC, "Incorrect dimensions for CC." + elif imageType == 'N': + assert I.size == self.nN, "Incorrect dimensions for N." + elif imageType == 'Fx': + if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') + elif imageType == 'Fy': + if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') + elif imageType == 'Fz': + if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') + elif imageType == 'Ex': + if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') + elif imageType == 'Ey': + if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') + elif imageType == 'Ez': + if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') + elif imageType[0] == 'E': + plotAll = len(imageType) == 1 + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + fig = plt.figure(figNum) + # Determine the subplot number: 131, 121 + numPlots = 130 if plotAll else len(imageType)/2*10+100 + pltNum = 1 + ex, ey, ez = self.r(I,'E','E','M') + if plotAll or 'Ex' in imageType: + ax_x = plt.subplot(numPlots+pltNum) + self.plotImage(ex, imageType='Ex', ax=ax_x, **options) + pltNum +=1 + if plotAll or 'Ey' in imageType: + ax_y = plt.subplot(numPlots+pltNum) + self.plotImage(ey, imageType='Ey', ax=ax_y, **options) + pltNum +=1 + if plotAll or 'Ez' in imageType: + ax_z = plt.subplot(numPlots+pltNum) + self.plotImage(ez, imageType='Ez', ax=ax_z, **options) + pltNum +=1 + return + elif imageType[0] == 'F': + plotAll = len(imageType) == 1 + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + fig = plt.figure(figNum) + # Determine the subplot number: 131, 121 + numPlots = 130 if plotAll else len(imageType)/2*10+100 + pltNum = 1 + fx, fy, fz = self.r(I,'F','F','M') + if plotAll or 'Fx' in imageType: + ax_x = plt.subplot(numPlots+pltNum) + self.plotImage(fx, imageType='Fx', ax=ax_x, **options) + pltNum +=1 + if plotAll or 'Fy' in imageType: + ax_y = plt.subplot(numPlots+pltNum) + self.plotImage(fy, imageType='Fy', ax=ax_y, **options) + pltNum +=1 + if plotAll or 'Fz' in imageType: + ax_z = plt.subplot(numPlots+pltNum) + self.plotImage(fz, imageType='Fz', ax=ax_z, **options) + pltNum +=1 + return + else: + raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") + + + if ax is None: + fig = plt.figure(figNum) + fig.clf() + ax = plt.subplot(111) + else: + assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" + fig = ax.figure + + if self.dim == 1: + if imageType == 'CC': + ph = ax.plot(self.vectorCCx, I, '-ro') + elif imageType == 'N': + ph = ax.plot(self.vectorNx, I, '-bs') + ax.set_xlabel("x") + ax.axis('tight') + elif self.dim == 2: + if imageType == 'CC': + C = I[:].reshape(self.n, order='F') + elif imageType == 'N': + C = I[:].reshape(self.n+1, order='F') + C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) + elif imageType == 'Fx': + C = I[:].reshape(self.nFx, order='F') + C = 0.5*(C[:-1, :] + C[1:, :] ) + elif imageType == 'Fy': + C = I[:].reshape(self.nFy, order='F') + C = 0.5*(C[:, :-1] + C[:, 1:] ) + elif imageType == 'Ex': + C = I[:].reshape(self.nEx, order='F') + C = 0.5*(C[:,:-1] + C[:,1:] ) + elif imageType == 'Ey': + C = I[:].reshape(self.nEy, order='F') + C = 0.5*(C[:-1,:] + C[1:,:] ) + + ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) + ax.axis('tight') + ax.set_xlabel("x") + ax.set_ylabel("y") + + elif self.dim == 3: + if direction == 'z': + + # get copy of image and average to cell-centres is necessary + if imageType == 'CC': + Ic = I[:].reshape(self.n, order='F') + elif imageType == 'N': + Ic = I[:].reshape(self.n+1, order='F') + Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) + elif imageType == 'Fx': + Ic = I[:].reshape(self.nFx, order='F') + Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) + elif imageType == 'Fy': + Ic = I[:].reshape(self.nFy, order='F') + Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) + elif imageType == 'Fz': + Ic = I[:].reshape(self.nFz, order='F') + Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) + elif imageType == 'Ex': + Ic = I[:].reshape(self.nEx, order='F') + Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) + elif imageType == 'Ey': + Ic = I[:].reshape(self.nEy, order='F') + Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) + elif imageType == 'Ez': + Ic = I[:].reshape(self.nEz, order='F') + Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) + + # determine number oE slices in x and y dimension + nX = np.ceil(np.sqrt(self.nCz)) + nY = np.ceil(self.nCz/nX) + + # allocate space for montage + nCx = self.nCx + nCy = self.nCy + + C = np.zeros((nX*nCx,nY*nCy)) + + for iy in range(int(nY)): + for ix in range(int(nX)): + iz = ix + iy*nX + if iz < self.nCz: + C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] + else: + C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan + + C = np.ma.masked_where(np.isnan(C), C) + xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] + yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] + # Plot the mesh + ph = ax.pcolormesh(xx, yy, C.T) + # Plot the lines + gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) + gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) + # Repeat and seperate with NaN + gxX = np.c_[gx, gx, gx+np.nan].ravel() + gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() + gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() + gyY = np.c_[gy, gy, gy+np.nan].ravel() + ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) + ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) + ax.axis('tight') + + if numbering: + pad = np.sum(self.hx)*0.04 + for iy in range(int(nY)): + for ix in range(int(nX)): + iz = ix + iy*nX + if iz < self.nCz: + ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, + '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') + + ax.set_title(imageType) + if showIt: plt.show() + return ph + + def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. + + :param bool nodes: plot nodes + :param bool faces: plot faces + :param bool centers: plot centers + :param bool edges: plot edges + :param bool lines: plot lines connecting nodes + :param bool showIt: call plt.show() + + .. plot:: examples/mesh/plot_grid_2D.py + :include-source: + + .. plot:: examples/mesh/plot_grid_3D.py + :include-source: + """ + if self.dim == 1: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + ax.hold(True) + ax.plot(xn, np.ones(np.shape(xn)), 'bs') + ax.plot(xc, np.ones(np.shape(xc)), 'ro') + ax.plot(xn, np.ones(np.shape(xn)), 'k--') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + if showIt: plt.show() + elif self.dim == 2: + fig = plt.figure(2) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + xs1 = self.gridFx + xs2 = self.gridFy + + ax.hold(True) + if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') + if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') + if faces: + ax.plot(xs1[:, 0], xs1[:, 1], 'g>') + ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + + # Plot the grid lines + if lines: + NN = self.r(self.gridN, 'N', 'N', 'M') + X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() + X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] + plt.plot(X, Y) + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + if showIt: plt.show() + elif self.dim == 3: + fig = plt.figure(3) + fig.clf() + ax = fig.add_subplot(111, projection='3d') + xn = self.gridN + xc = self.gridCC + xfs1 = self.gridFx + xfs2 = self.gridFy + xfs3 = self.gridFz + + xes1 = self.gridEx + xes2 = self.gridEy + xes3 = self.gridEz + + ax.hold(True) + if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) + if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) + if faces: + ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) + ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) + ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) + if edges: + ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) + ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) + ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) + + # Plot the grid lines + if lines: + NN = self.r(self.gridN, 'N', 'N', 'M') + X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() + Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() + X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() + Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() + X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() + Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() + Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() + X = np.r_[X1, X2, X3] + Y = np.r_[Y1, Y2, Y3] + Z = np.r_[Z1, Z2, Z3] + plt.plot(X, Y, 'b-', zs=Z) + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + ax.set_zlabel('x3') + if showIt: plt.show() diff --git a/SimPEG/visulize/__init__.py b/SimPEG/visulize/__init__.py new file mode 100644 index 00000000..485f9782 --- /dev/null +++ b/SimPEG/visulize/__init__.py @@ -0,0 +1,2 @@ +import vtk +#import mpl diff --git a/SimPEG/visulize/vtk/__init__.py b/SimPEG/visulize/vtk/__init__.py new file mode 100644 index 00000000..6d60ed5c --- /dev/null +++ b/SimPEG/visulize/vtk/__init__.py @@ -0,0 +1,2 @@ +from vtkTools import vtkTools +from vtkView import vtkView \ No newline at end of file diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py new file mode 100644 index 00000000..50387a14 --- /dev/null +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -0,0 +1,361 @@ +import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb +from SimPEG.utils import mkvc + + +class vtkTools(object): + """ + Class that interacts with VTK visulization toolkit. + + """ + + def __init__(self): + """ Initializes the VTK vtkTools. + + """ + + pass + + @staticmethod + def makeCellVTKObject(mesh,model): + """ + Make and return a cell based VTK object for a simpeg mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells + + Output: + :rtype: vtkRecilinearGrid object + :return: vtkObj + """ + + # 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. + # Asaign the spatial information. + vtkObj = vtk.vtkRectilinearGrid() + vtkObj.SetDimensions(xD,yD,zD) + vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) + vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) + vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1)) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + return vtkObj + + @staticmethod + def makeFaceVTKObject(mesh,model): + """ + Make and return a face based VTK object for a simpeg mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). + Property array must be order hstack(Fx,Fy,Fz) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def faceR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTFx = np.prod(mesh.nFx) + FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) + FyCellBlock = np.array([],dtype='int64') + FzCellBlock = np.array([],dtype='int64') + # Second direction + if mesh.dim >= 2: + nTFy = np.prod(mesh.nFy) + FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) + # Third direction + if mesh.dim == 3: + nTFz = np.prod(mesh.nFz) + FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) + # Cells -cell array + FCellArr = vtk.vtkCellArray() + FCellArr.SetNumberOfCells(np.sum(mesh.nF)) + FCellArr.SetCells(np.sum(mesh.nF)*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + # Cell type + FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) + # Cell location + FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF)*5,5,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(FCellType,FCellLoc,FCellArr) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() + return vtkObj + + @staticmethod + def makeEdgeVTKObject(mesh,model): + """ + Make and return a edge based VTK object for a simpeg mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). + Property array must be order hstack(Ex,Ey,Ez) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def edgeR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTEx = np.prod(mesh.nEx) + ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)]) + # Second direction + if mesh.dim >= 2: + nTEy = np.prod(mesh.nEy) + EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)]) + # Third direction + if mesh.dim == 3: + nTEz = np.prod(mesh.nEz) + EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) + # Cells -cell array + ECellArr = vtk.vtkCellArray() + ECellArr.SetNumberOfCells(np.sum(mesh.nE)) + ECellArr.SetCells(np.sum(mesh.nE)*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) + # Cell type + ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(np.sum(mesh.nE),dtype='uint8'),deep=1) + # Cell location + ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nE)*3,3,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(ECellType,ECellLoc,ECellArr) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + return vtkObj + + @staticmethod + def makeRenderWindow(ren): + renwin = vtk.vtkRenderWindow() + renwin.AddRenderer(ren) + iren = vtk.vtkRenderWindowInteractor() + iren.SetRenderWindow(renwin) + + return iren, renwin + + + @staticmethod + def closeRenderWindow(iren): + renwin = iren.GetRenderWindow() + renwin.Finalize() + iren.TerminateApp() + + del iren, renwin + + @staticmethod + def makeVTKActor(vtkObj): + """ Makes a vtk mapper and Actor""" + mapper = vtk.vtkDataSetMapper() + mapper.SetInput(vtkObj) + actor = vtk.vtkActor() + actor.SetMapper(mapper) + actor.GetProperty().SetColor(0,0,0) + actor.GetProperty().SetRepresentationToWireframe() + return actor + + @staticmethod + def makeVTKLODActor(vtkObj,clipper): + """Make LOD vtk Actor""" + selectMapper = vtk.vtkDataSetMapper() + selectMapper.SetInputConnection(clipper.GetOutputPort()) + selectMapper.SetScalarVisibility(1) + selectMapper.SetColorModeToMapScalars() + selectMapper.SetScalarModeToUseCellData() + selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange()) + + selectActor = vtk.vtkLODActor() + selectActor.SetMapper(selectMapper) + selectActor.GetProperty().SetEdgeColor(1,0.5,0) + selectActor.GetProperty().SetEdgeVisibility(0) + selectActor.VisibilityOn() + selectActor.SetScale(1.01, 1.01, 1.01) + return selectActor + + @staticmethod + def setScalar2View(vtkObj,scalarName): + """ Sets the sclar to view """ + useArr = vtkObj.GetCellData().GetArray(scalarName) + if useArr == None: + raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName)) + vtkObj.GetCellData().SetActiveScalars(scalarName) + + @staticmethod + def makeRectiVTKVOIThres(vtkObj): + """Make volume of interest and threshold for rectilinear grid.""" + cellCore = vtk.vtkExtractRectilinearGrid() + cellCore.SetInput(vtkObj) + cellCore.SetVOI(vtkObj.GetExtent()) + cellThres = vtk.vtkThreshold() + cellThres.AllScalarsOn() + cellThres.SetInputConnection(cellCore.GetOutputPort()) + cellThres.ThresholdByUpper(-1) + cellThres.Update() + return cellThres.GetOutput(), cellCore.GetOutput() + + @staticmethod + def makePlaneClipper(vtkObj): + """Makes a plane and clipper """ + plane = vtk.vtkPlane() + clipper = vtk.vtkClipDataSet() + clipper.SetInputConnection(vtkObj.GetProducerPort()) + clipper.SetClipFunction(plane) + clipper.InsideOutOff() + return clipper, plane + + @staticmethod + def makePlaneWidget(vtkObj,iren,plane,actor): + """Make an interactive planeWidget""" + + # Callback function + def movePlane(obj, events): + obj.GetPlane(intPlane) + intActor.VisibilityOn() + + # Associate the line widget with the interactor + planeWidget = vtk.vtkImplicitPlaneWidget() + planeWidget.SetInteractor(iren) + planeWidget.SetPlaceFactor(1.25) + planeWidget.SetInput(vtkObj) + planeWidget.PlaceWidget() + #planeWidget.AddObserver("InteractionEvent", movePlane) + planeWidget.SetScaleEnabled(0) + planeWidget.SetEnabled(1) + planeWidget.SetOutlineTranslation(0) + planeWidget.GetPlaneProperty().SetOpacity(0.1) + return planeWidget + + + @staticmethod + def startRenderWindow(iren): + """ Start a vtk rendering window""" + iren.Initialize() + renwin = iren.GetRenderWindow() + renwin.Render() + iren.Start() + + + # Simple write/read VTK xml model functions. + @staticmethod + def writeVTPFile(fileName,vtkPolyObject): + '''Function to write vtk polydata file (vtp).''' + polyWriter = vtk.vtkXMLPolyDataWriter() + polyWriter.SetInput(vtkPolyObject) + polyWriter.SetFileName(fileName) + polyWriter.Update() + + @staticmethod + def writeVTUFile(fileName,vtkUnstructuredGrid): + '''Function to write vtk unstructured grid (vtu).''' + Writer = vtk.vtkXMLUnstructuredGridWriter() + Writer.SetInput(vtkUnstructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTRFile(fileName,vtkRectilinearGrid): + '''Function to write vtk rectilinear grid (vtr).''' + Writer = vtk.vtkXMLRectilinearGridWriter() + Writer.SetInput(vtkRectilinearGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTSFile(fileName,vtkStructuredGrid): + '''Function to write vtk structured grid (vts).''' + Writer = vtk.vtkXMLStructuredGridWriter() + Writer.SetInput(vtkStructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def readVTSFile(fileName): + '''Function to read vtk structured grid (vts) and return a grid object.''' + Reader = vtk.vtkXMLStructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTUFile(fileName): + '''Function to read vtk structured grid (vtu) and return a grid object.''' + Reader = vtk.vtkXMLUnstructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTRFile(fileName): + '''Function to read vtk structured grid (vtr) and return a grid object.''' + Reader = vtk.vtkXMLRectilinearGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTPFile(fileName): + '''Function to read vtk structured grid (vtp) and return a grid object.''' + Reader = vtk.vtkXMLPolyDataReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py new file mode 100644 index 00000000..334561ba --- /dev/null +++ b/SimPEG/visulize/vtk/vtkView.py @@ -0,0 +1,137 @@ +import numpy as np, vtk +import SimPEG as simpeg +#import SimPEG.visulize.vtk.vtkTools as vtkSP # Always get an error for this import + +class vtkView(object): + """ + Class for storing and view of SimPEG models in VTK (visulization toolkit). + + Inputs: + :param mesh, SimPEG mesh. + :param propdict, dictionary of property models. + Can have these dictionary names: + 'cell' - cell model; 'face' - face model; 'edge' - edge model + The dictionary properties are given as dictionaries with: + {'NameOfThePropertyModel': np.array of the properties}. + The property array has to be ordered in compliance with SimPEG standards. + + :: + Example of usages. + + ToDo + + """ + + def __init__(self,mesh,propdict): + """ + """ + + self.name = 'VTK figure of SimPEG model' + self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] + self.limits = [0, 10000] + self._mesh = mesh + + # Set vtk object containers + self._cell = None + self._faces = None + self._edges = None + + self._readPropertyDictionary(propdict) + + # Setup hidden properties + self._ren = None + self._iren = None + self._renwin = None + self._core = None + self._viewobj = None + self._plane = None + self._clipper = None + self._widget = None + self._actor = None + self._lut = None + + def _readPropertyDictionary(self,propdict): + """ + Reads the property and assigns to the object + """ + import SimPEG.visulize.vtk.vtkTools as vtkSP + + # Test the property dictionary + if len(propdict) > 3: + raise(Exception,'Too many input items in the property dictionary') + for propitem in propdict.iteritems(): + if propitem[0] in ['cell','face','edge']: + if propitem[0] == 'cell': + self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) + if propitem[0] == 'face': + self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) + if propitem[0] == 'edge': + self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) + else: + raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) + + def Show(self,imageType='cell'): + """ + Open the VTK figure window and show the mesh. + + Inputs: + param: str imageType: type of image {'cell','face','edge'} + + """ + #vtkSP = simpeg.visulize.vtk.vtkTools + import SimPEG.visulize.vtk.vtkTools as vtkSP + + # Make a renderer + self._ren = vtk.vtkRenderer() + # Make renderwindow. Returns the interactor. + self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) + + + # Sort out the actor + if imageType == 'cell': + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell) + elif imageType == 'face': + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._face) + elif imageType == 'edge': + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._edge) + else: + raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) + + + global intPlane, intActor + self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) + intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper) + self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor) + # Callback function + self._plane = intPlane + self._actor = intActor + def movePlane(obj, events): + global intPlane, intActor + obj.GetPlane(intPlane) + intActor.VisibilityOn() + + self._widget.AddObserver("InteractionEvent",movePlane) + lut = vtk.vtkLookupTable() + lut.SetNumberOfColors(256) + lut.SetHueRange(0,0.66667) + lut.Build() + self._lut = lut + self._actor.GetMapper().SetLookupTable(lut) + + # Set renderer options + self._ren.SetBackground(.5,.5,.5) + self._ren.AddActor(self._actor) + + # Start the render Window + vtkSP.startRenderWindow(self._iren) + # Close the window when exited + vtkSP.closeRenderWindow(self._iren) + del self._iren, self._renwin + + + + + + + + diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb new file mode 100644 index 00000000..7209b10e --- /dev/null +++ b/notebooks/3DRenderingWithvtkTools.ipynb @@ -0,0 +1,59 @@ +{ + "metadata": { + "name": "3D rendering with vtkTools" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": "import numpy as np, vtk\nimport SimPEG as simpeg", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n# Make the vtk viewer object.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + } + ], + "metadata": {} + } + ] +} From b29cc6e08d24aa3aea0bac455371e8436cfcc590 Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Mon, 25 Nov 2013 17:44:26 -0800 Subject: [PATCH 02/16] Deleted local 'visulize' folder --- SimPEG/visulize/__init__.py | 2 - SimPEG/visulize/vtk/__init__.py | 2 - SimPEG/visulize/vtk/vtkTools.py | 361 ------------------- SimPEG/visulize/vtk/vtkView.py | 137 ------- notebooks/3DRenderingWithvtkTools.ipynb.orig | 115 ++++++ 5 files changed, 115 insertions(+), 502 deletions(-) delete mode 100644 SimPEG/visulize/__init__.py delete mode 100644 SimPEG/visulize/vtk/__init__.py delete mode 100644 SimPEG/visulize/vtk/vtkTools.py delete mode 100644 SimPEG/visulize/vtk/vtkView.py create mode 100644 notebooks/3DRenderingWithvtkTools.ipynb.orig diff --git a/SimPEG/visulize/__init__.py b/SimPEG/visulize/__init__.py deleted file mode 100644 index 485f9782..00000000 --- a/SimPEG/visulize/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -import vtk -#import mpl diff --git a/SimPEG/visulize/vtk/__init__.py b/SimPEG/visulize/vtk/__init__.py deleted file mode 100644 index 6d60ed5c..00000000 --- a/SimPEG/visulize/vtk/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -from vtkTools import vtkTools -from vtkView import vtkView \ No newline at end of file diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py deleted file mode 100644 index 50387a14..00000000 --- a/SimPEG/visulize/vtk/vtkTools.py +++ /dev/null @@ -1,361 +0,0 @@ -import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb -from SimPEG.utils import mkvc - - -class vtkTools(object): - """ - Class that interacts with VTK visulization toolkit. - - """ - - def __init__(self): - """ Initializes the VTK vtkTools. - - """ - - pass - - @staticmethod - def makeCellVTKObject(mesh,model): - """ - Make and return a cell based VTK object for a simpeg mesh and model. - - Input: - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells - - Output: - :rtype: vtkRecilinearGrid object - :return: vtkObj - """ - - # 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. - # Asaign the spatial information. - vtkObj = vtk.vtkRectilinearGrid() - vtkObj.SetDimensions(xD,yD,zD) - vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) - vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) - vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1)) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - return vtkObj - - @staticmethod - def makeFaceVTKObject(mesh,model): - """ - Make and return a face based VTK object for a simpeg mesh and model. - - Input: - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). - Property array must be order hstack(Fx,Fy,Fz) - - Output: - :rtype: vtkUnstructuredGrid object - :return: vtkObj - """ - - ## Convert simpeg mesh to VTK properties - # Convert mesh nodes to vtkPoints - vtkPts = vtk.vtkPoints() - vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - - # Define the face "cells" - # Using VTK_QUAD cell for faces (see VTK file format) - nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') - def faceR(mat,length): - return mat.T.reshape((length,1)) - # First direction - nTFx = np.prod(mesh.nFx) - FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) - FyCellBlock = np.array([],dtype='int64') - FzCellBlock = np.array([],dtype='int64') - # Second direction - if mesh.dim >= 2: - nTFy = np.prod(mesh.nFy) - FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) - # Third direction - if mesh.dim == 3: - nTFz = np.prod(mesh.nFz) - FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) - # Cells -cell array - FCellArr = vtk.vtkCellArray() - FCellArr.SetNumberOfCells(np.sum(mesh.nF)) - FCellArr.SetCells(np.sum(mesh.nF)*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) - # Cell type - FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) - # Cell location - FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF)*5,5,dtype='int64'),deep=1) - - ## Make the object - vtkObj = vtk.vtkUnstructuredGrid() - # Set the objects properties - vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(FCellType,FCellLoc,FCellArr) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - vtkObj.Update() - return vtkObj - - @staticmethod - def makeEdgeVTKObject(mesh,model): - """ - Make and return a edge based VTK object for a simpeg mesh and model. - - Input: - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). - Property array must be order hstack(Ex,Ey,Ez) - - Output: - :rtype: vtkUnstructuredGrid object - :return: vtkObj - """ - - ## Convert simpeg mesh to VTK properties - # Convert mesh nodes to vtkPoints - vtkPts = vtk.vtkPoints() - vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - - # Define the face "cells" - # Using VTK_QUAD cell for faces (see VTK file format) - nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') - def edgeR(mat,length): - return mat.T.reshape((length,1)) - # First direction - nTEx = np.prod(mesh.nEx) - ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)]) - # Second direction - if mesh.dim >= 2: - nTEy = np.prod(mesh.nEy) - EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)]) - # Third direction - if mesh.dim == 3: - nTEz = np.prod(mesh.nEz) - EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) - # Cells -cell array - ECellArr = vtk.vtkCellArray() - ECellArr.SetNumberOfCells(np.sum(mesh.nE)) - ECellArr.SetCells(np.sum(mesh.nE)*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) - # Cell type - ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(np.sum(mesh.nE),dtype='uint8'),deep=1) - # Cell location - ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nE)*3,3,dtype='int64'),deep=1) - - ## Make the object - vtkObj = vtk.vtkUnstructuredGrid() - # Set the objects properties - vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(ECellType,ECellLoc,ECellArr) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - return vtkObj - - @staticmethod - def makeRenderWindow(ren): - renwin = vtk.vtkRenderWindow() - renwin.AddRenderer(ren) - iren = vtk.vtkRenderWindowInteractor() - iren.SetRenderWindow(renwin) - - return iren, renwin - - - @staticmethod - def closeRenderWindow(iren): - renwin = iren.GetRenderWindow() - renwin.Finalize() - iren.TerminateApp() - - del iren, renwin - - @staticmethod - def makeVTKActor(vtkObj): - """ Makes a vtk mapper and Actor""" - mapper = vtk.vtkDataSetMapper() - mapper.SetInput(vtkObj) - actor = vtk.vtkActor() - actor.SetMapper(mapper) - actor.GetProperty().SetColor(0,0,0) - actor.GetProperty().SetRepresentationToWireframe() - return actor - - @staticmethod - def makeVTKLODActor(vtkObj,clipper): - """Make LOD vtk Actor""" - selectMapper = vtk.vtkDataSetMapper() - selectMapper.SetInputConnection(clipper.GetOutputPort()) - selectMapper.SetScalarVisibility(1) - selectMapper.SetColorModeToMapScalars() - selectMapper.SetScalarModeToUseCellData() - selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange()) - - selectActor = vtk.vtkLODActor() - selectActor.SetMapper(selectMapper) - selectActor.GetProperty().SetEdgeColor(1,0.5,0) - selectActor.GetProperty().SetEdgeVisibility(0) - selectActor.VisibilityOn() - selectActor.SetScale(1.01, 1.01, 1.01) - return selectActor - - @staticmethod - def setScalar2View(vtkObj,scalarName): - """ Sets the sclar to view """ - useArr = vtkObj.GetCellData().GetArray(scalarName) - if useArr == None: - raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName)) - vtkObj.GetCellData().SetActiveScalars(scalarName) - - @staticmethod - def makeRectiVTKVOIThres(vtkObj): - """Make volume of interest and threshold for rectilinear grid.""" - cellCore = vtk.vtkExtractRectilinearGrid() - cellCore.SetInput(vtkObj) - cellCore.SetVOI(vtkObj.GetExtent()) - cellThres = vtk.vtkThreshold() - cellThres.AllScalarsOn() - cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(-1) - cellThres.Update() - return cellThres.GetOutput(), cellCore.GetOutput() - - @staticmethod - def makePlaneClipper(vtkObj): - """Makes a plane and clipper """ - plane = vtk.vtkPlane() - clipper = vtk.vtkClipDataSet() - clipper.SetInputConnection(vtkObj.GetProducerPort()) - clipper.SetClipFunction(plane) - clipper.InsideOutOff() - return clipper, plane - - @staticmethod - def makePlaneWidget(vtkObj,iren,plane,actor): - """Make an interactive planeWidget""" - - # Callback function - def movePlane(obj, events): - obj.GetPlane(intPlane) - intActor.VisibilityOn() - - # Associate the line widget with the interactor - planeWidget = vtk.vtkImplicitPlaneWidget() - planeWidget.SetInteractor(iren) - planeWidget.SetPlaceFactor(1.25) - planeWidget.SetInput(vtkObj) - planeWidget.PlaceWidget() - #planeWidget.AddObserver("InteractionEvent", movePlane) - planeWidget.SetScaleEnabled(0) - planeWidget.SetEnabled(1) - planeWidget.SetOutlineTranslation(0) - planeWidget.GetPlaneProperty().SetOpacity(0.1) - return planeWidget - - - @staticmethod - def startRenderWindow(iren): - """ Start a vtk rendering window""" - iren.Initialize() - renwin = iren.GetRenderWindow() - renwin.Render() - iren.Start() - - - # Simple write/read VTK xml model functions. - @staticmethod - def writeVTPFile(fileName,vtkPolyObject): - '''Function to write vtk polydata file (vtp).''' - polyWriter = vtk.vtkXMLPolyDataWriter() - polyWriter.SetInput(vtkPolyObject) - polyWriter.SetFileName(fileName) - polyWriter.Update() - - @staticmethod - def writeVTUFile(fileName,vtkUnstructuredGrid): - '''Function to write vtk unstructured grid (vtu).''' - Writer = vtk.vtkXMLUnstructuredGridWriter() - Writer.SetInput(vtkUnstructuredGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def writeVTRFile(fileName,vtkRectilinearGrid): - '''Function to write vtk rectilinear grid (vtr).''' - Writer = vtk.vtkXMLRectilinearGridWriter() - Writer.SetInput(vtkRectilinearGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def writeVTSFile(fileName,vtkStructuredGrid): - '''Function to write vtk structured grid (vts).''' - Writer = vtk.vtkXMLStructuredGridWriter() - Writer.SetInput(vtkStructuredGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def readVTSFile(fileName): - '''Function to read vtk structured grid (vts) and return a grid object.''' - Reader = vtk.vtkXMLStructuredGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTUFile(fileName): - '''Function to read vtk structured grid (vtu) and return a grid object.''' - Reader = vtk.vtkXMLUnstructuredGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTRFile(fileName): - '''Function to read vtk structured grid (vtr) and return a grid object.''' - Reader = vtk.vtkXMLRectilinearGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTPFile(fileName): - '''Function to read vtk structured grid (vtp) and return a grid object.''' - Reader = vtk.vtkXMLPolyDataReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py deleted file mode 100644 index 334561ba..00000000 --- a/SimPEG/visulize/vtk/vtkView.py +++ /dev/null @@ -1,137 +0,0 @@ -import numpy as np, vtk -import SimPEG as simpeg -#import SimPEG.visulize.vtk.vtkTools as vtkSP # Always get an error for this import - -class vtkView(object): - """ - Class for storing and view of SimPEG models in VTK (visulization toolkit). - - Inputs: - :param mesh, SimPEG mesh. - :param propdict, dictionary of property models. - Can have these dictionary names: - 'cell' - cell model; 'face' - face model; 'edge' - edge model - The dictionary properties are given as dictionaries with: - {'NameOfThePropertyModel': np.array of the properties}. - The property array has to be ordered in compliance with SimPEG standards. - - :: - Example of usages. - - ToDo - - """ - - def __init__(self,mesh,propdict): - """ - """ - - self.name = 'VTK figure of SimPEG model' - self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] - self.limits = [0, 10000] - self._mesh = mesh - - # Set vtk object containers - self._cell = None - self._faces = None - self._edges = None - - self._readPropertyDictionary(propdict) - - # Setup hidden properties - self._ren = None - self._iren = None - self._renwin = None - self._core = None - self._viewobj = None - self._plane = None - self._clipper = None - self._widget = None - self._actor = None - self._lut = None - - def _readPropertyDictionary(self,propdict): - """ - Reads the property and assigns to the object - """ - import SimPEG.visulize.vtk.vtkTools as vtkSP - - # Test the property dictionary - if len(propdict) > 3: - raise(Exception,'Too many input items in the property dictionary') - for propitem in propdict.iteritems(): - if propitem[0] in ['cell','face','edge']: - if propitem[0] == 'cell': - self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'face': - self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'edge': - self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) - else: - raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) - - def Show(self,imageType='cell'): - """ - Open the VTK figure window and show the mesh. - - Inputs: - param: str imageType: type of image {'cell','face','edge'} - - """ - #vtkSP = simpeg.visulize.vtk.vtkTools - import SimPEG.visulize.vtk.vtkTools as vtkSP - - # Make a renderer - self._ren = vtk.vtkRenderer() - # Make renderwindow. Returns the interactor. - self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) - - - # Sort out the actor - if imageType == 'cell': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell) - elif imageType == 'face': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._face) - elif imageType == 'edge': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._edge) - else: - raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) - - - global intPlane, intActor - self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) - intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper) - self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor) - # Callback function - self._plane = intPlane - self._actor = intActor - def movePlane(obj, events): - global intPlane, intActor - obj.GetPlane(intPlane) - intActor.VisibilityOn() - - self._widget.AddObserver("InteractionEvent",movePlane) - lut = vtk.vtkLookupTable() - lut.SetNumberOfColors(256) - lut.SetHueRange(0,0.66667) - lut.Build() - self._lut = lut - self._actor.GetMapper().SetLookupTable(lut) - - # Set renderer options - self._ren.SetBackground(.5,.5,.5) - self._ren.AddActor(self._actor) - - # Start the render Window - vtkSP.startRenderWindow(self._iren) - # Close the window when exited - vtkSP.closeRenderWindow(self._iren) - del self._iren, self._renwin - - - - - - - - diff --git a/notebooks/3DRenderingWithvtkTools.ipynb.orig b/notebooks/3DRenderingWithvtkTools.ipynb.orig new file mode 100644 index 00000000..8d2c833d --- /dev/null +++ b/notebooks/3DRenderingWithvtkTools.ipynb.orig @@ -0,0 +1,115 @@ +{ + "metadata": { +<<<<<<< HEAD + "name": "3D rendering with vtkTools" +======= + "name": "" +>>>>>>> develop + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "import numpy as np, vtk\nimport SimPEG as simpeg", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 +======= + "input": [ + "import numpy as np, vtk\n", + "import SimPEG as simpeg" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 +>>>>>>> develop + }, + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 +======= + "input": [ + "#Make a mesh and model\n", + "x0 = np.zeros(3)\n", + "h1 = np.ones(20)*50\n", + "h2 = np.ones(10)*100\n", + "h3 = np.ones(5)*200\n", + "\n", + "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 6 +>>>>>>> develop + }, + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n# Make the vtk viewer object.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 +======= + "input": [ + "# Make a models that correspond to the cells, faces and edges.\n", + "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n", + "# Make the vtk viewer object.\n", + "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) \n", + " " + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 7 +>>>>>>> develop + }, + { + "cell_type": "code", + "collapsed": false, +<<<<<<< HEAD + "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 +======= + "input": [ + "# Show the image \n", + "vtkViewer.Show()\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": "*" +>>>>>>> develop + } + ], + "metadata": {} + } + ] +} From 4b292b574d68b2d4a63dd477b500743de95e9feb Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Tue, 26 Nov 2013 07:45:18 -0800 Subject: [PATCH 03/16] Updating vtkView --- SimPEG/visualize/vtk/vtkView.py | 60 ++++++++++++++++++++++++--------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index abedd0bf..f65c8e7b 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -14,7 +14,7 @@ class vtkView(object): :param mesh, SimPEG mesh. :param propdict, dictionary of property models. Can have these dictionary names: - 'cell' - cell model; 'face' - face model; 'edge' - edge model + 'C' - cell model; 'F' - face model; 'E' - edge model The dictionary properties are given as dictionaries with: {'NameOfThePropertyModel': np.array of the properties}. The property array has to be ordered in compliance with SimPEG standards. @@ -30,11 +30,16 @@ class vtkView(object): """ """ - # ToDo: Set the properties up so that there are set/get methods + # Error check inputs + if type(mesh).__name__ != 'TensorMesh': + raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh)) + + + # Set default values self.name = 'VTK figure of SimPEG model' - self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] - self.limits = [0, 1e12] - self.viewprop = {'cell':0} # Name of the tyep and Int order of the array or name of the vector. + self._extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] + self._limits = [0, 1e12] + self._viewprop = {'C':0} # Name of the tyep and Int order of the array or name of the vector. self._mesh = mesh @@ -57,6 +62,31 @@ class vtkView(object): self._actor = None self._lut = None + # Set/Get properties + @property + def extent(self): + return getattr(self,_extent) + + @extent.setter + def extent(self,value): + + # Error check + valnp = np.array(value) + if valnp.dtype != int or len(valnp) != 6: + raise Exception('.extent has to be list or nparray of 6 integers.') + + self._extent = valnp + + @property + def limits(self): + return getattr(self._limits) + + @limits.setter + def limits(self,value): + valnp = np.array(value) + if valnp.dtype != float or len(valnp) != 2: + raise Exception('.limits has to be list or nparray of 2 floats.') + def _readPropertyDictionary(self,propdict): """ Reads the property and assigns to the object @@ -67,15 +97,15 @@ class vtkView(object): if len(propdict) > 3: raise(Exception,'Too many input items in the property dictionary') for propitem in propdict.iteritems(): - if propitem[0] in ['cell','face','edge']: - if propitem[0] == 'cell': + if propitem[0] in ['C','F','E']: + if propitem[0] == 'C': self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'face': + if propitem[0] == 'F': self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) - if propitem[0] == 'edge': + if propitem[0] == 'E': self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) else: - raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) + raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) def Show(self): """ @@ -91,16 +121,16 @@ class vtkView(object): imageType = self.viewprop.keys()[0] # Sort out the actor - if imageType == 'cell': + if imageType == 'C': self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) - elif imageType == 'face': + elif imageType == 'F': extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._face,extent,self.limits) - elif imageType == 'edge': + elif imageType == 'E': extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edge,extent,self.limits) else: - raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) + raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType)) # Set the active scalar. if type(self.viewprop.values()[0]) == int: @@ -153,7 +183,7 @@ if __name__ == '__main__': mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0) # Make a models that correspond to the cells, faces and edges. - models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} + models = {'C':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'E':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} # Make the vtk viewer object. vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) # Show the image From 1ec3c5cd0936d0b014571fdee1aca78c291dc6c0 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Tue, 26 Nov 2013 21:07:30 -0800 Subject: [PATCH 04/16] Fixed vtkTools and vtkView to work for all cell,faces, edges. Changed index for cells, faces and edges to be C,F,E respectively. Made viewprop, limits, extent as @properties, with error checks in the set methods. Need to finish the limits error checks. --- SimPEG/visualize/vtk/vtkTools.py | 6 +- SimPEG/visualize/vtk/vtkView.py | 172 +++++++++++++++++++++++-------- 2 files changed, 134 insertions(+), 44 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index 9853d089..7b405d4c 100644 --- a/SimPEG/visualize/vtk/vtkTools.py +++ b/SimPEG/visualize/vtk/vtkTools.py @@ -61,6 +61,7 @@ class vtkTools(object): vtkObj.GetCellData().AddArray(vtkDoubleArr) vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() return vtkObj @staticmethod @@ -104,7 +105,7 @@ class vtkTools(object): # Cells -cell array FCellArr = vtk.vtkCellArray() FCellArr.SetNumberOfCells(mesh.nF) - FCellArr.SetCells(mesh.nF*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + FCellArr.SetCells(mesh.nF,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) # Cell type FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1) # Cell location @@ -166,7 +167,7 @@ class vtkTools(object): # Cells -cell array ECellArr = vtk.vtkCellArray() ECellArr.SetNumberOfCells(mesh.nE) - ECellArr.SetCells(mesh.nE*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) + ECellArr.SetCells(mesh.nE,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) # Cell type ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1) # Cell location @@ -186,6 +187,7 @@ class vtkTools(object): vtkObj.GetCellData().AddArray(vtkDoubleArr) vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() return vtkObj @staticmethod diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index f65c8e7b..dde60fd3 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -8,16 +8,16 @@ import SimPEG as simpeg class vtkView(object): """ - Class for storing and view of SimPEG models in VTK (visulization toolkit). + Class for storing and view of SimPEG models in VTK (visualization toolkit). Inputs: :param mesh, SimPEG mesh. :param propdict, dictionary of property models. Can have these dictionary names: - 'C' - cell model; 'F' - face model; 'E' - edge model - The dictionary properties are given as dictionaries with: + 'C' - cell model; 'F' - face model; 'E' - edge model; ('V' - vector field : NOT SUPPORTED) + The dictionary values are given as dictionaries with: {'NameOfThePropertyModel': np.array of the properties}. - The property array has to be ordered in compliance with SimPEG standards. + The property np.array has to be ordered in compliance with SimPEG standards. :: Example of usages. @@ -30,27 +30,7 @@ class vtkView(object): """ """ - # Error check inputs - if type(mesh).__name__ != 'TensorMesh': - raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh)) - - - # Set default values - self.name = 'VTK figure of SimPEG model' - self._extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] - self._limits = [0, 1e12] - self._viewprop = {'C':0} # Name of the tyep and Int order of the array or name of the vector. - self._mesh = mesh - - - # Set vtk object containers - self._cell = None - self._faces = None - self._edges = None - - self._readPropertyDictionary(propdict) - - # Setup hidden properties + # Setup hidden properties, used for the visualization self._ren = None self._iren = None self._renwin = None @@ -61,31 +41,122 @@ class vtkView(object): self._widget = None self._actor = None self._lut = None + # Set vtk object containers + self._cells = None + self._faces = None + self._edges = None + self._vectors = None # Not implemented + # Set default values + self.name = 'VTK figure of SimPEG model' + + + # Error check the input mesh + if type(mesh).__name__ != 'TensorMesh': + raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh)) + # Set the mesh + self._mesh = mesh + + # Read the property dictionary + self._readPropertyDictionary(propdict) + + + + # Set/Get properties @property def extent(self): - return getattr(self,_extent) + ''' Extent of the sub-domain of the model to view''' + if getattr(self,'_extent',None) is None: + self._extent = [0,self._mesh.nCx-1,0,self._mesh.nCy-1,0,self._mesh.nCz-1] + return self._extent @extent.setter def extent(self,value): + import warnings # Error check - valnp = np.array(value) + valnp = np.array(value,dtype=int) if valnp.dtype != int or len(valnp) != 6: raise Exception('.extent has to be list or nparray of 6 integers.') - + # Test the range of the values + loB = np.zeros(3,dtype=int) + upB = np.array(self._mesh.nCv - np.ones(3),dtype=int) + # Test the bounds + change = 0 + # Test for lower bounds, can't be smaller the 0 + #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) + tlb = valnp[::2] < loB + if tlb.any(): + valnp[::2][tlb] = loB[tlb] + change = 1 + warnings.warn('Lower bounds smaller then 0') + # Test for lower bounds, can't be larger then upB + #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) + tlub = valnp[::2] > upB + if tlub.any(): + valnp[::2][tlub] = upB[tlub] - 1 + change = 1 + warnings.warn('Lower bounds larger then uppermost bounds') + # Test for upper bounds, can't be larger the extent of the mesh + #tub = np.prod(valnp[1::2] > upB,axis=0,dtype=bool) + tub = valnp[1::2] > upB + if tub.any(): + valnp[1::2][tub] = upB[tub] + change = 1 + warnings.warn('Upper bounds greater then number of cells') + # Test if lower is smaller the upper + #tgt = np.prod(valnp[::2] < valnp[1::2],axis=0,dtype=bool) + tgt = valnp[::2] > valnp[1::2] + if tgt.any(): + valnp[1::2][tgt] = valnp[::2][tgt] + 1 + change = 1 + warnings.warn('Lower bounds greater the Upper bounds') + # Print a warning + if change: + warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist())) + + # Set extent self._extent = valnp @property def limits(self): - return getattr(self._limits) + ''' Lower and upper limits (cutoffs) of the values to view. ''' + return getattr(self,'_limits',None) + @limits.setter def limits(self,value): - valnp = np.array(value) - if valnp.dtype != float or len(valnp) != 2: - raise Exception('.limits has to be list or nparray of 2 floats.') + if value is None: + self._limits = None + else: + valnp = np.array(value) + if valnp.dtype != float or len(valnp) != 2: + raise Exception('.limits has to be list or numpy array of 2 floats.') + self._limits = valnp + + + @property + def viewprop(self): + ''' Controls the property that will be viewed.''' + + if getattr(self,'_viewprop',None) is None: + self._viewprop = {'C':0} # Name of the type and Int order of the array or name of the vector. + return self._viewprop + + @viewprop.setter + def viewprop(self,value): + if type(value) != dict: + raise Exception('{:s} has to be a python dictionary containing property type and name index. ') + if len(value) > 1: + raise Exception('Too many input items in the viewprop dictionary') + if value.keys()[0] not in ['C','F','E']: + raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) + if not(type(self.viewprop.values()[0]) is int or type(self.viewprop.values()[0]) is str): + raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be integer or a string with the index.') + + + self._viewprop = value def _readPropertyDictionary(self,propdict): """ @@ -94,18 +165,20 @@ class vtkView(object): import SimPEG.visualize.vtk.vtkTools as vtkSP # Test the property dictionary - if len(propdict) > 3: - raise(Exception,'Too many input items in the property dictionary') + if type(propdict) != dict: + raise Exception('{:s} has to be a python dictionary containing property models. ') + if len(propdict) > 4: + raise Exception('Too many input items in the property dictionary') for propitem in propdict.iteritems(): if propitem[0] in ['C','F','E']: if propitem[0] == 'C': - self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) + self._cells = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) if propitem[0] == 'F': - self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) + self._faces = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) if propitem[0] == 'E': - self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) + self._edges = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) else: - raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) + raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0])) def Show(self): """ @@ -120,15 +193,22 @@ class vtkView(object): self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) imageType = self.viewprop.keys()[0] + # Sort out the actor if imageType == 'C': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) + if self.limits is None: + self.limits = self._cells.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cells,self.extent,self.limits) elif imageType == 'F': + if self.limits is None: + self.limits = self._faces.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] - self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._face,extent,self.limits) + self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._faces,extent,self.limits) elif imageType == 'E': + if self.limits is None: + self.limits = self._edges.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] - self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edge,extent,self.limits) + self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edges,extent,self.limits) else: raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType)) @@ -174,6 +254,8 @@ class vtkView(object): if __name__ == '__main__': + + #Make a mesh and model x0 = np.zeros(3) h1 = np.ones(20)*50 @@ -183,8 +265,14 @@ if __name__ == '__main__': mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0) # Make a models that correspond to the cells, faces and edges. - models = {'C':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'E':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} + models = {'C':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,mesh.nF),'AllOnce':np.ones(mesh.nF)},'E':{'Test':np.arange(0,mesh.nE),'AllOnce':np.ones(mesh.nE)}} # Make the vtk viewer object. vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) + # Set the .viewprop for which model to view + vtkViewer.viewprop = {'F':'Test'} # Show the image vtkViewer.Show() + + # Set subset of the mesh to view (remove padding) + vtkViewer.extent = [4,14,0,7,0,3] + vtkViewer.Show() From 5fd4290927702e191b4900d09cf3c0e4dacc1759 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 09:11:02 -0800 Subject: [PATCH 05/16] Fixed errors. C,F,E works with limits, extent. Added .cmap for color mapping and .range to set the color range. Need to add colorbar, and more settings as it comes along. --- SimPEG/visualize/vtk/vtkTools.py | 8 ++-- SimPEG/visualize/vtk/vtkView.py | 78 ++++++++++++++++++++++---------- 2 files changed, 58 insertions(+), 28 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index 7b405d4c..c3d9f4a7 100644 --- a/SimPEG/visualize/vtk/vtkTools.py +++ b/SimPEG/visualize/vtk/vtkTools.py @@ -46,7 +46,7 @@ class vtkTools(object): vZ = mesh.vectorNz zD = mesh.nNz # Use rectilinear VTK grid. - # Asaign the spatial information. + # Assign the spatial information. vtkObj = vtk.vtkRectilinearGrid() vtkObj.SetDimensions(xD,yD,zD) vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) @@ -257,8 +257,7 @@ class vtkTools(object): cellThres = vtk.vtkThreshold() cellThres.AllScalarsOn() cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(limits[0]) - cellThres.ThresholdByLower(limits[1]) + cellThres.ThresholdBetween(limits[0],limits[1]) cellThres.Update() return cellThres.GetOutput(), cellCore.GetOutput() @@ -273,8 +272,7 @@ class vtkTools(object): cellThres = vtk.vtkThreshold() cellThres.AllScalarsOn() cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(limits[0]) - cellThres.ThresholdByLower(limits[1]) + cellThres.ThresholdBetween(limits[0],limits[1]) cellThres.Update() return cellThres.GetOutput(), cellCore.GetOutput() diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index dde60fd3..666cde57 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -1,6 +1,6 @@ -import numpy as np +import numpy as np, matplotlib as mpl try: - import vtk + import vtk, vtk.util.numpy_support as npsup #import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import except Exception, e: print 'VTK import error. Please ensure you have VTK installed to use this visualization package.' @@ -64,13 +64,37 @@ class vtkView(object): # Set/Get properties + @property + def cmap(self): + ''' Colormap to use in vtkView. Colormap is a matplotlib cmap(cm) array, has to be uint8(use flag bytes=True during cmap generation).''' + if getattr(self,'_cmap',None) is None: + # Set default + self._cmap = mpl.cm.hsv(np.arange(0.,1.,0.05),bytes=True) + return self._cmap + @cmap.setter + def cmap(self,value): + if value.min() > 0 or value.max() < 255 or value.shape[1] != 4 or value.dtype != np.uint8: + raise Exception('Input not an allowed array.\n Use matplotlib.cm to generate an array of size [nrColors,4] and dtype = uint8(flag bytes=True).') + self._cmap = value + + @property + def range(self): + ''' Range of the colors in vtkView.''' + if getattr(self,'_range',None) is None: + self._range = np.array(self._getActiveVTKobj().GetArray(self.viewprop.values()[0]).GetRange()) + return self._range + @range.setter + def range(self,value): + if type(value) not in [tuple, list, np.ndarray] or len(value) != 2 or np.array(value).dtype is not np.dtype('float'): + raise Exception('Input not in correct format. \n Has to be a list, tuple or np.arry of 2 floats.') + self._range = np.array(value) + @property def extent(self): ''' Extent of the sub-domain of the model to view''' if getattr(self,'_extent',None) is None: self._extent = [0,self._mesh.nCx-1,0,self._mesh.nCy-1,0,self._mesh.nCz-1] return self._extent - @extent.setter def extent(self,value): @@ -85,28 +109,24 @@ class vtkView(object): # Test the bounds change = 0 # Test for lower bounds, can't be smaller the 0 - #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) tlb = valnp[::2] < loB if tlb.any(): valnp[::2][tlb] = loB[tlb] change = 1 warnings.warn('Lower bounds smaller then 0') # Test for lower bounds, can't be larger then upB - #tlb = np.prod(valnp[::2] < loB,axis=0,dtype=bool) tlub = valnp[::2] > upB if tlub.any(): valnp[::2][tlub] = upB[tlub] - 1 change = 1 warnings.warn('Lower bounds larger then uppermost bounds') # Test for upper bounds, can't be larger the extent of the mesh - #tub = np.prod(valnp[1::2] > upB,axis=0,dtype=bool) tub = valnp[1::2] > upB if tub.any(): valnp[1::2][tub] = upB[tub] change = 1 warnings.warn('Upper bounds greater then number of cells') # Test if lower is smaller the upper - #tgt = np.prod(valnp[::2] < valnp[1::2],axis=0,dtype=bool) tgt = valnp[::2] > valnp[1::2] if tgt.any(): valnp[1::2][tgt] = valnp[::2][tgt] + 1 @@ -123,8 +143,6 @@ class vtkView(object): def limits(self): ''' Lower and upper limits (cutoffs) of the values to view. ''' return getattr(self,'_limits',None) - - @limits.setter def limits(self,value): if value is None: @@ -143,7 +161,6 @@ class vtkView(object): if getattr(self,'_viewprop',None) is None: self._viewprop = {'C':0} # Name of the type and Int order of the array or name of the vector. return self._viewprop - @viewprop.setter def viewprop(self,value): if type(value) != dict: @@ -158,6 +175,20 @@ class vtkView(object): self._viewprop = value + def _getActiveVTKobj(self): + """ + Finds the active VTK object. + """ + + if self.viewprop.keys()[0] is 'C': + vtkCellData = self._cells.GetCellData() + elif self.viewprop.keys()[0] is 'F': + vtkCellData = self._faces.GetCellData() + elif self.viewprop.keys()[0] is 'E': + vtkCellData = self._edges.GetCellData() + + return vtkCellData + def _readPropertyDictionary(self,propdict): """ Reads the property and assigns to the object @@ -192,9 +223,17 @@ class vtkView(object): # Make renderwindow. Returns the interactor. self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) - imageType = self.viewprop.keys()[0] - + + # Set the active scalar. + if type(self.viewprop.values()[0]) == int: + actScalar = self._getActiveVTKobj().GetArrayName(self.viewprop.values()[0]) + elif type(self.viewprop.values()[0]) == str: + actScalar = self.viewprop.values()[0] + else : + raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be interger or a string.') + self._getActiveVTKobj().SetActiveScalars(actScalar) # Sort out the actor + imageType = self.viewprop.keys()[0] if imageType == 'C': if self.limits is None: self.limits = self._cells.GetCellData().GetArray(self.viewprop.values()[0]).GetRange() @@ -211,15 +250,7 @@ class vtkView(object): self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edges,extent,self.limits) else: raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType)) - - # Set the active scalar. - if type(self.viewprop.values()[0]) == int: - actScalar = self._vtkobj.GetCellData().GetArrayName(self.viewprop.values()[0]) - elif type(self.viewprop.values()[0]) == str: - actScalar = self.viewprop.values()[0] - else : - raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be interger or a string.') - self._vtkobj.GetCellData().SetActiveScalars(actScalar) + #self._vtkobj.GetCellData().SetActiveScalars(actScalar) # Set up the plane, clipper and the user interaction. global intPlane, intActor self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) @@ -235,10 +266,11 @@ class vtkView(object): self._widget.AddObserver("InteractionEvent",movePlane) lut = vtk.vtkLookupTable() - lut.SetNumberOfColors(256) - lut.SetHueRange(0,0.66667) + lut.SetNumberOfColors(len(self.cmap)) + lut.SetTable(npsup.numpy_to_vtk(self.cmap)) lut.Build() self._lut = lut + self._actor.GetMapper().SetScalarRange(self.range) self._actor.GetMapper().SetLookupTable(lut) # Set renderer options From 4cdafc12c1b5a29a568b3bfc433fb9b9c21db090 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 09:29:26 -0800 Subject: [PATCH 06/16] Added a colorbar to the figure. --- SimPEG/visualize/vtk/vtkView.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index 666cde57..8143f6e2 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -189,6 +189,17 @@ class vtkView(object): return vtkCellData + def _getActiveArrayName(self): + """ + Finds the name of the active array. + """ + actArr = self.viewprop.values()[0] + if type(actArr) is str: + activeName = actArr + elif type(actArr) is int: + activeName = self._getActiveVTKobj().GetArrayName(actArr) + return activeName + def _readPropertyDictionary(self,propdict): """ Reads the property and assigns to the object @@ -270,12 +281,23 @@ class vtkView(object): lut.SetTable(npsup.numpy_to_vtk(self.cmap)) lut.Build() self._lut = lut + scalarBar = vtk.vtkScalarBarActor() + scalarBar.SetLookupTable(lut) + scalarBar.SetTitle(self._getActiveArrayName()) + scalarBar.GetPositionCoordinate().SetCoordinateSystemToNormalizedViewport() + scalarBar.GetPositionCoordinate().SetValue(0.1,0.01) + scalarBar.SetOrientationToHorizontal() + scalarBar.SetWidth(0.8) + scalarBar.SetHeight(0.17) + self._actor.GetMapper().SetScalarRange(self.range) self._actor.GetMapper().SetLookupTable(lut) # Set renderer options self._ren.SetBackground(.5,.5,.5) self._ren.AddActor(self._actor) + self._ren.AddActor2D(scalarBar) + self._renwin.SetSize(450,450) # Start the render Window vtkSP.startRenderWindow(self._iren) From 9cca076651fe0487f93271ef66772cdd3598128e Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 09:52:25 -0800 Subject: [PATCH 07/16] Updated the example at the bottom of the vtkView. __name__ == __main__ --- SimPEG/visualize/vtk/vtkView.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index 8143f6e2..3b52b992 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -312,14 +312,18 @@ if __name__ == '__main__': #Make a mesh and model x0 = np.zeros(3) - h1 = np.ones(20)*50 - h2 = np.ones(10)*100 - h3 = np.ones(5)*200 + h1 = np.ones(60)*50 + h2 = np.ones(60)*100 + h3 = np.ones(50)*200 mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0) # Make a models that correspond to the cells, faces and edges. - models = {'C':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,mesh.nF),'AllOnce':np.ones(mesh.nF)},'E':{'Test':np.arange(0,mesh.nE),'AllOnce':np.ones(mesh.nE)}} + t = np.ones(mesh.nC) + t[10000:50000] = 100 + t[100000:120000] = 100 + t[100000:120000] = 50 + models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t, 'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,mesh.nF),'AllOnce':np.ones(mesh.nF)},'E':{'Test':np.arange(0,mesh.nE),'AllOnce':np.ones(mesh.nE)}} # Make the vtk viewer object. vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) # Set the .viewprop for which model to view @@ -330,3 +334,17 @@ if __name__ == '__main__': # Set subset of the mesh to view (remove padding) vtkViewer.extent = [4,14,0,7,0,3] vtkViewer.Show() + + # Change viewing property + vtkViewer.viewprop = {'C':'Model'} + # Set the color range + # Reset extent. + vtkViewer.extent = [-1,1000,-1,1000,-1,1000] + vtkViewer.range = [0.,100.] + vtkViewer.Show() + # Change color scale, has to be set to bytes=True. + vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True) + vtkViewer.Show() + # Set limits of values to view + vtkViewer.limits = [5.0,100.0] + vtkViewer.Show() \ No newline at end of file From c9195f0651962fb2c332c55f14b38696c7862c46 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 12:51:16 -0800 Subject: [PATCH 08/16] Made a notebook that works with updated version of vtkView --- notebooks/3DRenderingWithvtkTools.ipynb | 115 ------------------------ 1 file changed, 115 deletions(-) delete mode 100644 notebooks/3DRenderingWithvtkTools.ipynb diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb deleted file mode 100644 index 8d2c833d..00000000 --- a/notebooks/3DRenderingWithvtkTools.ipynb +++ /dev/null @@ -1,115 +0,0 @@ -{ - "metadata": { -<<<<<<< HEAD - "name": "3D rendering with vtkTools" -======= - "name": "" ->>>>>>> develop - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ - { - "cells": [ - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "import numpy as np, vtk\nimport SimPEG as simpeg", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 1 -======= - "input": [ - "import numpy as np, vtk\n", - "import SimPEG as simpeg" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 -======= - "input": [ - "#Make a mesh and model\n", - "x0 = np.zeros(3)\n", - "h1 = np.ones(20)*50\n", - "h2 = np.ones(10)*100\n", - "h3 = np.ones(5)*200\n", - "\n", - "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 6 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n# Make the vtk viewer object.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 -======= - "input": [ - "# Make a models that correspond to the cells, faces and edges.\n", - "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n", - "# Make the vtk viewer object.\n", - "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) \n", - " " - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 7 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 -======= - "input": [ - "# Show the image \n", - "vtkViewer.Show()\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": "*" ->>>>>>> develop - } - ], - "metadata": {} - } - ] -} From e0e64a6f9e720c0d5caf39384f7204c850eedbf1 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 12:52:55 -0800 Subject: [PATCH 09/16] Commit the new notebook for vtkView --- notebooks/VisualizeWithvtkView-updated..ipynb | 142 ++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 notebooks/VisualizeWithvtkView-updated..ipynb diff --git a/notebooks/VisualizeWithvtkView-updated..ipynb b/notebooks/VisualizeWithvtkView-updated..ipynb new file mode 100644 index 00000000..6dabae72 --- /dev/null +++ b/notebooks/VisualizeWithvtkView-updated..ipynb @@ -0,0 +1,142 @@ +{ + "metadata": { + "name": "VisualizeWithvtkView-updated." + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import SimPEG as simpeg, matplotlib as mpl" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The history saving thread hit an unexpected error (OperationalError('disk I/O error',)).History will not be written to the database.\n", + "Warning: mumps solver not available." + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Simple notebook of how to use vtkView to visualize SimPEG models. It will pop-up external vtk windows." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make a mesh and model\n", + "x0 = np.zeros(3)\n", + "h1 = np.ones(60)*50\n", + "h2 = np.ones(60)*100\n", + "h3 = np.ones(50)*200\n", + "\n", + "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "\n", + "# Make a models that correspond to the cells, faces and edges.\n", + "t = np.ones(mesh.nC)\n", + "t[10000:50000] = 100\n", + "t[100000:120000] = 100\n", + "t[100000:120000] = 50\n", + "# Make models called 'Test' for all with a range. \n", + "models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t},'F':{'Test':np.arange(0,mesh.nF)},'E':{'Test':np.arange(0,mesh.nE)}}\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make the vtk viewer object.\n", + "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)\n", + "# Set the .viewprop for which model to view\n", + "vtkViewer.viewprop = {'F':'Test'}\n", + "# Show the image\n", + "vtkViewer.Show()\n", + "\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Set subset of the mesh to view (remove padding)\n", + "vtkViewer.extent = [4,14,0,7,0,3]\n", + "vtkViewer.Show()\n", + "\n", + "# Change viewing property \n", + "vtkViewer.viewprop = {'C':'Model'}\n", + "# Set the color range\n", + "# Reset extent. Error check will reset the limits correctly.\n", + "vtkViewer.extent = [-1,1000,-1,1000,-1,1000]\n", + "# Set the range\n", + "vtkViewer.range = [0.,100.]\n", + "# Show\n", + "vtkViewer.Show()\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:116: UserWarning: Lower bounds smaller then 0\n", + " warnings.warn('Lower bounds smaller then 0')\n", + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:128: UserWarning: Upper bounds greater then number of cells\n", + " warnings.warn('Upper bounds greater then number of cells')\n", + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:137: UserWarning: Changed given extent from [-1, 1000, -1, 1000, -1, 1000] to [0, 59, 0, 59, 0, 49]\n", + " warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))\n" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Change color scale, has to be set to bytes=True.\n", + "vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)\n", + "vtkViewer.Show()\n", + "# Set limits of values to view \n", + "vtkViewer.limits = [5.0,100.0]\n", + "vtkViewer.Show()" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file From 60c1fc7c426ab128f1b2e3dc9e6f09d03c9c1440 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Wed, 27 Nov 2013 13:04:48 -0800 Subject: [PATCH 10/16] New notebook for vtkView example --- notebooks/VisualizeWithvtkView-updated.ipynb | 142 +++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 notebooks/VisualizeWithvtkView-updated.ipynb diff --git a/notebooks/VisualizeWithvtkView-updated.ipynb b/notebooks/VisualizeWithvtkView-updated.ipynb new file mode 100644 index 00000000..f4e8a3a7 --- /dev/null +++ b/notebooks/VisualizeWithvtkView-updated.ipynb @@ -0,0 +1,142 @@ +{ + "metadata": { + "name": "VisualizeWithvtkView-updated" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import SimPEG as simpeg, matplotlib as mpl" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The history saving thread hit an unexpected error (OperationalError('disk I/O error',)).History will not be written to the database.\n", + "Warning: mumps solver not available." + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Simple notebook of how to use vtkView to visualize SimPEG models. It will pop-up external vtk windows." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make a mesh and model\n", + "x0 = np.zeros(3)\n", + "h1 = np.ones(60)*50\n", + "h2 = np.ones(60)*100\n", + "h3 = np.ones(50)*200\n", + "\n", + "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "\n", + "# Make a models that correspond to the cells, faces and edges.\n", + "t = np.ones(mesh.nC)\n", + "t[10000:50000] = 100\n", + "t[100000:120000] = 100\n", + "t[100000:120000] = 50\n", + "# Make models called 'Test' for all with a range. \n", + "models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t},'F':{'Test':np.arange(0,mesh.nF)},'E':{'Test':np.arange(0,mesh.nE)}}\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make the vtk viewer object.\n", + "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)\n", + "# Set the .viewprop for which model to view\n", + "vtkViewer.viewprop = {'F':'Test'}\n", + "# Show the image\n", + "vtkViewer.Show()\n", + "\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Set subset of the mesh to view (remove padding)\n", + "vtkViewer.extent = [4,14,0,7,0,3]\n", + "vtkViewer.Show()\n", + "\n", + "# Change viewing property \n", + "vtkViewer.viewprop = {'C':'Model'}\n", + "# Set the color range\n", + "# Reset extent. Error check will reset the limits correctly.\n", + "vtkViewer.extent = [-1,1000,-1,1000,-1,1000]\n", + "# Set the range\n", + "vtkViewer.range = [0.,100.]\n", + "# Show\n", + "vtkViewer.Show()\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:116: UserWarning: Lower bounds smaller then 0\n", + " warnings.warn('Lower bounds smaller then 0')\n", + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:128: UserWarning: Upper bounds greater then number of cells\n", + " warnings.warn('Upper bounds greater then number of cells')\n", + "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:137: UserWarning: Changed given extent from [-1, 1000, -1, 1000, -1, 1000] to [0, 59, 0, 59, 0, 49]\n", + " warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))\n" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Change color scale, has to be set to bytes=True.\n", + "vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)\n", + "vtkViewer.Show()\n", + "# Set limits of values to view \n", + "vtkViewer.limits = [5.0,100.0]\n", + "vtkViewer.Show()" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file From 1b58820b5e040071c93f6c5319de83b2779007e6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 27 Nov 2013 14:30:27 -0700 Subject: [PATCH 11/16] Delete 3DRenderingWithvtkTools.ipynb.orig --- notebooks/3DRenderingWithvtkTools.ipynb.orig | 115 ------------------- 1 file changed, 115 deletions(-) delete mode 100644 notebooks/3DRenderingWithvtkTools.ipynb.orig diff --git a/notebooks/3DRenderingWithvtkTools.ipynb.orig b/notebooks/3DRenderingWithvtkTools.ipynb.orig deleted file mode 100644 index 8d2c833d..00000000 --- a/notebooks/3DRenderingWithvtkTools.ipynb.orig +++ /dev/null @@ -1,115 +0,0 @@ -{ - "metadata": { -<<<<<<< HEAD - "name": "3D rendering with vtkTools" -======= - "name": "" ->>>>>>> develop - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ - { - "cells": [ - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "import numpy as np, vtk\nimport SimPEG as simpeg", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 1 -======= - "input": [ - "import numpy as np, vtk\n", - "import SimPEG as simpeg" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 -======= - "input": [ - "#Make a mesh and model\n", - "x0 = np.zeros(3)\n", - "h1 = np.ones(20)*50\n", - "h2 = np.ones(10)*100\n", - "h3 = np.ones(5)*200\n", - "\n", - "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 6 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n# Make the vtk viewer object.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 -======= - "input": [ - "# Make a models that correspond to the cells, faces and edges.\n", - "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n", - "# Make the vtk viewer object.\n", - "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) \n", - " " - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 7 ->>>>>>> develop - }, - { - "cell_type": "code", - "collapsed": false, -<<<<<<< HEAD - "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 -======= - "input": [ - "# Show the image \n", - "vtkViewer.Show()\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": "*" ->>>>>>> develop - } - ], - "metadata": {} - } - ] -} From f9dd00ca8851a4c55472bb64f72e627198891dac Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 27 Nov 2013 15:32:50 -0700 Subject: [PATCH 12/16] Delete InteractiveView.py --- SimPEG/InteractiveView.py | 335 -------------------------------------- 1 file changed, 335 deletions(-) delete mode 100644 SimPEG/InteractiveView.py diff --git a/SimPEG/InteractiveView.py b/SimPEG/InteractiveView.py deleted file mode 100644 index ba4f1b63..00000000 --- a/SimPEG/InteractiveView.py +++ /dev/null @@ -1,335 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -from utils import mkvc - - -class TensorView(object): - """ - Provides viewing functions for TensorMesh - - This class is inherited by TensorMesh - """ - def __init__(self): - pass - - def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): - """ - Mesh.plotImage(I) - - Plots scalar fields on the given mesh. - - Input: - - :param numpy.array I: scalar field - - Optional Input: - - :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz - :param int figNum: number of figure to plot to - :param matplotlib.axes.Axes ax: axis to plot to - :param str direction: slice dimensions, 3D only ('x', 'y', 'z') - :param bool numbering: show numbering of slices, 3D only - :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_image_2D.py - :include-source: - - .. plot:: examples/mesh/plot_image_3D.py - :include-source: - """ - assert type(I) == np.ndarray, "I must be a numpy array" - assert type(numbering) == bool, "numbering must be a bool" - assert direction in ["x", "y","z"], "direction must be either x,y, or z" - - - if imageType == 'CC': - assert I.size == self.nC, "Incorrect dimensions for CC." - elif imageType == 'N': - assert I.size == self.nN, "Incorrect dimensions for N." - elif imageType == 'Fx': - if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') - elif imageType == 'Fy': - if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') - elif imageType == 'Fz': - if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') - elif imageType == 'Ex': - if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') - elif imageType == 'Ey': - if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') - elif imageType == 'Ez': - if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') - elif imageType[0] == 'E': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - ex, ey, ez = self.r(I,'E','E','M') - if plotAll or 'Ex' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(ex, imageType='Ex', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Ey' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(ey, imageType='Ey', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Ez' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(ez, imageType='Ez', ax=ax_z, **options) - pltNum +=1 - return - elif imageType[0] == 'F': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - fx, fy, fz = self.r(I,'F','F','M') - if plotAll or 'Fx' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(fx, imageType='Fx', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Fy' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(fy, imageType='Fy', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Fz' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(fz, imageType='Fz', ax=ax_z, **options) - pltNum +=1 - return - else: - raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") - - - if ax is None: - fig = plt.figure(figNum) - fig.clf() - ax = plt.subplot(111) - else: - assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" - fig = ax.figure - - if self.dim == 1: - if imageType == 'CC': - ph = ax.plot(self.vectorCCx, I, '-ro') - elif imageType == 'N': - ph = ax.plot(self.vectorNx, I, '-bs') - ax.set_xlabel("x") - ax.axis('tight') - elif self.dim == 2: - if imageType == 'CC': - C = I[:].reshape(self.n, order='F') - elif imageType == 'N': - C = I[:].reshape(self.n+1, order='F') - C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) - elif imageType == 'Fx': - C = I[:].reshape(self.nFx, order='F') - C = 0.5*(C[:-1, :] + C[1:, :] ) - elif imageType == 'Fy': - C = I[:].reshape(self.nFy, order='F') - C = 0.5*(C[:, :-1] + C[:, 1:] ) - elif imageType == 'Ex': - C = I[:].reshape(self.nEx, order='F') - C = 0.5*(C[:,:-1] + C[:,1:] ) - elif imageType == 'Ey': - C = I[:].reshape(self.nEy, order='F') - C = 0.5*(C[:-1,:] + C[1:,:] ) - - ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) - ax.axis('tight') - ax.set_xlabel("x") - ax.set_ylabel("y") - - elif self.dim == 3: - if direction == 'z': - - # get copy of image and average to cell-centres is necessary - if imageType == 'CC': - Ic = I[:].reshape(self.n, order='F') - elif imageType == 'N': - Ic = I[:].reshape(self.n+1, order='F') - Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) - elif imageType == 'Fx': - Ic = I[:].reshape(self.nFx, order='F') - Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) - elif imageType == 'Fy': - Ic = I[:].reshape(self.nFy, order='F') - Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) - elif imageType == 'Fz': - Ic = I[:].reshape(self.nFz, order='F') - Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) - elif imageType == 'Ex': - Ic = I[:].reshape(self.nEx, order='F') - Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) - elif imageType == 'Ey': - Ic = I[:].reshape(self.nEy, order='F') - Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) - elif imageType == 'Ez': - Ic = I[:].reshape(self.nEz, order='F') - Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) - - # determine number oE slices in x and y dimension - nX = np.ceil(np.sqrt(self.nCz)) - nY = np.ceil(self.nCz/nX) - - # allocate space for montage - nCx = self.nCx - nCy = self.nCy - - C = np.zeros((nX*nCx,nY*nCy)) - - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] - else: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan - - C = np.ma.masked_where(np.isnan(C), C) - xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] - yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] - # Plot the mesh - ph = ax.pcolormesh(xx, yy, C.T) - # Plot the lines - gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) - gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) - # Repeat and seperate with NaN - gxX = np.c_[gx, gx, gx+np.nan].ravel() - gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() - gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() - gyY = np.c_[gy, gy, gy+np.nan].ravel() - ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) - ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) - ax.axis('tight') - - if numbering: - pad = np.sum(self.hx)*0.04 - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, - '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') - - ax.set_title(imageType) - if showIt: plt.show() - return ph - - def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. - - :param bool nodes: plot nodes - :param bool faces: plot faces - :param bool centers: plot centers - :param bool edges: plot edges - :param bool lines: plot lines connecting nodes - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_grid_2D.py - :include-source: - - .. plot:: examples/mesh/plot_grid_3D.py - :include-source: - """ - if self.dim == 1: - fig = plt.figure(1) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - ax.hold(True) - ax.plot(xn, np.ones(np.shape(xn)), 'bs') - ax.plot(xc, np.ones(np.shape(xc)), 'ro') - ax.plot(xn, np.ones(np.shape(xn)), 'k--') - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - if showIt: plt.show() - elif self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - xs1 = self.gridFx - xs2 = self.gridFy - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') - if faces: - ax.plot(xs1[:, 0], xs1[:, 1], 'g>') - ax.plot(xs2[:, 0], xs2[:, 1], 'g^') - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] - plt.plot(X, Y) - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - if showIt: plt.show() - elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') - xn = self.gridN - xc = self.gridCC - xfs1 = self.gridFx - xfs2 = self.gridFy - xfs3 = self.gridFz - - xes1 = self.gridEx - xes2 = self.gridEy - xes3 = self.gridEz - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) - if faces: - ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) - ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) - ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) - if edges: - ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) - ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) - ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() - Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() - Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() - X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() - X = np.r_[X1, X2, X3] - Y = np.r_[Y1, Y2, Y3] - Z = np.r_[Z1, Z2, Z3] - plt.plot(X, Y, 'b-', zs=Z) - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - ax.set_zlabel('x3') - if showIt: plt.show() From 2a79e31d1f011707e6867d4024a51df98d892753 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 27 Nov 2013 15:33:18 -0700 Subject: [PATCH 13/16] Delete VisualizeWithvtkView-updated..ipynb --- notebooks/VisualizeWithvtkView-updated..ipynb | 142 ------------------ 1 file changed, 142 deletions(-) delete mode 100644 notebooks/VisualizeWithvtkView-updated..ipynb diff --git a/notebooks/VisualizeWithvtkView-updated..ipynb b/notebooks/VisualizeWithvtkView-updated..ipynb deleted file mode 100644 index 6dabae72..00000000 --- a/notebooks/VisualizeWithvtkView-updated..ipynb +++ /dev/null @@ -1,142 +0,0 @@ -{ - "metadata": { - "name": "VisualizeWithvtkView-updated." - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ - { - "cells": [ - { - "cell_type": "code", - "collapsed": false, - "input": [ - "import SimPEG as simpeg, matplotlib as mpl" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "The history saving thread hit an unexpected error (OperationalError('disk I/O error',)).History will not be written to the database.\n", - "Warning: mumps solver not available." - ] - }, - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "\n" - ] - } - ], - "prompt_number": 1 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Simple notebook of how to use vtkView to visualize SimPEG models. It will pop-up external vtk windows." - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "# Make a mesh and model\n", - "x0 = np.zeros(3)\n", - "h1 = np.ones(60)*50\n", - "h2 = np.ones(60)*100\n", - "h3 = np.ones(50)*200\n", - "\n", - "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", - "\n", - "# Make a models that correspond to the cells, faces and edges.\n", - "t = np.ones(mesh.nC)\n", - "t[10000:50000] = 100\n", - "t[100000:120000] = 100\n", - "t[100000:120000] = 50\n", - "# Make models called 'Test' for all with a range. \n", - "models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t},'F':{'Test':np.arange(0,mesh.nF)},'E':{'Test':np.arange(0,mesh.nE)}}\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "# Make the vtk viewer object.\n", - "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)\n", - "# Set the .viewprop for which model to view\n", - "vtkViewer.viewprop = {'F':'Test'}\n", - "# Show the image\n", - "vtkViewer.Show()\n", - "\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "# Set subset of the mesh to view (remove padding)\n", - "vtkViewer.extent = [4,14,0,7,0,3]\n", - "vtkViewer.Show()\n", - "\n", - "# Change viewing property \n", - "vtkViewer.viewprop = {'C':'Model'}\n", - "# Set the color range\n", - "# Reset extent. Error check will reset the limits correctly.\n", - "vtkViewer.extent = [-1,1000,-1,1000,-1,1000]\n", - "# Set the range\n", - "vtkViewer.range = [0.,100.]\n", - "# Show\n", - "vtkViewer.Show()\n" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stderr", - "text": [ - "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:116: UserWarning: Lower bounds smaller then 0\n", - " warnings.warn('Lower bounds smaller then 0')\n", - "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:128: UserWarning: Upper bounds greater then number of cells\n", - " warnings.warn('Upper bounds greater then number of cells')\n", - "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:137: UserWarning: Changed given extent from [-1, 1000, -1, 1000, -1, 1000] to [0, 59, 0, 59, 0, 49]\n", - " warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))\n" - ] - } - ], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "# Change color scale, has to be set to bytes=True.\n", - "vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)\n", - "vtkViewer.Show()\n", - "# Set limits of values to view \n", - "vtkViewer.limits = [5.0,100.0]\n", - "vtkViewer.Show()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 - } - ], - "metadata": {} - } - ] -} \ No newline at end of file From 51c00e79bfe61ef4e3738286d7b4bde9e7176d28 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 2 Dec 2013 16:46:15 -0800 Subject: [PATCH 14/16] Update README.md --- README.md | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 172cf62b..e107fea3 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,13 @@ -simpeg -====== +![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png) Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. +The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features: + +* modular with respect to the spacial discretization, optimization routine, and geophysical problem +* built with the inverse problem in mind +* provides a framework for geophysical and hydrogeologic problems +* supports 1D, 2D and 3D problems +* designed for large-scale inversions [![Build Status](https://travis-ci.org/simpeg/simpeg.png)](https://travis-ci.org/simpeg/simpeg) From a2259b78112c72c330cef20350aaebb6b0407284 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 11:48:17 -0800 Subject: [PATCH 15/16] updates to Newton Root Finding. --- SimPEG/inverse/Optimize.py | 27 ++++++++++++++++++--------- SimPEG/tests/test_optimizers.py | 10 ++++++++++ 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 7e934e13..6edacabd 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -795,12 +795,23 @@ class NewtonRoot(object): setKwargs(self, **kwargs) def root(self, fun, x): + """root(fun, x) + + Function Should have the form:: + + def evalFunction(x, return_g=False): + out = (f,) + if return_g: + out += (g,) + return out if len(out) > 1 else out[0] + + """ if self.comments: print 'Newton Method:\n' self._iter = 0 while True: - [r,J] = fun(x); + r, J = fun(x, return_g=True) if self.solveTol == 0: Jinv = Solver(J) dh = - Jinv.solve(r) @@ -812,13 +823,12 @@ class NewtonRoot(object): muLS = 1. LScnt = 1 xt = x + dh - rt, Jt = fun(xt) # TODO: get rid of Jt + rt = fun(xt, return_g=False) if self.comments: print '\tLinesearch:\n' # Enter Linesearch while True and self.doLS: - if self.comments: - print '\t\tResid: %e\n'%norm(rt) + if self.comments: print '\t\tResid: %e\n'%norm(rt) if norm(rt) <= norm(r) or norm(rt) < self.tol: break @@ -827,10 +837,9 @@ class NewtonRoot(object): print '.' if LScnt > self.maxLS: print 'Newton Method: Line search break.' - root = NaN - return + return None xt = x + muLS*dh - rt, Jt = fun(xt) # TODO: get rid of Jt + rt = fun(xt, return_g=False) x = xt self._iter += 1 @@ -854,7 +863,7 @@ if __name__ == '__main__': print 'test the newtonRoot finding.' - fun = lambda x: (np.sin(x), sdiag(np.cos(x))) + fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), sdiag( np.cos(x) ) ) x = np.array([np.pi-0.3, np.pi+0.1, 0]) - pnt = NewtonRoot(comments=False).root(fun,x) + pnt = NewtonRoot(comments=True).root(fun,x) print pnt diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/tests/test_optimizers.py index 6fff75d3..654d1804 100644 --- a/SimPEG/tests/test_optimizers.py +++ b/SimPEG/tests/test_optimizers.py @@ -50,5 +50,15 @@ class TestOptimizers(unittest.TestCase): print 'x_true: ', x_true self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) + def test_NewtonRoot(self): + fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), sdiag( np.cos(x) ) ) + x = np.array([np.pi-0.3, np.pi+0.1, 0]) + xopt = inverse.NewtonRoot(comments=False).root(fun,x) + x_true = np.array([np.pi,np.pi,0]) + print 'Newton Root Finding' + print 'xopt: ', xopt + print 'x_true: ', x_true + self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) + if __name__ == '__main__': unittest.main() From ccb4d4052de7b59ace47eeb8a10af81454a4a61a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 4 Dec 2013 15:31:11 -0800 Subject: [PATCH 16/16] Change callHooks structure. --- SimPEG/inverse/Inversion.py | 42 ++++++++++----------------------- SimPEG/inverse/Optimize.py | 40 ++++++++++++------------------- SimPEG/tests/test_optimizers.py | 2 +- SimPEG/utils/__init__.py | 38 +++++++++++++++++++++++++---- 4 files changed, 62 insertions(+), 60 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 160678aa..a2ef2513 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -85,25 +85,19 @@ class BaseInversion(object): if self.stoppingCriteria(): break self.printDone() + self.finish() + return self.m + @callHooks('startup') def startup(self, m0): """ **startup** is called at the start of any new run call. - If you have things that also need to run on startup, you can create a method:: - - def _startup*(self, x0): - pass - - Where the * can be any string. If present, _startup* will be called at the start of the default startup call. - You may also completely overwrite this function. - :param numpy.ndarray x0: initial x :rtype: None :return: None """ - callHooks(self,'startup',m0) if not hasattr(self.reg, '_mref'): print 'Regularization has not set mref. SimPEG will set it to m0.' @@ -115,43 +109,25 @@ class BaseInversion(object): self.phi_d_last = np.nan self.phi_m_last = np.nan + @callHooks('doStartIteration') def doStartIteration(self): """ **doStartIteration** is called at the end of each run iteration. - If you have things that also need to run at the end of every iteration, you can create a method:: - - def _doStartIteration*(self): - pass - - Where the * can be any string. If present, _doStartIteration* will be called at the start of the default doStartIteration call. - You may also completely overwrite this function. - :rtype: None :return: None """ - callHooks(self,'doStartIteration') - self._beta = self.getBeta() + @callHooks('doEndIteration') def doEndIteration(self): """ **doEndIteration** is called at the end of each run iteration. - If you have things that also need to run at the end of every iteration, you can create a method:: - - def _doEndIteration*(self): - pass - - Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call. - You may also completely overwrite this function. - :rtype: None :return: None """ - callHooks(self,'doEndIteration') - # store old values self.phi_d_last = self.phi_d self.phi_m_last = self.phi_m @@ -213,6 +189,14 @@ class BaseInversion(object): """ printStoppers(self, self.stoppers) + @callHooks('finish') + def finish(self): + """finish() + + **finish** is called at the end of the optimization. + """ + pass + @timeIt def evalFunction(self, m, return_g=True, return_H=True): """evalFunction(m, return_g=True, return_H=True) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 6edacabd..80a03632 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -155,6 +155,7 @@ class Minimize(object): doEndIteration(xt) printDone() + finish() return xc """ self.evalFunction = evalFunction @@ -175,6 +176,7 @@ class Minimize(object): self.doEndIteration(xt) self.printDone() + self.finish() return self.xc @@ -188,6 +190,7 @@ class Minimize(object): def parent(self, value): self._parent = value + @callHooks('startup') def startup(self, x0): """ **startup** is called at the start of any new minimize call. @@ -198,19 +201,10 @@ class Minimize(object): xc = x0 _iter = _iterLS = 0 - If you have things that also need to run on startup, you can create a method:: - - def _startup*(self, x0): - pass - - Where the * can be any string. If present, _startup* will be called at the start of the default startup call. - You may also completely overwrite this function. - :param numpy.ndarray x0: initial x :rtype: None :return: None """ - callHooks(self,'startup',x0) self._iter = 0 self._iterLS = 0 @@ -222,6 +216,7 @@ class Minimize(object): self.x_last = x0 @count + @callHooks('doStartIteration') def doStartIteration(self): """doStartIteration() @@ -230,7 +225,8 @@ class Minimize(object): :rtype: None :return: None """ - callHooks(self,'doStartIteration') + pass + def printInit(self, inLS=False): """ @@ -244,6 +240,7 @@ class Minimize(object): name = self.name if not inLS else self.nameLS printTitles(self, self.printers if not inLS else self.printersLS, name, pad) + @callHooks('printIter') def printIter(self, inLS=False): """ **printIter** is called directly after function evaluations. @@ -252,8 +249,6 @@ class Minimize(object): parent.printIter function and call that. """ - callHooks(self,'printIter',inLS) - pad = ' '*10 if inLS else '' printLine(self, self.printers if not inLS else self.printersLS, pad=pad) @@ -270,6 +265,11 @@ class Minimize(object): stoppers = self.stoppers if not inLS else self.stoppersLS printStoppers(self, stoppers, pad='', stop=stop, done=done) + + def finish(self): + pass + + def stoppingCriteria(self, inLS=False): if self._iter == 0: self.f0 = self.f @@ -277,6 +277,7 @@ class Minimize(object): return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) @timeIt + @callHooks('projection') def projection(self, p): """projection(p) @@ -288,7 +289,6 @@ class Minimize(object): :rtype: numpy.ndarray :return: p, projected search direction """ - callHooks(self,'projection',p) return p @timeIt @@ -402,6 +402,7 @@ class Minimize(object): return p, False @count + @callHooks('doEndIteration') def doEndIteration(self, xt): """doEndIteration(xt) @@ -411,21 +412,10 @@ class Minimize(object): self.xc must be updated in this code. - - If you have things that also need to run at the end of every iteration, you can create a method:: - - def _doEndIteration*(self, xt): - pass - - Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call. - You may also completely overwrite this function. - :param numpy.ndarray xt: tested new iterate that ensures a descent direction. :rtype: None :return: None """ - callHooks(self,'doEndIteration',xt) - # store old values self.f_last = self.f self.x_last, self.xc = self.xc, xt @@ -630,7 +620,7 @@ class ProjectedGradient(Minimize, Remember): if self.debug: print 'doEndIteration.ProjGrad, f_current_decrease: ', f_current_decrease if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max - if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD + if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingPG class BFGS(Minimize, Remember): diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/tests/test_optimizers.py index 654d1804..0253852c 100644 --- a/SimPEG/tests/test_optimizers.py +++ b/SimPEG/tests/test_optimizers.py @@ -32,7 +32,7 @@ class TestOptimizers(unittest.TestCase): self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) def test_ProjGradient_quadraticBounded(self): - PG = inverse.ProjectedGradient() + PG = inverse.ProjectedGradient(debug=True) PG.lower, PG.upper = -2, 2 xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0])) x_true = np.array([2.,2.]) diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 622557d6..19aecec6 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -23,7 +23,7 @@ def hook(obj, method, name=None, overwrite=False, silent=False): If name is None, the name of the method is used. """ - if name is None: + if name is None: name = method.__name__ if name == '': raise Exception('Must provide name to hook lambda functions.') @@ -87,11 +87,39 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'): print pad + stopper['str'] % (l<=r,l,r) print pad + "%s%s%s" % ('-'*25,done,'-'*25) -def callHooks(obj, match, *args, **kwargs): - for method in [posible for posible in dir(obj) if ('_'+match) in posible]: - if getattr(obj,'debug',False): print (match+' is calling self.'+method) - getattr(obj,method)(*args, **kwargs) +def callHooks(match): + """ + Use this to wrap a funciton:: + @callHooks('doEndIteration') + def doEndIteration(self): + pass + + This will call everything named _doEndIteration* at the beginning of the function call. + """ + def callHooksWrap(f): + @wraps(f) + def wrapper(self,*args,**kwargs): + + for method in [posible for posible in dir(self) if ('_'+match) in posible]: + if getattr(self,'debug',False): print (match+' is calling self.'+method) + getattr(self,method)(*args, **kwargs) + + return f(self,*args,**kwargs) + + extra = """ + If you have things that also need to run in the method %s, you can create a method:: + + def _%s*(self, ... ): + pass + + Where the * can be any string. If present, _%s* will be called at the start of the default %s call. + You may also completely overwrite this function. + """ % (match, match, match, match) + + wrapper.__doc__ += extra + return wrapper + return callHooksWrap class Counter(object):