RTKConsortium / RTK

Reconstruction Toolkit
Apache License 2.0
244 stars 144 forks source link

Add a 'geometry tolerance' variable for pre-compiled applications #504

Closed MarkGardnerUSyd closed 2 years ago

MarkGardnerUSyd commented 2 years ago

I create geometry files using Matlab and then run the pre-compiled rtk applications (rtkfdk.exe, rtkadmmtotalvariation.exe etc) on a windows terminal (yes I know it's weird). The use of two different operating systems causes small variations in the geometry file (), giving me errors such as:

File: C:\RTK-2.3.0\RTK-2.3.0\src\rtkThreeDCircularProjectionGeometryXMLFileReader.cxx
Line: 162
Description: ITK ERROR: Matrix and parameters are not consistent.
Read matrix from geometry file: 
-1237.5 -0 862.666 -160000
0 -1500 -0 0
0.48163 -0 0.876374 -1000

Computed matrix from parameters:
-1237.5 0 862.667 -160000
0 -1500 0 0
0.481631 0 0.876374 -1000

I believe this error is displayed when the difference between the two matrices is larger than a certain value, however this value is smaller than the computation error that occurs between different operating systems.

If you were to add an optional variable in the compiled applications that specifies the maximum allowed difference between the two geometry files, then I could increase this value so that this error no longer appears. An example could be something like:

rtkfdk --geometryTolerance e-5

Hope this makes sense.

MichaelColonel commented 2 years ago

It could be easily done. In example of yours the inconsistency is about 0.001 (value in first row, third column) i would say it's quite significant!

SimonRit commented 2 years ago

Thanks for the report. Can you provide the parameter values which were used to build this projection matrix? It would be good to understand the root cause of this behavior difference between OS.

MarkGardnerUSyd commented 2 years ago

We have a Matlab Script that creates the geometry file:

 StartAngle = -180;           %First projection angle
 EndAngle = 180;
 ProjNum = 895;
 angs = linspace(StartAngle,EndAngle,ProjNum) - 90;
 x_offset = -160;

 for jj = 1:numel(angs)
    G(:,jj) = [angs(jj);0;0;0;0;x_offset;0;SID;SDD];
    Pmats(:,:,jj) = gtop(G(:,jj));
 end

 WriteRTKGeometryXML(Pmats,geoFile);

where the function gtop is:

function [P] = gtop(g)
%Takes physically meaningful parameters g to make projection matrix P
%% the P part

GantryAngle = g(1);OutOfPlaneAngle=g(2);InPlaneAngle=g(3);SourceOffsetX=g(4);SourceOffsetY=g(5);DetOffsetX=g(6);DetOffsetY=g(7);SID=g(8);SDD=g(9);

GA = GantryAngle*((2*pi)/(360));                 
OPA = OutOfPlaneAngle*((2*pi)/(360));              
IPA = InPlaneAngle*((2*pi)/(360)); 
SOX = SourceOffsetX;
SOY = SourceOffsetY;
DOX = DetOffsetX;
DOY = DetOffsetY;

%% 
% R1 = [cos(-IPA),-sin(-IPA),0,0;sin(-IPA),cos(-IPA),0,0;0,0,1,0;0,0,0,1];
% R2 = [1,0,0,0;0,cos(-OPA),-sin(-OPA),0;0,sin(-OPA),cos(-OPA),0;0,0,0,1];
% R3 = [cos(-GA),0,sin(-GA),0;0,1,0,0;-sin(-GA),0,cos(-GA),0;0,0,0,1];
% 
% S1 = [1,0,SOX-DOX;0,1,SOY-DOY;0,0,1];
% S2 = [-SID,0,0,0;0,-SID,0,0;0,0,1,-SID];
% S3 = [1,0,0,-SOX;0,1,0,-SOY;0,0,1,0;0,0,0,1];
%
% P = S1*S2*S3*R1*R2*R3;

P = [   [ - sin(GA)*(cos(OPA)*(DOX - SOX) + SDD*sin(IPA)*sin(OPA)) - SDD*cos(GA)*cos(IPA), sin(OPA)*(DOX - SOX) - SDD*cos(OPA)*sin(IPA),   SDD*cos(IPA)*sin(GA) - cos(GA)*(cos(OPA)*(DOX - SOX) + SDD*sin(IPA)*sin(OPA)), SDD*SOX + SID*(DOX - SOX)];...
        [   SDD*cos(GA)*sin(IPA) - sin(GA)*(cos(OPA)*(DOY - SOY) + SDD*cos(IPA)*sin(OPA)), sin(OPA)*(DOY - SOY) - SDD*cos(IPA)*cos(OPA), - cos(GA)*(cos(OPA)*(DOY - SOY) + SDD*cos(IPA)*sin(OPA)) - SDD*sin(GA)*sin(IPA), SDD*SOY + SID*(DOY - SOY)];...
        [                                                                cos(OPA)*sin(GA),                                    -sin(OPA),                                                                cos(GA)*cos(OPA),                      -SID]];

and where WriteRTKGeometryXML prints the geometry matrices to the xml file

MichaelColonel commented 2 years ago

Could it be that the digits precision in Matlab has low value? It looks like the precision is 6, according to your matrix element data.

SimonRit commented 2 years ago

I agree with @MichaelColonel, this seems to be a problem in your WriteRTKGeometryXML, not an OS-specific issue with RTK.

SimonRit commented 2 years ago

BTW, it's possible to use RTK's python code from Matlab if you need to write an XML file from Matlab. I can share an example if you're interested.

MarkGardnerUSyd commented 2 years ago

If you could share an example of using RTK's python code from Matlab to write an XML file from Matlab that would be very helpful. Thanks

SimonRit commented 2 years ago

Here is an example (hope it works, I cannot test it right now)

itk = py.importlib.import_module('itk');
geometry = itk.RTK.ThreeDCircularProjectionGeometry.New();                                                                                                                                                                               
geometry.AddProjection(geometry, 1000., 1536., 28.)
xmlWriter = itk.RTK.ThreeDCircularProjectionGeometryXMLFileWriter.New();
xmlWriter.SetObject(xmlWriter,geometry);
xmlWriter.SetFilename(xmlWriter,'geometry.xml');
xmlWriter.WriteFile(xmlWriter);
SimonRit commented 2 years ago

Any update on this issue?

MarkGardnerUSyd commented 2 years ago

Sorry for the delay. Had to get the itk-rtk python module working.

The attached code works well. Thanks for helping with this issue.

Mark

SimonRit commented 2 years ago

Thanks for the feedback. I don't see any code attached, can you copy paste it please? Can you close the issue if you consider it solved?

MarkGardnerUSyd commented 2 years ago

Sorry I meant that the code you attached above works.