I have a Python script which reads two DICOM images. I want to subtract one image from the other, and then divide the resulting difference by a scalar. In other words:
output_image = (input_image2 - input_image1) / some_float
I am trying to use vtkImageMathematics to accomplish this, but I am running in to problems.
First, some cells in input_image1 contain values that are larger than input_image2. When I use vtkImageMathematics.SetOperationToSubtract(), the values of these cells in the output_image appears to have overflowed the maximum allowable value the cell’s data type. I think the cell data type must be unsigned short, because, for example, 30 - 40 results in 65525. Is there a way to control the behavior of vtkImageMathematics.SetOperationToSubtract() such that negative values become 0?
Second, how to I apply another arithmetic operation to the output_image?
Here is the subtraction:
imageMath = vtk.vtkImageMathematics()
imageMath .SetOperationToSubtract()
imageMath .SetInput1Data(input_image2)
imageMath.SetInput2Data(input_image1.GetOutput())
imageMath.Update()
How do I apply division by a scalar? the vtkImageMathematics object doesn’t seem to operate on itself, it seems to require an input from some other object. Also, there is no “SetOperationToDivisionByK” or anything similar. Here is what I unsuccessfully tried:
imageMath2= vtk.vtkImageMathematics()
imageMath2.SetInput1Data(imageMath.GetOutput())
K = -0.89 # for example
imageMath2.SetConstantK(1/K)
imageMath2.SetOperationToMultiplyByK()
imageMath2.Update()
Hello again, Andras. I have tried using numpy to subtract one vtkImageData from another. But I am having trouble. In order to use numpy, I need to first convert the vtkImageData to a numpy array, then perform the arithmetic, then finally convert the results back to a vtkImageData.
The converting is not straightforward. Does it sound to you like I am trying to go about this correctly?
There is no need to “convert back” from numpy array to vtkImageData, just use this function to get read-write access to all the input and output vtkImageData.
Notes:
Voxel indices in the numpy array are reversed (k,j,i vs. i,j,k).
When you set values in the output image, you don’t want numpy to reallocate the array but instead set values in the existing array, so use the [:] operator. Example: narrayOutput[:] = narrayInput*2+15
If the output image is used in a pipeline: call vimageOutput.GetPointData().GetScalars().Modified() when all your modifications are completed.