computeCOR has strange behaviour

[AnderBiguri]

Works great sometimes, does't others.

[ghost]

Do you have any idea about why it doesn't work exactly?

I think the problem comes from getting center slice ( data=squeeze(data(slice,:,:)); technically it is not the center slice of your projection data). Then s2 and angles aux are getting the same matrix size of [x y] and on the other hand angle_grid, det_grid and test_data get [x 3*y], so interpolation doesn't work anymore. I think this only works if you use angle and projection data numbers can fit the matrix size. Also i experienced another unexpected error:

Error using griddedInterpolant The grid must be created from grid vectors which contain unique points.

Error in interp2>makegriddedinterp (line 229) F = griddedInterpolant(varargin{:});

Error in interp2 (line 137) F = makegriddedinterp(X, Y, V, method,extrap);

there shouldn't be any duplicated points in the input of interp2 function. If you let me know about what you think about this issue and overall problem of cor maybe i can help to correct it.

Thanks for the toolbox. Best regards.

[AnderBiguri]

@3elions I added the slice option for testing, the default slice is indeed the middle slice of the projections. Can you add a full example when you get that error? I do not get the error and my experiences are that when I run the code in simulated data, it detects the COR with quite a decent accuracy.

But when using real data, strange behavior arises, for example in one dataset the COR is correct only if the DSD/DSO multiplication in the end of the code is removed, which makes no sense.. I will go back to this next week but I am open to any discussion of this, as I really want to fix this.

[AnderBiguri]

@3elions Hi again. I have been looking at the interpolation you mentioned and I can not see the problems you describe. The input is indeed [x 3y], but there are no duplicate points, unless you have angles going over 2pi in length, but this would be (AFAIK) a quite unusual way of storing your data (and if so, you get a "strictly monotonically increasing" error) . s2 and aux are smaller matrix size, but that is OK, as they are the points where we want to sample the data, it does not need to be the same size as the data, it can be smaller without any problem, even just a single point.

I have tested several cases I can not reproduce the error you showed, could you please show me how to reproduce it?

And with the slice comment, did you meant that for an even sized detector, the center slice would be between pixels, i.e. we would need to interpolate to get it? that would make sense, I will try to add that to the code.

Please let me know if I am wrong in any of the things I said, I have been looking at this problem without finding a solution for too long and I may have missed an obvious step.

[ghost]

@AnderBiguri Hi,

Sorry for the late reply. I tried with my own binned scan data which is 256x256x400, and ndetector=256x256, ddetector=0.272mm. And the scan obtained from -pi to +pi,so angles cover 2pi with 0.0157 radian step size (matrix size is 1x400). So you should able to get the same error with this input. I really want to help because it would be very helpful insight if we can determine the central ray in cone beam.

For center slice, i added the images of center slices i got. So i would say the first one is not the center slice of projections. Maybe i am missing a point but in cone beam your projection data is (x,y,z). So i did't get how you get middle slice with squeeze on x domain. It should be done on z domain i think. Maybe i am wrong but in reference code the center slice is the very centre of projection slices. And also the difference between our inputs is might be the source of my problem with center slice.

data=squeeze(data(slice,:,:)); --- 256x400 data=squeeze(data(:,:,slice)); --- 256x256

I will do some more trials withdifferent data sets and code modifications during the week. I will keep you updated.

[AnderBiguri]

@3elions

Hum, not really, the data is stored as data(v,u,projection) (there is no x,y,z here, as its not the image, its the projection data that we need here). When we mean the center slice, its the row of the detector that is in the same plane as the source and the origin of the image (center point).

To ensure your projection data is in the correct order, visualize with plotProj(data,angles) and see if you see it moving left to right (as in the gif in the README in TIGRE Githubs page).

[ghost]

@AnderBiguri

I know, i meant (u,v,projection). Then i just misunderstood what you mean with center slice, because of middle slice definition. But the error is not related with that slice misunderstanding, i used your order for test. I tried with squeeze on u and v, and ended up with the same error. So i think because of angle array but these are correct angles.

[AnderBiguri]

@3elions Could you post your full geometry definition in TIGRE code so I can test it? I am unable to get that error.

[ghost]

@AnderBiguri

geo.DSD = 15.0044; % Distance Source Detector (mm) geo.DSO = 140.0079; % Distance Source Origin (mm) % Detector parameters geo.nDetector=[256; 256]; % number of pixels (px) geo.dDetector=[0.272; 0.272]; % size of each pixel (mm) geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector (mm) % Image parameters geo.nVoxel=[256;256;256]; % number of voxels (vx) geo.sVoxel=[7.479;7.479;7.479]; % total size of the image (mm) geo.dVoxel=geo.sVoxel./geo.nVoxel; % size of each voxel (mm) % Offsets geo.offOrigin =[0;0;0]; % Offset of image from origin (mm)
geo.offDetector=[0; 0]; % Offset of Detector (mm) % Auxiliary geo.accuracy=0.5; % Accuracy of FWD proj (vx/sample) % Load data and generate projections % define angles angles=deg2rad(-180:(360/(400-1)):180);

[AnderBiguri]

@3elions I assume you mixed DSD and DSO there.

I found why this is happening. Apparently the code does not accept repeated angles (not that these would be wrong to use in reconstruction). In your code, angles(1) and angles(end) are a different number, but they are the same angle (+pi and -pi), thus the code errors because it finds 2 equal angles, and The grid must be created from grid vectors which contain unique points. In your code, if you remove the first angle, or define the angles as angles=deg2rad(-180:(360/(400-1)):(180-360/(400-1))); this does not happen. I generally use angles=linspace(0,2*pi-2*pi/nangles,nangles)-pi;.

This said, I should add some error checking to computeCOR so it detects this and skips the duplicated angles.

[ghost]

@AnderBiguri Thanks for the help, this might the reason but i think your suggestion won't solve the problem, because all angles in the grid except first one are still repeating after all. When i tried with angles=deg2rad(-180:(360/(400-1)):(180-360/(400-1))); i got another interpolation related error:

The grid vectors do not define a grid of points that match the given values.

Error in interp2>makegriddedinterp (line 229) F = griddedInterpolant(varargin{:});

Error in interp2 (line 137) F = makegriddedinterp(X, Y, V, method,extrap);

Error in computeCOR (line 45) test = interp2(angle_grid, det_grid, test_data, angles_aux, s2, 'linear', 0);

It is probably very specifically related with my dataset than problem of code. I will try another projection dataset with non-repeating angles.

[AnderBiguri]

@3elions The following code works:

geo.DSO = 15.0044; % Distance Source Detector (mm)
geo.DSD = 140.0079; % Distance Source Origin (mm)
% Detector parameters
geo.nDetector=[256; 256];   % number of pixels (px)
geo.dDetector=[0.272; 0.272]; % size of each pixel (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector (mm)
% Image parameters
geo.nVoxel=[256;256;256]; % number of voxels (vx)
geo.sVoxel=[7.479;7.479;7.479]; % total size of the image (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel; % size of each voxel (mm)
% Offsets
geo.offOrigin =[0;0;0]; % Offset of image from origin (mm)
geo.offDetector=[0; 0]; % Offset of Detector (mm)
% Auxiliary
geo.accuracy=0.5; % Accuracy of FWD proj (vx/sample)
% Load data and generate projections
% define angles
angles=deg2rad(-180:(360/(400-1)):(180-360/(400-1)));

cube=headPhantom(geo.nVoxel.');
% generate projections
geo.COR=0.4;
tic
projections=Ax(cube,geo,angles,'interpolated');
toc

cor=computeCOR(projections,geo,angles);

And gives the correct COR (0.3976 when I run it).

You were perhaps still inputting 400 projections? Note that the change to the angles I made removes 1 of the projections and the angles are now 399. For real data, if you do indeed have a repeated angle, please instead of redefining the angles, just do not input the repeated angle and projection, for example by doing:

cor=computeCOR(projections(:,:,1:end-1),geo,angles(1:end-1));

[ghost]

@AnderBiguri i agree, for phantom it works perfectly, but for real data it doesn't. Maybe you could check it with my data? as well to check the error.

[AnderBiguri]

@3elions

Just tested your data. Few things here:

1) Your data is not in the right order. TIGRE needs proj(v,u,angles) and your data is proj(u,v,angles). You can see this because when you use plotProj(projection,angles), the rotation is not shown vertically (as in the gif in the main page of TIGRE), but horizontally instead. This will certainly make TIGRE not work. To fix your data, permute it: permute(projection,[2 1 3]). In the specific case of your data, this will show the image "on top", but ultimately only leads to having a reconstruction upside down, and its not a big deal. use flipud(permute(projection,[2 1 3])) if you want the "proper TIGRE way" way of having the data.

2) while this looks like micro-tomography, thus you will have COR errors, your data as a lot of mechanical vibration, and that will be a bigger problem than the COR. I assume this random wiggling is detector vibration. This is also very important to fix before the COR, as the COR calculation assumes that there is no random vibration, just the rotation axis is a bit shifted.

3) Finally, as I said in the last comment, if you are dealing with real data that has one angle repeated, use cor=computeCOR(projections(:,:,1:end-1),geo,angles(1:end-1));

for your code it gives me -0.0065. But as said, this doesn't fix much due to point 2).

[ghost]

@AnderBiguri Okay, thanks a lot for the assitance and comments, and sorry for wasting your time with wrong order. You are certainly right on point two but i have my own reconstruction program for detect and compansate this random wobbling motion during backprojection. So i think can correct wobbling motion in projection then calculate COR.

[AnderBiguri]

@3elions No worries at all, I am happy to help. If you know the values in mm (or whatever metric units you are using), you can add them to TIGRE in geo.offDetector, for reconstruction. Unfortunately this doesnt work (yet) in computeCOR , so its another thing that I need to add to this function! Alternatively you can always compute corrected projections using interpolation and go from there.