Hello, I am using VTK for building a function of cropping freehand. My code is below (source: qSlicerSegmentEditorScissorsEffect.cxx)
# Step 1: Drawing a 2D contour on the screen
# Step 2: Mapping display space points to world positions
# Step 3: Take 2d contour as polydata line, and extrude surfaces from the near clipping plane to the far clipping plane
closedSurfaceStrips = vtk.vtkCellArray()
# Create cells by specifying a count of total points to be inserted,
# and then adding points one at a time using method InsertCellPoint()
closedSurfaceStrips.InsertNextCell(numberOfPoints * 2 + 2)
for i in range(numberOfPoints * 2):
closedSurfaceStrips.InsertCellPoint(i)
closedSurfaceStrips.InsertCellPoint(0)
closedSurfaceStrips.InsertCellPoint(1)
# Front cap
closedSurfacePolys = vtk.vtkCellArray()
closedSurfacePolys.InsertNextCell(numberOfPoints)
for i in range(numberOfPoints):
closedSurfacePolys.InsertCellPoint(i * 2)
# Back cap
closedSurfacePolys.InsertNextCell(numberOfPoints)
for i in range(numberOfPoints):
closedSurfacePolys.InsertCellPoint(i * 2 + 1)
# Construct polydata
closedSurfacePolyData = vtk.vtkPolyData()
closedSurfacePolyData.SetPoints(closedSurfacePoints)
closedSurfacePolyData.SetStrips(closedSurfaceStrips)
closedSurfacePolyData.SetPolys(closedSurfacePolys)
# Step 4: Take the polydata from step 3, and do a 3D stencil filter
brushPolyDataNormals = vtk.vtkPolyDataNormals()
brushPolyDataNormals.AutoOrientNormalsOn()
worldToModifierLabelmapIjkTransformer = vtk.vtkTransformPolyDataFilter()
worldToModifierLabelmapIjkTransform = vtk.vtkTransform()
worldToModifierLabelmapIjkTransformer.SetTransform(worldToModifierLabelmapIjkTransform)
worldToModifierLabelmapIjkTransformer.SetInputConnection(brushPolyDataNormals.GetOutputPort())
brushPolyDataToStencil = vtk.vtkPolyDataToImageStencil()
brushPolyDataToStencil .SetOutputOrigin(0, 0, 0)
brushPolyDataToStencil.SetOutputSpacing(1, 1, 1)
brushPolyDataToStencil.SetInputConnection(worldToModifierLabelmapIjkTransformer.GetOutputPort())
brushPolyDataNormals.SetInputData(closedSurfacePolyData)
brushPolyDataNormals.Update()
brushModel_ModifierLabelmapIjk = worldToModifierLabelmapIjkTransformer.GetOutput() # vtkPolyData object
boundsIjk = brushModel_ModifierLabelmapIjk.GetBounds()
brushPolyDataToStencil.SetOutputWholeExtent(
math.floor(boundsIjk[0]) - 1, math.ceil(boundsIjk[1]) + 1,
math.floor(boundsIjk[2]) - 1, math.ceil(boundsIjk[3]) + 1,
math.floor(boundsIjk[4]) - 1, math.ceil(boundsIjk[5]) + 1,
)
# Step 5: Convert from Image Stencil Data to Image Data
stencilToImage = vtk.vtkImageStencilToImage()
stencilToImage.SetInputData(self.brushPolyDataToStencil.GetOutput())
stencilToImage.SetInsideValue(1)
stencilToImage.SetOutsideValue(0)
stencilToImage.SetOutputScalarType(vtk.VTK_SHORT) # vtk.VTK_SHORT: [-32768->32767]
stencilToImage.Update()
# Step 6: Create a Image Data extends origin image data some properties, all scalar values equal 0
baseImage = vtk.vtkImageData()
baseImage.SetExtent(self.imageData.GetExtent())
baseImage.SetOrigin(self.imageData.GetOrigin())
baseImage.SetSpacing(self.imageData.GetSpacing())
baseImage.SetDirectionMatrix(self.imageData.GetDirectionMatrix())
baseImage.AllocateScalars(self.imageData.GetScalarType(), 1)
# Step 7: Filter scalar values inside, all inside values equal 1, and outside equal 0
modifierImage = stencilToImage.GetOutput()
baseExtent = baseImage.GetExtent()
updateExtent = [v for v in baseExtent]
modifierExtent = modifierImage.GetExtent()
for idx in range(3):
if modifierExtent[idx * 2] > updateExtent[idx * 2]:
updateExtent[idx * 2] = modifierExtent[idx * 2]
if modifierExtent[idx * 2 + 1] < updateExtent[idx * 2 + 1]:
updateExtent[idx * 2 + 1] = modifierExtent[idx * 2 + 1]
baseImageModified = False
for z in range(updateExtent[4], updateExtent[5] + 1):
for y in range(updateExtent[2], updateExtent[3] + 1):
for x in range(updateExtent[0], updateExtent[1] + 1):
if modifierImage.GetScalarComponentAsFloat(x, y, z, 0) > baseImage.GetScalarComponentAsFloat(x, y, z, 0):
baseImage.SetScalarComponentFromFloat(x, y, z, 0, modifierImage.GetScalarComponentAsFloat(x, y, z, 0))
baseImageModified = True
if baseImageModified:
baseImage.Modified()
# Step 8: Apply mask volume
maskToStencil = vtk.vtkImageToImageStencil()
maskToStencil.ThresholdByLower(0)
maskToStencil.SetInputData(baseImage)
maskToStencil.Update()
stencil = vtk.vtkImageStencil()
stencil.SetInputData(self.imageData)
stencil.SetStencilConnection(maskToStencil.GetOutputPort())
stencil.Update()
# Step 9: Render the new volume
newVolume = stencil.GetOutput()
When my origin image data have another spacing, for example: (0.703125, 0.703125, 1.25) after applying a mask volume, the cropped segment will deviate.What is the effect of spacing when I crop?
