UCL / PETPVC

Partial Volume Correction in PET
Apache License 2.0
51 stars 14 forks source link

IterativeYang > 2 regions question #40

Closed ncalvertuk closed 5 years ago

ncalvertuk commented 5 years ago

Hi,

I have been attempting to implement the Iterative Yang method in MATLAB, and I have been unable to replicate some of the results in PETPVC.

I have been attempting to debug this myself, and I have edited some of the IY code in PETPVC to output some extra bits of information so that I can cross-reference with my MATLAB code.

What I have found is that, for the example I am currently looking at, my MATLAB code and PETPVC match up until I blur 'ImageYang' with the Gaussian psf (line 258 of petpvcIterativeYangPVCImageFilter.txx). One of my lines of debug code checks the sum of the image counts, and after the application of the Gaussian blur (in PETPVC) I find that the image counts has almost doubled. I have edited the code to write out this image so that I can compare to my MATLAB code. This is caused by the doubling of the counts in my background region (which has been segmented in the mask file) after the blur is applied. This makes no sense to me, as the blurring should preserve counts?

After doing a bit more testing, I found that this appears to only occur when I have more than 2 regions in my mask.

I was thinking this might be a problem with either my image or my mask (i.e. I've done something stupid!), however I am not sure why this would have this effect on the blurring stage of the algorithm (and not before).

Any help would be gratefully received!

Some extra info: I last cloned the PETPVC git repo on 2nd November, and my ITK version is 4.11.0. My image is a SPECT image (128 x 128 x 128 voxels, 4.42 mm^3 voxel size) of 6 spheres, which was converted from .dcm to .nii using MATLAB. My mask was also generated in MATLAB and is a .nii file with 7 regions (including the background). The mask voxel values are either 1 or 0 (i.e. not a probabilistic mask), and the sum along the 4th dimension is 1 for each voxel.

Thanks, Nick

bathomas commented 5 years ago

Hi,

Are you calling pvc_iy or pvc_diy?

ncalvertuk commented 5 years ago

Hi Ben,

I've been calling it like so:

./petpvc -i img.nii -m mask.nii -o PVC_Out.nii -pvc IY -d -x 15 -y 15 -z 15

-Nick

ncalvertuk commented 5 years ago

When using the pvc_diy function I am not seeing the same doubling of the counts when I blur ImageYang (line 222 of petpvcDiscreteIYPVCImageFilter.txx). It only seems to happen when I use the pvc_iy method (which I am using as I am interested in having a probabilistic mask).

bathomas commented 5 years ago

Strange! Can you post the debug output from:

./petpvc -i img.nii -m mask.nii -o PVC_Out.nii -pvc IY -d -x 15 -y 15 -z 15

Just to clarify, your mask image is 128 x 128 x 128 x 7, where each 3D volume is one region?

ncalvertuk commented 5 years ago

Hi Ben,

Yes, that's correct. The image is a set of 6 spheres, so each 128 x 128 x 128 x 1 is a sphere region and the final one is the background region.

The debug output, including my new lines is as follows:

./../pvc_iy Scan1.nii ROIs.nii PVC_Out.nii -d -x 15 -y 15 -z 15 -i 1 -d Start fuzziness calculation sum in region 1 = 6 sum in region 2 = 8 sum in region 3 = 23 sum in region 4 = 49 sum in region 5 = 92 sum in region 6 = 184 sum in region 7 = 2.09679e+06 matrix: 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1

Iteration:
1 Summed imageExtractedRegion[0] = 118 SummedCounts[0] = 118

Sum = 118 , Mean = 19.6667 , Size = 6 Summed imageExtractedRegion[1] = 60 SummedCounts[1] = 60

Sum = 60 , Mean = 7.5 , Size = 8 Summed imageExtractedRegion[2] = 1028 SummedCounts[2] = 1028

Sum = 1028 , Mean = 44.6957 , Size = 23 Summed imageExtractedRegion[3] = 1279 SummedCounts[3] = 1279

Sum = 1279 , Mean = 26.102 , Size = 49 Summed imageExtractedRegion[4] = 1214 SummedCounts[4] = 1214

Sum = 1214 , Mean = 13.1957 , Size = 92 Summed imageExtractedRegion[5] = 2927 SummedCounts[5] = 2927

Sum = 2927 , Mean = 15.9076 , Size = 184 Summed imageExtractedRegion[6] = 1.35071e+07 SummedCounts[6] = 1.35071e+07

Sum = 1.35071e+07 , Mean = 6.44179 , Size = 2.09679e+06 19.6667 7.5 44.6957 26.102 13.1957 15.9076 6.44179 vecRegMeansCurrent: 19.6667 7.5 44.6957 26.102 13.1957 15.9076 6.44179 sum of ImageYang: 1.35137e+07 PSF = [40.5758, 40.5758, 40.5758] sum of blurred ImageYang: 2.70208e+07 sum of IY./conv(IY,gauss): 2.09714e+06 size of IY./conv(IY,gauss): [128, 128, 128]

bathomas commented 5 years ago

There is another application called pvc_simulate which just applies the smoothing to an image. Can you try smoothing Scan1.nii with this and seeing if you get the same problem.

Is the image matrix cropped/zoomed tightly around the spheres? I'm wondering if it could be boundary problems, or as the PSF is quite large whether the Discrete Gaussian filter is somehow failing.

ncalvertuk commented 5 years ago

I don't think I am getting the same problem with the pvc_simulate. The image is 128 x 128 x 128 and the voxels are 4.42 mm so the FOV is quite large around the spheres which take up a relatively small area in the centre of the image.

I have found that it appears to affect the final mask region only. When I switched the order of the mask file so the background is the first region and one of the spheres is the last (i.e. 7th) region, then it is this sphere that doubles.

I am currently updating my version of ITK (4.11 to 4.13) and I will also do a fresh install of PETPVC once this is done to see if there is some issue. I'll let you know how I get on!

Nick

ncalvertuk commented 5 years ago

Updating ITK & reinstall PETPVC seemed to solve the problem. I guess something went wrong during one of the installations that I did not see...

Thanks! Nick

bathomas commented 5 years ago

Thanks. I would hope that one of the iterative Yang tests would have failed during make test, but I guess we'll never know.

Ben