- load the image: . The image will be displayed correctly (but reinit is necessary)
- right click on image -> reinit
Revisions and Commits
|Restricted Differential Revision|
|Resolved||kislinsk||T25682 mitkBaseGeometry.cpp: Internal ITK matrix inversion error, cannot proceed.|
|Open||None||T26165 How to handle erroneous DICOMS/input data?|
- Mentioned In
- T27033: Clean up stale remote branches
rMITK47af5550e874: Check if matrix is rotation. Safe calculate dominant axes.
rMITK6b4ea15d8155: Check if matrix is rotation. Safe calculate dominant axes.
- Mentioned Here
- T22254: Render windows show *not* the corresponding orientation after reinit with permuted axes
rMITK1ca1a7b0574e: Compose sliced geometries of planes with matching normal
Log output with debug version:
!5.547! WARNING: EnsurePerpendicularNormal(): Plane geometry was _/askew/_, so here it gets rectified by substituting the 3rd column of the indexToWorldMatrix with an appropriate normal vector: [ 0.921861 -0.918454 -0.746066 ], original vector zed was: [ -0.943643 -3.64367e-09 -1.16599 ]. MitkWorkbench: ../../../../MITK/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp:273: virtual void mitk::SlicedGeometry3D::InitializePlanes(const mitk::BaseGeometry*, mitk::PlaneGeometry::PlaneOrientation, bool, bool, bool): Assertion `(standardPlaneNormal - (top ? 1.0 : -1.0) * worldPlaneNormal).GetSquaredNorm() < 0.000001' failed.
@full @schererj @kislinsk : how did you convert these files exactly? If I load the DICOM with a current MITK Workbench -> no problem, re-init fine. Save as nifti, quit, restart, load nifti, re-init -> no problem as well.
So we should investigate two things: 1. what is different in your conversion process? 2. are "your" nifti files somehow invalid and should be rejected on load or should we able to handle them as well.
I converted them with a simple two-liner: mitk::IOUtil::Load<mitk::Image>(); mitk::IOUtil::Save();
Peter and I checked some output NIFTI images and they are valid 3d+t images. Reinit did also work. So for now I suspect that there is another processing step involved? The images above do not match my output files regarding name and file extension I think.
I have the issue again with other data data from a cooperation partner. Also Slicer has issues with the original DICOM data set:
Geometric issues were found with 1 of the series. Please use caution. reg inside examine Loading with imageIOName: GDCM Window/level found in DICOM tags (center=212.0, width=487.0) has been applied to volume 25: CINE_segmented_SAX_InlineVF Irregular volume geometry detected, but maximum error is within tolerance (maximum error of 4.00539e-05 mm, tolerance threshold is 0.001 mm). Loading with imageIOName: GDCM Window/level found in DICOM tags (center=212.0, width=487.0) has been applied to volume 25: CINE_segmented_SAX_InlineVF_1 Irregular volume geometry detected, but maximum error is within tolerance (maximum error of 4.00539e-05 mm, tolerance threshold is 0.001 mm).
Looks like the Slicer warning is important, but still it succeeds in loading the data. The numerical error is rather low. Found this discussion which also links to the python code Slicer is using and which emits the warning. Could be of some use for further investigations. https://discourse.slicer.org/t/dicom-scalar-volume-load-irregular-geometry-warning-overly-stringent/3761/14
had a further look into it:
- looked into the matrix vectors. Looks like the first and the second vector is normalized, but the third is not (https://academo.org/demos/3d-vector-plotter/)
- tried to transpose the matrix
- tried to change the last column by sign
- the rectification leads to a normal vector that is the same like the first vector.
- we can set the matrix manually to space directions: (1.9,0,0) (0,1.9,0) (0,0,12) (pixel spacing). But that's not a real solution as we usually can't alter the data.
- newer data tends to work better
- newer in terms of date
- newer in terms of acquisition machine software version
- working and not working data stems from the same machine (Siemens)
The meaning of the dominant axes is also documented here:
It is basically a mapping from the axes of the image coordinate axes to world coordinate axes. It determines which axis of the image coordinate system is the "axial" axis, which one is the "sagittal" and which is the "coronal". Since the "index-to-world" transformation is not just a permutation of the indices and a scaling, but you can rotate the volume in any way, it is not always obvious which axis is which. That's why the "max" call.
If you get the same index twice, that's certainly wrong. It can happen if the plane is rotated by *exactly* 45 degrees around one or two world axes. This case is 'undecided', one choice is not better than the other. You should, however, assign different axes to the different dimensions as 'dominant'. When looking for the dominant axis, you should also check if the value has already been "taken".
Hope that helps.
Does anyone know why the inverse is always calculated and not simply set to the transpose? After all, we are talking about rotation matrices.
In this context, the EnsurePerpendicularNormal method (https://phabricator.mitk.org/source/mitk/browse/master/Modules/Core/src/DataManagement/mitkPlaneGeometry.cpp$37) seems rather complex and also does not catch all cases since it only checks if axis 3 is perpendicular to 1&2 but not if 1&2 are perpendicular to each other. I would prefer to simply check if the matrix R is a rotation matrix by checking
- det(R) == 1?
- R*R^T = ID?
I'm also not sure if it is fine to simply substitute one vector if this is not the case. It does not work in all cases and we also do not want to change the data, right? It is further more unclear why this is checked only when calculating the standard planes.
I would propose to check if R is a valid rotation matrix immediately when loading an image and output a warning if this is not the case. For example in cases of minor deviations from condition 2., as it is the case in our test data, it should be up to the user to ignore this. It's a question of how you want to set your epsilon.
Such an "imperfect" rotation matrix in conjunction with "almost" 45° rotations can of course yield the above mentioned error during the dominant axis calculation, but this can easily be caught. I will implement a proposal.
While I changed the matrix check, I left it at the current position in in the PlaneGeometry.
We should discuss if it is sensible to move this ckec and also it's applications to a different location, e.g. to SetIndexToWorldTransform of BaseGeometry.
As I see, the inverse matrix calculation came in with this commit of mine:
which refers to this wiki page from Slicer as 'inspiration':
But that page only determines the dominant axes per orientation and does not actually uses inverse matrix calculation.
Most likely I took that code from NifTK that has a viewer that can display several images side-by-side and in image space rather than world space.
Now I'm curious if it would crash with the attached image.