title: Converting VTK UnstructuredGrid to StructuredGrid
body: <p>The background to my problem is that I have a 3D structure saved in a .vtk file that I need to manipulate (dilate, erode, etc.). The following code snippets are designed to be run sequentially, i.e. if you run them one after the other, there should be no problems (apart from those I mention!).</p>
  
  <p>I'm very new to VTK, so apologies for any very basic mistakes!</p>
  
  <h1>Problem</h1>
  
  <p>My problem stems from a problem with SimpleITK, wherein it is unable to read UnstructuredGrid or PolyData:</p>
  
  <pre><code>In [1]: import SimpleITK as sitk
  In [2]: img_vtk = sitk.ReadImage(file_vtk)
  Traceback (most recent call last):
  
    File "&lt;ipython-input-52-435ce999db50&gt;", line 1, in &lt;module&gt;
      img_vtk = sitk.ReadImage(file_vtk)
  
    File "/usr/local/lib/python3.5/dist-packages/SimpleITK/SimpleITK.py", line 8614, in ReadImage
      return _SimpleITK.ReadImage(*args)
  
  RuntimeError: Exception thrown in SimpleITK ReadImage: /tmp/SimpleITK/Code/IO/src/sitkImageReaderBase.cxx:97:
  sitk::ERROR: Unable to determine ImageIO reader for "/data/ROMPA_MRIandSeg/09S/Analysis/1_model/clip_dilate.vtk"
  </code></pre>
  
  <p>SimpleITK can, however, read StructuredGrid, so I tried to solve this by reading using VTK and converting.</p>
  
  <pre><code>import vtk
  reader = vtk.vtkGenericDataObjectReader() # Using generic to allow it to match either Unstructured or PolyData
  reader.SetFileName(file_vtk)
  reader.Update()
  output = reader.GetOutput()
  </code></pre>
  
  <p>However, from that point on, every method I've tried seems to have failed.</p>
  
  <h1>Proposed Solutions</h1>
  
  <h2>Conversion to numpy, then conversion to sitk image</h2>
  
  <p>I attempted to convert it to a numpy array (), then interpolate a regular grid, with a dummy variable of 1 to specify the values on the structure.</p>
  
  <pre><code>from vtk.utils import numpy_support
  import scipy.interpolate
  import numpy as np
  
  nparray = numpy_support.vtk_to_numpy(output.GetPointData().GetArray(0))
  
  output_bounds = output.GetBounds()
  x_grid = range(math.floor(output_bounds[0]),math.ceil(output_bounds[1]),1)
  y_grid = range(math.floor(output_bounds[2]),math.ceil(output_bounds[3]),1)
  z_grid = range(math.floor(output_bounds[4]),math.ceil(output_bounds[5]),1)
  grid = list()
  for x in x_grid:
     for y in y_grid:
        for z in z_grid:
           grid.append((x,y,z))
  dummy = np.array([1 for i in range(nparray.shape[0])])
  npgrid = scipy.interpolate.griddata(nparray,dummy,grid,fill_value=0)
  
  npgrid.reshape(len(x_grid),len(y_grid),len(z_grid))
  img = sitk.GetImageFromArray(npgrid)
  sitk.WriteImage(img,file_out)
  </code></pre>
  
  <p>However, when I load this in ParaView, a bounding box is displayed for the output, but a contour of the output is empty.</p>
  
  <h2>Using <code>ShepardMethod</code></h2>
  
  <p>I attempted to interpolate using the built-in <code>ShepardMethod</code>, after converting the UnstructuredGrid to PolyData (as I'd mostly seen ShepardMethod being applied to PolyData):</p>
  
  <pre><code>bounds = output.GetBounds()
  spacings = [1.0,1.0,1.0] # arbitrary spacing
  dimensions = [0,0,0]
  for i,spacing in enumerate(spacings):
      dimensions[i] = int(math.ceil((bounds[i*2 + 1]-bounds[i*2])/spacing))
  
  vtkPoints = vtk.vtkPoints()
  for i in range(0,nparray.shape[0]):
     x=nparray[i,0]
     y=nparray[i,1]
     z=nparray[i,2]
     p=[x,y,z]
     vtkPoints.InsertNextPoint(p)
  poly = vtk.vtkPolyData()
  poly.SetPoints(vtkPoints)    
  
  shepard = vtk.vtkShepardMethod()
  shepard.SetInputData(poly)
  shepard.SetSampleDimensions(dimensions)
  shepard.SetModelBounds(output.GetBounds())
  shepard.Update()
  shepard_data = shepard.GetOutput().GetPointData().GetArray(0)
  
  shepard_numpy = numpy_support.vtk_to_numpy(shepard_data)
  shepard_numpy = shepard_numpy.reshape(dimensions[0],dimensions[1],dimensions[2])
  shepard_img = sitk.GetImageFromArray(shepard_numpy)
  sitk.WriteImage(shepard_img,file_out)
  </code></pre>
  
  <p>As with the numpy effort above, this provided a bounding box in ParaView. Applying a contour provided a structure of two triangles, i.e. next to nothing seems to have been successfully written. Alternatively, I attempted to write the output directly using VTK.</p>
  
  <pre><code>shepard_data = shepard.GetOutput()
  shepard_grid = vtk.vtkImageToStructuredGrid()
  shepard_grid.SetInputData(shepard_data)
  shepard_grid.Update()
  
  writer = vtk.vtkStructuredGridWriter()
  writer.SetFileName(file_out)
  writer.SetInputData(shepard_grid.GetOutput())
  writer.Write()
  </code></pre>
  
  <p>This produced the same output as before.</p>
  
  <h2>Using <code>ProbeFilter</code></h2>
  
  <p>I tried the above using <code>ProbeFilter</code> instead (with both conversion to <code>numpy</code> and writing directly). Unfortunately, the output was the same as above.</p>
  
  <pre><code>mesh = vtk.vtkStructuredGrid()
  mesh.SetDimensions(dimensions)
  
  probe = vtk.vtkProbeFilter()
  probe.SetInputData(mesh)
  probe.SetSourceData(output)
  probe.Update()
  probe_out = probe.GetOutput()
  
  writer = vtk.vtkStructuredGridWriter()
  writer.SetFileName(file_out)
  writer.SetInputData(probe.GetOutput())
  writer.Write()
  
  probe_data = probe.GetOutput().GetPointData().GetArray(0)
  probe_numpy = numpy_support.vtk_to_numpy(probe_data)
  probe_numpy = probe_numpy.reshape(dimensions[0],dimensions[1],dimensions[2])
  probe_img = sitk.GetImageFromArray(probe_numpy)
  </code></pre>
  
  <p>sitk.WriteImage(probe_img,file_out)</p>
  
  <p>However, this seemed to produce no viable output (<code>vtkStructuredGridWriter</code> produced an empty file, and <code>probe_numpy</code> was empty).</p>
  
  <h2>Changing ParaView output</h2>
  
  <p>My original data comes from a structuredGrid .vtk file, that I open using ParaView, and then clip to remove structures that aren't required in the mesh. Saving the output saves an unstructuredGrid, and I have been unable to figure out whether I can change that, and avoid this mess in the first place!</p>