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.