Crop patch out of an image

Beginner Questions
,
1

I’m trying to use itk (not sitk) to crop out a patch of another image. I have two images A and B that are overlapping, but don’t share the same size, spacing or image origin.

Image A is basically a 224x224 square patch. I’d like to crop out the same patch of image B. I’ll then resample the patch from image B to also make it 224x224.

I tried to do it like that, but it seems like I have misunderstood the transformation to physical points:

def crop_xray(xray, target):
    target_region = target.GetLargestPossibleRegion()
    physical_size = target.TransformIndexToPhysicalPoint(
        target_region.GetSize()
    )
    physcial_index = target.TransformIndexToPhysicalPoint(
        target_region.GetIndex()
    )

    xray_size = xray.TransformPhysicalPointToIndex(physical_size)
    xray_index = xray.TransformPhysicalPointToIndex(physcial_index)

    # Extract the patch from image A
    xray_region = itk.ImageRegion[3](xray_index, xray_size)
    extractFilter = itk.ExtractImageFilter.New(Input=xray)
    extractFilter.SetExtractionRegion(xray_region)
    extractedPatch = extractFilter.GetOutput()
    extractedPatch.Update()
    return extractedPatch
2

You should divide target_region.GetSize() by target.GetSpacing() to get physical_size. Then multiply by xray.GetSpacing() to get xray_size.

The rest of the code look OK. Why not use RegionOfInterest filter?

3

Thank you for the fast reply. I see why i would need to divide the size by the spacing. I just assumed that’s what all the TransformPhysical… functions do. I guess I didn’t get the concept.

The xray_size makes sense now. The patch seems to have a size of [36.37523560968537, 36.37523560968537, 1.0] which corresponds to what I can measure in slicer. However, I still get a type error despite casting this size to int:

_itkImageRegionPython.itkImageRegion3_swiginit(self, _itkImageRegionPython.new_itkImageRegion3(*args))
                                                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: Expecting a sequence of int (or long)

If I understood you correctly my xray_index is correct. I converted that to int too.

RegionOfInterest sounds like a better approach. I’ve tried to use it for the last 2 hours now but it always failed in various places. I’m still a bit overwhelmed by ITK. Do I have to learn C++ to better understand the documentation? There are no examples in python for certain functions and I struggle to use them correctly.

4

Does this example help with syntax?

5

Actually

smallSize = itk.Size[Dimension]()
smallSize.Fill(10)

index = itk.Index[Dimension]()
index.Fill(0)

confused me even more. I didn’t get index in my example either. I see how my patch is 36.37 by 36.37 wide. But I didn’t get what the starting point for my patch is. In the mean time I tried to used target.GetOrigin() as a starting point to my patch which I then convert to the xray image. However, I now get the error that my patch is out of the image.

I’ll try to study your example more. Maybe it starts to make sense. But thank you anyway.

6

How about:

patch_size = itk.Size[Dimension]()
for d in range(Dimension):
  patch_size[d] = xray_size[d]

and then you can use patch_size for region size.

7

I think I found an/the issue. In some cases the patch to crop is not entirely inside the xray. My workaround will be to check that before, then apply the code as suggested by you. In the end I probably have to pad the image to 224x224 again.

Or is there a smarter work around for that?

8

Why not use resample filter? It does what you want. Use your region of interest to construct the target image. And you could pass a boundary condition to resample filter.

2 Likes
9

Perfect that could help a lot. Thank you very much for your help.

10

I’m sorry I still can’t get it to work. I tried to incorporate the suggestions in the last couple dates. I still keep getting the error “Requested region is (at least partially) outside the largest possible region.”

Is there anymore hints you could give me? I hoped the resample filter will get rid of this.

def crop_xray(image_B, image_A):
    # Get the starting point and size of Image A
    start_point = np.asarray(image_A.GetOrigin(), dtype=int)
    size = np.array(image_A.GetLargestPossibleRegion().GetSize())

    # Convert the starting point and size of Image A to points in Image B
    end_point = start_point + (size - 1)

    start_point_in_B = image_A.TransformIndexToPhysicalPoint(
        start_point.tolist()
    )
    end_point_in_B = image_A.TransformIndexToPhysicalPoint(end_point.tolist())

    region = itk.ImageRegion[3](
        np.array(
            np.asarray(start_point_in_B) * image_B.GetSpacing(), dtype=int
        ).tolist(),
        np.array(
            np.asarray(end_point_in_B) * image_B.GetSpacing(), dtype=int
        ).tolist(),
    )

    # Use the starting and end points in Image B to extract the corresponding region from Image B
    image_B.SetRequestedRegion(region)
    crop_filter = itk.RegionOfInterestImageFilter[image_B, image_B].New()
    crop_filter.SetInput(image_B)
    crop_filter.SetRegionOfInterest(region)
    interpolator = itk.LinearInterpolateImageFunction.New(image_B)
    resampled = itk.resample_image_filter(
        crop_filter,
        interpolator=interpolator,
        size=[224, 224, 1],
        output_spacing=image_A.GetSpacing(),
        output_origin=image_A.GetOrigin(),
    )
    return resampled
11

I suggested bypassing region of interest filter entirely, and the parameters you are giving to RegionOfInterestImageFilter should be given to resample filter directly.

image_B.SetRequestedRegion(region) is the most likely cause of “Requested region is (at least partially) outside the largest possible region.”

1 Like
12

That was it! Works now :slight_smile: Thanks a lot.

1 Like
13

Maybe share the working code, preferably the full version, in case someone else has the same question.

1 Like
14

Sure:

def crop_xray(image_B, image_A):
    # Get the starting point and size of Image A
    start_point = np.asarray(image_A.GetOrigin(), dtype=int)
    size = np.array(image_A.GetLargestPossibleRegion().GetSize())

    # Convert the starting point and size of Image A to points in Image B
    end_point = start_point + (size - 1)

    start_point_in_B = image_A.TransformIndexToPhysicalPoint(
        start_point.tolist()
    )
    end_point_in_B = image_A.TransformIndexToPhysicalPoint(end_point.tolist())

    region = itk.ImageRegion[3](
        np.array(
            np.asarray(start_point_in_B) * image_B.GetSpacing(), dtype=int
        ).tolist(),
        np.array(
            np.asarray(end_point_in_B) * image_B.GetSpacing(), dtype=int
        ).tolist(),
    )

    # Use the starting and end points in Image B to extract the corresponding region from Image B
    image_B.SetRequestedRegion(region)
    interpolator = itk.LinearInterpolateImageFunction.New(image_B)
    resampled = itk.resample_image_filter(
        image_B,
        interpolator=interpolator,
        size=[224, 224, 1],
        output_spacing=image_A.GetSpacing(),
        output_origin=image_A.GetOrigin(),
    )
    return resampled
1 Like