Empty volume after resampling, struggling with the origin

Beginner Questions
1

Hi!
I have some problems with resampling the volume and determining the new origin.
I have a volume (57x57x100) and then I permute the axis,change the direction and set the spacing using:

img = sitk.ReadImage(nifti_path)    
img = sitk.PermuteAxes(img, [1,2,0]) 
rai_direction = (1,0,0,0,-1,0,0,0,-1)
img.SetDirection(rai_direction)
spacing = (res,res,res)
img.SetSpacing(spacing)

I want to rotate around the top middle part of the volume and I also want to translate this point. Therefore, I make sure that this point is at 0,0,0 by using:

origin = img.GetOrigin()
out_origin = origin[0]-(img.GetSize()[0] * img.GetSpacing()[0])/2+0.5*img.GetSpacing()[0], origin[1], origin[2]+(img.GetSize()[2] * img.GetSpacing()[2])/2-0.5*img.GetSpacing()[0]
img.SetOrigin(out_origin)    
rot_center = (0,0,0)

Then I want to apply an translation and rotation:

transform = sitk.VersorRigid3DTransform()
 translation = sitk.TranslationTransform(3, (pos[0][0]*1000, pos[0][1]*1000, pos[0][2]*1000))
transform.SetTranslation(translation.GetOffset())
 rotation = sitk.VersorTransform([orientation[0][0], orientation[0][1], orientation[0][2], orientation[0][3]], rot_center)
transform.SetRotation(rotation.GetVersor())
transform.SetCenter(rotation.GetCenter())

Then I want to resample the volume:

size = img.GetSize()
    bounds = list()
    for x in [0, size[0]-1]:
        for y in [0, size[1]-1]:
            for z in [0, size[2]-1]:
                bounds.append(img.TransformIndexToPhysicalPoint((x, y, z)))

    trans_bounds = list()
    for b in bounds:
        trans_bounds.append(transform.TransformPoint((b[0], b[1], b[2])))
    trans_bounds = np.array(trans_bounds)

    xmin = np.min(trans_bounds.T[0])
    xmax = np.max(trans_bounds.T[0])

    ymin = np.min(trans_bounds.T[1])
    ymax = np.max(trans_bounds.T[1])

    zmin = np.min(trans_bounds.T[2])
    zmax = np.max(trans_bounds.T[2])

    output_origin = (xmax, ymax, zmin)
    res_all = (0.5,0.5,0.5)
    res = res_all[0]
    output_size = np.array([int(np.round((xmax - xmin) / res)), int(np.round((ymax - ymin) / res)), int(np.round((zmax - zmin) / res))])
    output_size = output_size.astype(int)

    resampleFilter = sitk.ResampleImageFilter()
    resampleFilter.SetTransform(transform.GetInverse())
    resampleFilter.SetInterpolator(sitk.sitkLinear)

    resampleFilter.SetSize(output_size.tolist())
    resampleFilter.SetOutputOrigin(output_origin)
    resampleFilter.SetOutputSpacing(res_all)
    resampleFilter.SetOutputDirection(img.GetDirection())
    resampleFilter.SetOutputPixelType(sitk.sitkUInt8)
    resampleFilter.SetDefaultPixelValue(0.0)
    registeredImg = resampleFilter.Execute(img)

However, the resampling parts give me some problems. Because I change the direction form (1 0 0 0 1 0 0 0 1) to (1 0 0 0 -1 0 0 0 -1) and then do the resampling as stated above, I got an empty volume. I think I should change ‘output_origin’ in the resampling part so that the volume does not fall out of the view. But what should this be? I don’t understand how this works. Or is there another way? Hope you can help me.

Kind regards,
Renate

2

You could use Slicer to explore how resampling works. These two modules would be of interest:

Slicer is good for experimentation. Once you have established a procedure to accomplish what you want, it is then easier to write code to accomplish it.

1 Like