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