Add a finetune step to the find_origin_slice method for the CatPhan to improve localisation accuracy

[keithoffer]

We've had some issues at Peter Mac where the origin slice of the CatPhan is detected slightly incorrectly. This can result in two symptoms:

1. The slice thickness is significantly under reported as the slice ramp isn't entirely in the detected origin slice
2. The other CatPhan modules aren't detected when offsetting, causing an error and PyLinac failing to analyse the scan.

I've played with a few solutions, and this was the one that seemed to work the best. It adds an extra 'fine tune' step at the end of the normal localisation of the image slice looking for the three markers embedded in the edge of the phantom. If it finds them in any slices near where it thought the center was, it updates the location of the origin slice by taking the median of the detected slices (same idea as currently done using the HU inserts). If none are found that meet the criterion, it defaults back to what it would have done before.

I couldn't work out how to download the run the unit tests for PyLinac, so I regression tested against ~600 CatPhan scans at Peter Mac across fanbeam and CBCT, CatPhan model 604 and 504. Overall this notably improved the detection across the set - you're welcome to some of these scans if you want.

It's not a free lunch, so it will make analysis take slightly longer - but hopefully that's a fair trade off to make the localisation more accurate and robust.

[crcrewso]

Thank you for this solution, I've been looking at other ones and havent liked the results as much! I'll test against our datasets, this might be the solution blocking our implimentation of this module in QA Track.

[crcrewso]

This should only plot when debugging,

[keithoffer]

Yep, looks like I accidentally pushed the version I was using when testing .... I'll try and fix this up.

[keithoffer]

Fixed

[jrkerns]

Thanks for the contribution Keith. I'll run this against the test suite. I removed it from the codebase because we are cloning the repo directly quite often now in our CI/CD pipelines and the download was taking far too long.

I can see why fine-tuning might be needed for point #2, but I can't follow #1. The ramp is quite long and an offset of a slice or two shouldn't affect the slice thickness. This makes me worry there is a bug or something.

[keithoffer]

I can send you some example image sets that show the issue, but when it happens you can see it in the final analysed image - sometimes the ramp is towards the edge of the ROI. It happens relatively rarely, but we hit it every now and then as we take so many CatPhan scans due to how many machines we have. Example picture below.

[jrkerns]

I see. That looks like it's several slices off. The HU-finding algorithm appears to be skewed here. I've not seen this phantom before. what model is it? Is that the 20cm housing around it?

[jrkerns]

All right, it's been a while since I looked at the 604 =). I will add some comments.

[jrkerns]

One dataset failed here because the HU values are quite skewed. Elekta in particular does not have very accurate HU values. I changed the threshold to np.percentile(slice.image.array, 98) and it worked better.

[keithoffer]

Sure, sounds reasonable. Had a quick test on some of my datasets and it seemed to still work like that so I'll change it. I've not got access to any Elekta sets, all my testing are Varian CBCT's or fanbeam scans on a Phillips scanner.

[jrkerns]

The 3 markers only appear to be for the 604. The other catphans have 0 or 1 AFAICT. This means this value should probably be a class attribute so in can be controlled per catphan.

[crcrewso]

There are similar markers in other newer CATPHAN phantoms, its just that they might be centered on other modules, some with a 4th marker on the underside. What if we had a 'refinable' property and a 'refine' method for phantoms where this technique is possible

[keithoffer]

Ah you're right. I asked if other models had three and was told they do, but looking at some old scans we have for the 504 you're right in that it only has one. I've not sure how robust finding only one peak would be though. I'll see if I can isolate my CatPhan 504 scans and see how it looks.

[jrkerns]

You shouldn't need this and below line as that check will be performed anyway below if you don't return here.

[keithoffer]

Yep, you're right. Adjusted.

[jrkerns]

Non-empty lists coerce to True so you can simply say if slices_where_markers_found:

[keithoffer]

Ah right. I normally prefer to manually check the length as I think the code is clearer, but I'll change it.

[jrkerns]

From what I can tell this appears to be unique to the 604. The 604 makes some of the HU ROIs longer than the others. I don't know why.

This appears to skew the HU detection center.

[crcrewso]

@jrkerns I don't think I accurately described this issue previously with issue #289 but this is exactly what I was seeing there too. I've used this patch to analyze a few of our outstanding test datasets and it in principle works, I think it would make more sense to set the defaults to one slice less than of the distance between any two modules.

[jrkerns]

Well, well. This is what happens when you don't have a good memory. Changing the defaults would work in theory. I can try that as an alternative and see if it's as stable as Keith's solution.

[jrkerns]

@crcrewso This may have been a 500 IQ call. I played around with the default for 604 and things seemed to work out nicely. @keithoffer Either 1) In the source code change this to ~4 https://github.com/jrkerns/pylinac/blob/master/pylinac/ct.py#L2081 or 2) Overload into a new class like so:

class MyCatphan604(CatPhan604):
    modules = {
        CTP404CP604: {"offset": 4},  # tweak this
        CTP486: {"offset": -80},
        CTP528CP604: {"offset": 42},
        CTP515: {"offset": -40},
    }

[crcrewso]

@jrkerns at this point I'm confused. I find @keithoffer's solution elegant for confirming the centre slice of the CTP404 module. Is there a technical or philosophical issue with migrating to this method to refine the location of the 0 slice? What am I missing?

[keithoffer]

I've updated the pull request based on a few of the suggestions.

If we can get a good result tweaking the offset then I think that's preferred, as that'd not require the extra complexity and time spent doing the fine tune step. However I do think the finetune will probably out perform it. I'll have play with our datasets and see how it looks. I'm not really too concerned what solution we go with as long as it works.

[jrkerns]

@crcrewso Elegance is in the eye of the beholder. A one-line change seems far more preferable than many. I think we can all agree that the 604 center-finding needs tweaking. I don't have as many datasets of the 604 but tweaking the default offset for the 604 404 module gets me within 1 slice of the center slice as given by Keith's solution. Considering that these modules are ~2cm wide or more, a single slice should not make a significant difference whereas 2-3 seems to be outside the comfort zone.

@keithoffer I would definitely be interested to see if there were a significant difference between an offset and the marker method. After adjusting by 4mm (setting the module offset to 4) the existing tests passed and the visual results seem pretty centered. I'd be happy to test a few of your scans to confirm. I'm not disagreeing with your solution, and it appears to work. I'd like to see if the simpler solution works first, and if not we can go up from there.

[crcrewso]

@jrkerns

a single slice should not make a significant difference

I think I can now see where our differing understandings are coming from. Since the Zero slice is used to determine the location of other modules, including the high resolution one, a single slice can and does make a huge difference to the analysis when trying to let a program like QATrack+ automatically analyze a dataset without direct interaction.

Additionally, by directly using internal markers to find the zero slice, its possible to analyze different slice thicknesses with the same settings, an additional requirement for our use cases.

Would you be open refactoring so that the center finding algorithm can be overridden on a per phantom basis?

[jrkerns]

@crcrewso I'm not against a refinement method, but the method is already overridable.

class MyCatPhan(CatPhan604):
     def find_origin_slice(self):
           default_hu_slice = super().find_origin_slice()
           # do whatever here to get a refined value
           return my_new_slice

Adding a refine method would literally look like this in the codebase:

class CatPhan604:
     def find_origin_slice(self):
           slice_num = super().find_origin_slice()
           return self.refine(slice_num)

    def refine(self, slice_num: int) -> int:
         """Default implementation where no refinement is done"""
          return slice_num

And then I would direct the user to override the refine method:

class MyCatphan(CatPhan604):
     def refine(self, slice_num):
           # do stuff here
           return my_refined_slice_num

It's the same difference for the user, just a different method to override.

by directly using internal markers to find the zero slice, its possible to analyze different slice thicknesses with the same settings

I'm not following how this is different.

a single slice can and does make a huge difference to the analysis

If you have a dataset demonstrating this I'm happy to add it to the test suite as I don't have too many. When you change the offset of one module it does not affect the other modules. The origin-finding algorithm is just a slice-finding algorithm; it just happens to usually correlate to the center of the 404 slice. Whatever slice it finds, each module is then offset by the amount (in mm) defined in the class or what the user has overridden. So changing the module offset, even of the CTP404 will not result in a different slice number for the 528. You can print out the <phantom>.ctp528.slice_num and change the offset for the 404 to see.

The default implementation for the 528 is to combine 3 slices together, taking the max value. This trades some potential accuracy for robustness (changing the slices to 1 from 3 resulted in going from 0.43lp/mm to 0.42 for me on the demo 604 dataset). I originally added it to mitigate pitch and yaw errors, but it also works for small longitudinal offsets.

The takeaway here is that besides me needing to change the 404 offset, if you believe the existing origin-finding implementation is consistent then use the offsets to refine the module positions. If you believe the existing implementation is inconsistent then override the find_origin_slice method.

[crcrewso]

@jrkerns

It's the same difference for the user

No it's not. I'm trying to use your data clinically, without relying on teaching others how to use Juniper or similar. I'm trying to maintain the bare minimum of out of band code in those systems so that upgrades can be done without my support.

The above refine methods look like a great solution, however I would like to see the default behaviour as, for any phantom with reliable definable fiducial markers near the zero slice, to refine off of them, and to wherever possible base offsets off of the manufacturers manual, in millimeters, not slices.

[jrkerns]

I'm still confused as to what you actually want. Can you write up a boilerplate code example of what you think the user should be doing when they want to refine the slice-finding algorithm?

[jrkerns]

@keithoffer Were you able to run scans on your local datasets?

[keithoffer]

Not yet, I've been occupied with commissioning AlignRT on our two linacs. I should have more time in the next two weeks to get back to this.

[jrkerns]

Hey, so we had a few RadMachine customers have 604 datasets that had slightly lower slice thickness than expected. After evaluating the datasets, this issue was the cause, specifically for 604 datasets (due to the longer HU plugs). I picked up this PR and while it mostly worked, there were two things I couldn't get around: a handful of clinics use CT sim BBs on their catphan for setup and they obviously appear quite bright and caused the resulting profile to sometimes be off by +/-90 degrees (the bbs were on the left and right sides). The second was that even though this worked for the 604, some other phantoms don't appear to be well-aligned to the markers. E.g.:

I don't know why; maybe it's an old design. Even though the 604s don't seem to be affected by this problem I also wanted something I could apply to other phantom designs.

I tried:

• A constant offset. This unfortunately is inconsistent from my own making. Using the RadMachine Jig for the catphan causes a few slices to get dropped from consideration for the HU module origin.
• This collapsed profile approach. Again, it did work for the vast majority of cases other than the ones I mentioned.
• Using the angles between the wires. This ended up working for all but one case of an extremely noisy and thin-slice scan.

For the sake of time for the RM customers who were having issues w/ this I opted for the last approach in the upcoming pylinac release. I did not get any meaningful differences in results other than the customers' slice thickness, which all improved (e.g. 1.5mm -> 1.9mm for a 2mm slice) due to not being at the edge of the wire section.

I'd like to leave this PR open in the event this latest approach does not work robustly. However, the new refine_origin_slice method is available and can be overloaded to do this or another approach as desired. The release notes also cover this in less detail: https://pylinac.readthedocs.io/en/latest/changelog.html#ct

Let me know if this works for your use cases.

[keithoffer]

Sorry for never getting back to this.

So if I run version 3.19.0 of PyLinac through my dataset looking at slice thickness I do get a definite improvement - but still some cases that look like incorrect central slice thickness detection to me. I'd suggest just closing this pull request as it was based on a flawed assumption (all models of the CatPhan have those pins) and the code would need to be entirely re-written anyway with the changes in 3.19.0.

If you're interested, I should be able to get permission to shoot you the set of images still failing for me - although the dataset is a little messy, a mix of different scan settings / different acquisition types / different machines.