# Need Help with Image Arithmetic (vtkImageMathematics, subtraction, division by a scalar)

**URL:** https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716
**Category:** Support
**Created:** [September 6, 2019, 10:22pm UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716 "2019-09-06T22:22:10Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![MichelleKline](https://discourse.vtk.org/user_avatar/discourse.vtk.org/michellekline/32/2155_2.png) [@MichelleKline](https://discourse.vtk.org/u/MichelleKline)
#### Post date: [September 6, 2019, 10:22pm UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716/1 "2019-09-06T22:22:10Z")

</div>

Hi!

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()

Thanks in advance for any feedback.

Michelle Kline

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [September 7, 2019, 1:13am UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716/2 "2019-09-07T01:13:53Z")

</div>

In Python, I would perform these operations using numpy. See description of the API here: [https://blog.kitware.com/improved-vtk-numpy-integration/](https://blog.kitware.com/improved-vtk-numpy-integration/)

---

<div class="post-metadata">

### Author: ![MichelleKline](https://discourse.vtk.org/user_avatar/discourse.vtk.org/michellekline/32/2155_2.png) [@MichelleKline](https://discourse.vtk.org/u/MichelleKline)
#### Post date: [September 8, 2019, 11:43pm UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716/3 "2019-09-08T23:43:37Z")

</div>

Thank you, Andras. I will try to do my image arithmetic using numpy instead.

- Michelle

---

<div class="post-metadata">

### Author: ![MichelleKline](https://discourse.vtk.org/user_avatar/discourse.vtk.org/michellekline/32/2155_2.png) [@MichelleKline](https://discourse.vtk.org/u/MichelleKline)
#### Post date: [September 10, 2019, 1:00pm UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716/4 "2019-09-10T13:00:41Z")

</div>

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?

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [September 10, 2019, 4:48pm UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716/5 "2019-09-10T16:48:40Z")

</div>

You can use this function to get read-write access to voxels of a vtkImageData object as a numpy array:

```python
def vtk_image_to_numpy_array(vimage):
    import vtk.util.numpy_support
    narray_shape = tuple(reversed(vimage.GetDimensions()))
    narray = vtk.util.numpy_support.vtk_to_numpy(vimage.GetPointData().GetScalars()).reshape(narray_shape)
    return narray

```

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.

---

<div class="post-metadata">

### Author: ![MichelleKline](https://discourse.vtk.org/user_avatar/discourse.vtk.org/michellekline/32/2155_2.png) [@MichelleKline](https://discourse.vtk.org/u/MichelleKline)
#### Post date: [September 10, 2019, 7:26pm UTC](https://discourse.vtk.org/t/need-help-with-image-arithmetic-vtkimagemathematics-subtraction-division-by-a-scalar/1716/6 "2019-09-10T19:26:04Z")

</div>

Thanks, Andras! I understand, and this is very helpful.

- Michelle
