How should one set the output origin when resampling an image?

1

Hello,

I would like to know how to set the output origin when resampling a volume to a different spacing. I understand from Fundamental Concepts — SimpleITK 2.4.0 documentation that in SITK the origin is in the middle of the first (corner) voxel, not in the corner of the first corner voxel. I just want to make sure the origin of the resampled volume should be new_origin = original_origin - original_spacing/2 + new_spacing/2. Like that it would be in the center of the new corner voxel.

Can you confirm or explain why it should not be the case?
Does this function look like a correct wat of resampling to you?

def resample_volume(
        volume, interpolator=sitk.sitkBSpline, new_spacing=[1.0, 1.0, 1.0]
    ):
        original_spacing = volume.GetSpacing()
        original_size = volume.GetSize()
        new_size = [
            int(round(osz * ospc / nspc))
            for osz, ospc, nspc in zip(original_size, original_spacing, new_spacing)
        ]
        return sitk.Resample(
            image1=volume,
            size=new_size,
            transform=sitk.Transform(),
            interpolator=interpolator,
            outputOrigin=(
                    # Shifting the origin from center of olf voxel to corner of volume
                    np.array(volume.GetOrigin())-np.array(original_spacing)/2
                    # Shifting from corner to center of the new voxel of a new dimension
                    )+np.array(new_spacing)/2,
            outputSpacing=new_spacing,
            outputDirection=volume.GetDirection(),
            defaultPixelValue=0,
            outputPixelType=volume.GetPixelID(),
        )

Thank you.
Best regards,
Sebastien

2

Hello @sebquet,

It looks like you want to keep the same origin as your original image, so no need to do any computation, just set the origin outputOrigin = volume.GetOrigin(). For a detailed discussion of resampling, see the Transforms_and_Resampling notebook in the toolkit’s Jupyter notebooks repository.

If you want to create an isotropic (same sized voxels) volume, see the Results_Visualization notebook, make_isotropic function. That function also supports “standardizing the axes”, so that the resulting volume has an identity cosine matrix.

3

Hello,

Your proposed method does not take into account the DirectionCosine matrix of the images.

In the following code the continuous index is computed and than it’s transformed into the physical coordinates:

3 Likes
4

Hello @zivy,

I may have omitted to say that I want to retain the image original location.

Hello @blowekamp,

Thank you so much for providing this solution which is more robust that my proposed solution. This is exactly what I was looking for.

Thank you.

2 Likes
5

For reference, this question was also asked and discussed on the Slicer forum:

In short, I raised the issue that by shifting the origin outside of the original image (outside the bounding box of voxel centers - using a code like @blowekamp provided above) the image is not just interpolated but voxel values near the image boundary are computed by extrapolation. Extrapolation may introduce additional artifacts, therefore the image origin should only be shifted after carefully considering of all the consequences.

1 Like