SuperElastix / SimpleElastix

Multi-lingual medical image registration library
http://simpleelastix.github.io
Apache License 2.0
513 stars 149 forks source link

Problems with resampling prior to registration #377

Open ctorti opened 4 years ago

ctorti commented 4 years ago

Hi all,

After getting tripped up when trying to register a 3D CT image to a 3D MR image (Issue #370), I decided to give myself a better chance at success by first trying to succeed at registering two 3D MR images of similar (but different) resolutions and voxel sizes.

The ultimate aim is to transform a set of contours that belong to one MR image to another. The contour points are in a DICOM-RTSTRUCT file. The images themselves are DICOMs.

Following is some info about the two images:

FixedIm Dtype = 32-bit float MovingIm Dtype = 32-bit float

FixedIm Size = (212, 256, 30) MovingIm Size = (192, 256, 35)

FixedIm Spacing = (0.8984375, 0.8984375, 5.0) MovingIm Spacing = (0.9375, 0.9375, 5.0)

FixedIm Origin = (-104.92378234863281, -152.4906463623047, -9.22148609161377) MovingIm Origin = (-115.14457702636719, -124.83148193359375, -10.8151273727417)

I managed to succeed in registering the images as follows:

ElastixImFilt = sitk.ElastixImageFilter() ElastixImFilt.SetFixedImage(FixedIm) ElastixImFilt.SetMovingImage(MovingIm) ElastixImFilt.SetParameterMap(sitk.GetDefaultParameterMap('affine')) ElastixImFilt.LogToConsoleOn() ElastixImFilt.LogToFileOn() RegIm = ElastixImFilt.Execute()

The results seem fine:

FixedIm Dtype = 32-bit float MovingIm Dtype = 32-bit float RegIm Dtype = 32-bit float

FixedIm Size = (212, 256, 30) MovingIm Size = (192, 256, 35) RegIm Size = (212, 256, 30)

FixedIm Spacing = (0.8984375, 0.8984375, 5.0) MovingIm Spacing = (0.9375, 0.9375, 5.0) RegIm Spacing = (0.8984375, 0.8984375, 5.0)

FixedIm Origin = (-104.92378234863281, -152.4906463623047, -9.22148609161377) MovingIm Origin = (-115.14457702636719, -124.83148193359375, -10.8151273727417) RegIm Origin = (-104.92378234863281, -152.4906463623047, -9.22148609161377)

FixedIm Min, Max = 0.0, 2015.0 MovingIm Min, Max = 0.0, 1870.0 RegIm Min, Max = -136.6031036376953, 1769.6104736328125

Using matplotlib I plotted FixedIm, MovingIm and RegIm using an interactive plot with sliders that I got from (I beleive) SimpleITK (see attached). The registration seemed to work as expected. I even overlaid the contours onto FixedIm and RegIm and all seems well:

20200515_MovingIm_reg_to_FixedIm_mid-stacks

What's more tricky however, is to plot the ROIs over MovingIm. I attempted to do so as follows:

I read the DICOM-RTSTRUCT file and created a txt file containing the contour points in the format required for elastix, i.e.

point number_of_points (or) number_of_indeces point1_x point1_y point1_z point2_x point2_y point2_z ...

Then I transformed the contour points:

TransformixImFilt = sitk.TransformixImageFilter() TransformixImFilt.SetTransformParameterMap(ElastixImFilt.GetTransformParameterMap()) TransformixImFilt.SetFixedPointSetFileName(FixedContourPtsFname) # Need to explicitely tell elastix that the image is 3D since by default it will # only transform to 2D: TransformixImFilt.SetMovingImage(MovingIm) TransformixImFilt.Execute()

elastix automatically created a text file called "outputpoints.txt" - oddly with a completely different format/structure than the structure that was required of the input points...

I parsed the generated txt file and generated an array of points with the following structure (called TransformedContourPts):

[ [], [], [], [ [x1, y1, z1], [x2, y2, z2], ... ] ... ]

where each row corresponds to a slice in the 3D image. If the item is [] it simply means that the corresponding slice did not have any contours (and hence, no contour points). For any slice that does have a contour, there is an array of x, y, z points that make up that contour.

The above structure is the same that I used for FixedIm (called FixedContourPts). That way, when matplotlib plots slice i, it plots the points from row i in FixedContourPts.

But then the next issue hit me - the slice locations (z positions) for MovingIm are not the same as the slice locations for FixedIm. So aside from getting some practice in creating the point set, transforming them and parsing the output (which will hopefully come in handy soon enough), it was a path to nowhere.

Then another idea hit me: If I resample MovingIm using FixedIm as the reference image, then I can overlay the points in TransformedContourPts) onto ResampledMovingIm, since the resampling process aligns the slices in space.

Following is how I repeated the registration but this time using a resampled MovingIm:

ElastixImFilt = sitk.ElastixImageFilter() ElastixImFilt.SetFixedImage(FixedIm) # Resample MovingIm: Resampler = sitk.ResampleImageFilter() Resampler.SetReferenceImage(FixedIm) ResampledMovingIm = Resampler.Execute(MovingIm) ElastixImFilt.SetMovingImage(ResampledMovingIm) ElastixImFilt.SetParameterMap(sitk.GetDefaultParameterMap('affine')) ElastixImFilt.LogToConsoleOn() ElastixImFilt.LogToFileOn() ElastixImFilt.Execute() RegIm = ElastixImFilt.GetResultImage()

Following are a couple of screenshots - 2D slices of FixedIm, ResampledMovingIm and RegIm at slice number 14 within the stacks:

20200515_ResampledMovingIm_reg_to_FixedIm_slice14

and at slice number 7 within the stacks:

20200515_ResampledMovingIm_reg_to_FixedIm_slice7 Slice number 7.

Looking at the slices it's clear that something has gone wrong with ResampledMovingIm. Yet despite this the registration looks fine...

Does anyone have any ideas as to whether this is expected or not, and what I can do to remedy it?

If anyone has any general suggestions on how to best or most easily achieve the main goal of transforming contour points using the same transformation applied to register two 3D images, that would be much appreciated.

Thanks, Cristiano

kaspermarstal commented 4 years ago

Hi Cristiano

Thank you for your detailed explanation. What columns do you parse from outputpoints.txt? There should be a column with transformed indices (not points in world coordinates) that can you use to create your TransformedContourPts structure if I understand your problem correctly. Looks like you need the z image index to insert the contour at the correct slice.

ctorti commented 4 years ago

Hi Kasper,

Thanks for your support. I was only parsing the OutputPoint values. I've not been able to make much sense of the index values (InputIndex, OutputIndexFixed, OutputIndexMoving) since the x and y indeces exceed the indeces of the actual images. E.g.:

FixedIm.GetSize() = (212, 256, 30) MovingIm.GetSize() = (192, 256, 35)

yet

Point 0 ; InputIndex = [ 252 359 14 ] ; InputPoint = [ 126.344037 171.830275 -54.501356 ] ; OutputIndexFixed = [ 248 353 26 ] ; OutputPoint = [ 126.638100 184.489069 3.193988 ] ; Deformation = [ 0.294064 12.658794 57.695343 ] ; OutputIndexMoving = [ 220 356 4 ]

i.e. the x and y indices in InputIndex > the x and y dims of MovingIm.

Why would that be the case?

Also, much in the same way that I can't make sense of the z coordinates in OutputPoint (I'm getting values ~3 mm when the z coord of InputPoint is ~ -54 mm, and I know that there is only a z-displacement on the order of ~1 mm), similarly, I can't make any sense of the z indeces (I'm getting values ~4, whereas the z index of InputIndex is ~14, and I know that the indeces should be about a slice apart or so).

If I understood you correctly, it looks like you've suggested to use the z index from OutputIndexMoving to indicate which slice the point given by OutputPoint belongs.

The problem is, as indicated above, the difference between the z indeces of the input and output differ by ~ 10, so the point gets allocated to a completely different depth within the image.

I wondered if perhaps I could specify the fixed image for Transformix (as well as the moving image), i.e.:

TransformixImFilt.SetFixedImage(FixedIm)

but that command was not allowed.

I've attached the inputpoints.txt and outputpoints.txt files in case they are of any help.

inputpoints.txt outputpoints.txt

Taking a logical approach I see the following possibilities:

  1. I have not properly formatted inputpoints.txt for Elastix
  2. I have not properly used Transformix
  3. The transform parameter map is incorrect
  4. Transformix has not properly transformed the points for some reason
  5. I have misinterpreted what the points and indeces in outputpoints.txt correspond to

Re. 1: Since the format required was quite simple I think it's unlikely, but I've attached the file in the hope that if it is incorrect, someone might be able to highlight the error.

Re. 2: The steps used were copied from elsewhere so I think this is unlikely that I've not used Transformix correctly.

Re. 3: I've attached the parameter map here in the hope that if there is an issue it will be obvious to someone.

TransformParameters.0.txt

Re. 4: I'll assume that Transformix works properly when fed the correct inputs.

Which leads to what I suspect to be the more likely possibility: That I've misinterpreted what the points and indices in outputpoints.txt mean.

I've searched the documentation for SimpleElastix and Elastix and did a web search but haven't been able to find any guidance on this. Any advice would be much appreciated.