# Ray cast interpolation returns an unexpected output

**URL:** <https://discourse.itk.org/t/ray-cast-interpolation-returns-an-unexpected-output/255>\
**Category:** Algorithms\
**Created:** [October 3, 2017, 10:45am UTC](https://discourse.itk.org/t/ray-cast-interpolation-returns-an-unexpected-output/255 "2017-10-03T10:45:55Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![christophe.vandijck](https://discourse.itk.org/user_avatar/discourse.itk.org/christophe.vandijck/32/593_2.png) [@christophe.vandijck](https://discourse.itk.org/u/christophe.vandijck)\
**Post date:** [October 3, 2017, 10:45am UTC](https://discourse.itk.org/t/ray-cast-interpolation-returns-an-unexpected-output/255/1 "2017-10-03T10:45:55Z")

</div>

Hi fellow ITK-users,

I’m trying to create several DRR images from a CT volume following several (oblique) directions. These DRRs often give unexpected results (empty images, the actual intensities not where I expected them in the image etc.).

I created a simple python script, similar to the cpp test function, to test a simple case: a gradient image where I evaluate the ray cast interpolator in one point. This point should result in a non-zero value, but it doesn’t. I checked with the original cpp test function, which remarkably also returns a zero-valued integration.

I added the python script below. Am I missing some crucial understanding on what the RayCastInterpolateImageFunction actually evaluates?

Thanks in advance,  
Christophe

```
import itk

def main():
    image_type = itk.Image[itk.UC, 3]

    image = image_type.New()
    start = (0, 0, 0)
    size = (30, 30, 30)

    region = itk.ImageRegion[3]()
    region.SetIndex(start)
    region.SetSize(size)
    image.SetRegions(region)
    image.Allocate()

    origin = (0, 0, 0)
    spacing = (1, 1, 1)

    image.SetOrigin(origin)
    image.SetSpacing(spacing)

    for slice in range(0, size[2]):
        for row in range(0, size[1]):
            for col in range(0, size[0]):
                val = slice+row+col
                image.SetPixel((col, row, slice), val)

    # writer = itk.ImageFileWriter[type(image)].New()
    # writer.SetFileName('test_image.mhd')
    # writer.SetInput(image)
    # writer.Update()

    ray_caster_type = itk.RayCastInterpolateImageFunction[image_type, itk.D]
    interp = ray_caster_type.New()

    interp.SetInputImage(image)
    print(interp)

    focus = (15, 15, 100)
    interp.SetFocalPoint(focus)

    transform_type = itk.TranslationTransform[itk.D, 3]
    transform = transform_type.New()
    interp.SetTransform(transform)

    interpolate_type = itk.LinearInterpolateImageFunction[image_type, itk.D]
    aux_interpolator = interpolate_type.New()
    interp.SetInterpolator(aux_interpolator)

    interp.SetThreshold(1.0)

    query = (15, 15, 15)
    integral = interp.Evaluate(query)
    print('Intergral = ', integral)

if __name__ == ' __main__':
    main()
```

---

<div class="post-metadata">

**Author:** ![matt.mccormick](https://discourse.itk.org/user_avatar/discourse.itk.org/matt.mccormick/32/7_2.png) [@matt.mccormick](https://discourse.itk.org/u/matt.mccormick)\
**Post date:** [October 5, 2017, 7:46pm UTC](https://discourse.itk.org/t/ray-cast-interpolation-returns-an-unexpected-output/255/2 "2017-10-05T19:46:52Z")

</div>

Hi Christophe,

Yes, your code is correct. After some investigation, the point `15,15,15` interrogates a position that is effected by a few outstanding bugs. This patch corrects the behavior:

[http://review.source.kitware.com/#/c/22683/](http://review.source.kitware.com/#/c/22683/)

Thanks,  
Matt

---

<div class="post-metadata">

**Author:** ![christophe.vandijck](https://discourse.itk.org/user_avatar/discourse.itk.org/christophe.vandijck/32/593_2.png) [@christophe.vandijck](https://discourse.itk.org/u/christophe.vandijck)\
**Post date:** [October 6, 2017, 8:38am UTC](https://discourse.itk.org/t/ray-cast-interpolation-returns-an-unexpected-output/255/3 "2017-10-06T08:38:10Z")

</div>

Hi Matt,

Thanks for the quick reply. Will check the patch out!

Christophe