[Insight-developers] HDF5 Transform I/O Again -- definitions of fixed/moving parameters

[Johnson, Hans J]

Brian and others,

1) Several conversations are going on in e-mail with different
participants getting partial information.  Would you please consolidate
all these comments into a wiki page, or at least one document?


Let me be clear here:  I actually would like to have more generalized
transforms (at least up to affine transforms) natively supported in ITK
images themselves, and sincerely hope that your team finds a way to
implement this successfully in ITKv4.  I've several places where this
feature would simplify my book-keeping tremendously and it would improve
both efficiency and reduce smoothing.

=== Historical reference, NOT an endorsement of what we historical did ===
 Also, let us remember that historically IO has been a set of convenience
external modules to the core ITK framework.  In fact originally ITK
explicitly excluded IO from the toolkit in part as an effort to avoid
contaminating the core ITK framework with compromises from the largely
inconsistent IO conventions.  The addition of IO to ITK has been a
tremendous burden on the development team in trying to accommodate the
various interpretations and irreconcilable conventions adopted by the
community of users (be glad that we are not addressing the 100 different
flavors of the Analyze file format, it has been a nightmare to deal with).
 Right or wrong, historically when the team I work with comes across one
of these conflicts of interpretation, we try to match a "what would DICOM
do" least common denominator approach approach.

The NIFTI-1 IO context is a great one to use for discussion purposes, but
lets also consider the implications of all the other formats as well.  In
the very near future the NIFTI-2 format will be alive (probably this
summer) and that may add new wrinkles to this discussion.   The comments
below are based on what we INTENDED the IO mechanisms to do.  This was not
always adhered to or completely implemented.  I think that in the context
of ITKv4 we should focus on what we should be doing, and not on one
specific mis-implementation.  Most of the coding changes being discussed
will not be difficult to implement, but I want to think carefully about if
they are going to be impossible to support.

See more exhaustive (exhausting :) comments below.

Hans

On 4/2/11 7:46 AM, "brian avants" <stnava at gmail.com> wrote:

>hans
>
>> This makes sense to me, but ITK images do not have the means to
>>interpret
>> a general affine transform in the way described below.  It does not fit
>> the current ITK model (it does fit the "SPM" model of images and fits
>> nicely with matlab).
>
>was it jim or you that said the sform should translate to an itk
>transform ...
>perhaps that is the best/easiest route .... give nifti i/o a
>sform=>transform
>function.

It could be done, but that would change the fundamental behavior of ITK
images in the entire toolkit and for all file format types.  This is a big
conceptual change even if it is a potentially small coding change with no
performance penalties up to an affine matrix defintion.

Given the scenario where a nifti image with an generalized affine
transform is read in, what is the correct behavior for writing it to all
the other medical image formats that only support rigid transforms?

I am being VERY cautious here because I helped go through the incredibly
painful process of adding image orientation (I.e. The direction matrix),
and that process took more than 5 years and there are still occasional
bugs found related to adding that notion.  Admittedly this change *should*
not be as painful.   One scary part that we may encounter again is that
originally algorithms were written without rotation taken into account,
and when rotation was added derivatives and grid spacings were not/are not
being properly computed (I.e. Bspline Initializer is still not quite right
even after many years and several attempts).  How are we going to assure
that algorithmic optimizations were not made under the assumption of rigid
registrations (I think this is a very infrequent, but difficult task to
identify)?


>
>> Proposal:
>> 0) Under the premise that failure is better than ambiguity.
>
>agreed.
>
>> 1) Images that are read by ITK should be able to be written to the same
>> format without loss of basic physical space relationships as represented
>> by the ITK toolkit.
>
>ok --- even if it means a loss of information b/c itk does not
>represent something correctly,
>e.g. intent codes, sforms, etc?

Here I think we'll dis-agree a bit (which is not a bad thing).  We need to
think about the union of file format intent codes (NRRD, NIFTI, Analyze,
MetaIO, DICOM) and realize that it quickly becomes an intractable problem
to consistently address all possible information consistently in ITK
proper. I believe that IO should be decoupled from ITK proper, and that
where possible IO should map components to the proper itk::Image ivars,
and everything else gets shoved into the meta-data dictionary and becomes
the application programmers responsibility.


The crux of the problem is that the "correct representation" is often
problem/application/lab/programmer specific.  We should provide tools and
policies for helping the application developer deal with this, but within
limited resources and  influence I don't believe that it is possible to
resolve all the problems.

>
>i'd weigh the long-term impact on people who are new to itk (and
>maintainability) vs short-term difficulty of adding these
>representations.

I think we also need to weigh the overall maintainability of ITK for the
next 10 years.  It is my belief that consistent and pervasive
implementation of something as simple as "intent code" would be a MAJOR
(yes all caps MAJOR) development effort to implement in a way that would
work correct without explicit problem/application/lab/programmer
intervention.

>
>> 2) ITK only writes images in the qform method, as that most closely
>>aligns
>> with the dicom standard, of which the ITK image physical space also
>> conforms.
>
>ok.  is the idea to avoid resampling?
Resampling was never really part of the thought process, resampling is an
algorithmic consideration and we were only concerned with data
representation with respect to IO.  The situation was that given several
options of how to do this, and understanding that ITK images can ONLY have
rigid registrations, and that we are lazy programmers and there was
documentation of how to do this.  Ultimately the qform only definition
seemed like the perfect fit because it un-abmiguously represents an ITK
physical space completely for both writing and reading without any extra
computations needed decompose the sform matrix or to verify that it is
rigid conforming.

>
>> 3) Therefore we should change the default IO factory methods of ITK
>> behavior so that only images with only the qform matrix are read.
>>Reading
>> any other nifti image will require special handling by the application
>> developer.  We SHOULD provide examples and helper functionality to
>> accomplish this, but that should not be the default behavior of the IO
>> mechanism.
>>
>
>this is certainly a minimal requirement.   right now, these problems
>are mostly silent
>and thus frustrating.
I was actually hoping you would complain more loudly about this.  ITKv4 is
our only opportunity for the next 10 years to change this.  I agree 100%
about the silent failures, but we (mostly I) made early mistakes and had
to maintain backwards compatibility.  What ever the final decision is, I'd
like to implement item #0 more rigorously.

>
>> One of the proposals was to use only the sform, and abandon the qform.
>> This was investigated once, but ITK images can not support the general
>> affine matrix for physical space definitions.
>>
>
>you wrote previously:
>
> "In the end we choose to only interpret one, and somewhat arbitrarily
>choose that if both exists then only the qform would be respected.
>The sform only is used if the qform is ascent."
>
>compare to chris's comment:
>
>" I would personally suggest that by default you use the s-form when it is
>available, and only use the q-form if the s-form is not available. "
>
>and also
>
>"If the qform code is "NIFTI_XFORM_SCANNER_ANAT" and the sform code is
>"NIFTI_XFORM_ALIGNED_ANAT", "NIFTI_XFORM_TALAIRACH", or
>"NIFTI_XFORM_MNI_152", then you should certainly use the latter. This is
>also a consensus of all the other NIfTI implementations.
>
>I think this would be easy to implement, even if you include quats to
>represent the data internally. You can use mat44_to_quatern to get the
>sform to look like a qform."
>
>clearly this is a readily available solution in itk --- and all we
>lose is the shearing/scaling of the affine sform.  really, the change
>would be the prioritization.
>
>just food for thought --- i will send you some comments from phil on
>other aspects of nii+itk.

This raises two areas of ambiguity:
1) What do you do if the s-form and q-form are both rigid, but do not have
consistent physical space representations.
2) What if the q-form has a rigid (ITK compliant) and the s-form has a
generalized affine.  What is the correct behavior?


I agree that we should pay more close attention to to the qform and sform
codes.

ITKv3 images have only one data representation, and that is consistent
with "NIFTI_XFORM_SCANNER_ANAT".  The other 3 codes are not compatible
with ITKv3 images and would require the use of an ITKv3 Transform outside
the image to properly represent it.

In ITKv3 we SHOULD throw a very loud exception if only the sform code is
present and  any of "NIFTI_XFORM_ALIGNED_ANAT", "NIFTI_XFORM_TALAIRACH", or
"NIFTI_XFORM_MNI_152" because they do no map to the ITK data model.

In ITKv3 the way this could be conceptually handled is to use the image IO
mechanism to read in the image and qform, and then have a transform IO
mechanism read in the same file and return an ITK affine matrix that was
the the "sform with the qform rotation removed"


I think it is a VERY VERY VERY bad idea to loose the sheering/scaling of
the affine sform matrix as a default behavior.  It would be much better to
fail (unless of course the application programmer explicitly request this
behavior).


===========================================================================
===============

So.... Enough historical ranting.... Let me through out some ideas for
ITKv4.


Let's pursue the idea of having ITK image data representation support an
affine matrix that defines the voxel-lattice to physical space mapping.

Currently the ITK image representation has the following public interface
to ivars for supporting this mapping (3D for this example):

O=Origin  -- a 3D point
D=Direction -- A rigid orthonormal direction cosign matrix that is
represented as a 3x3 matrix, but always assumed to be rigid currently.
S=Spacing  --- Stored as an array of 3 values, but interpreted as a
diagonal matrix=
P=Computed Physical space
I=VoxelIndex -- a vector that may contain continuous index points into the
voxel lattice

P=D*S*I+O

There are also private cache variables that precompute the "D*S" matrix --
but hopefully we don't have to touch that.

##1   The problem becomes, given an general affine matrix, can it be
uniquely decomposed into a D, S, O components in such a way that all
algorithms that currently ask for myImage->GetSpacing(),
myImage->GetOrigin(), myImage->GetDirection() continue to have compatible
physical space usage?

It may be as simple as adding one new member function to itk::Image
/**
* The use of this function will decompose the 3x3Matrix, and make calls to
* "SetSpacing, SetOrigin, SetDirection" with appropriate decomposed values.
*/
SetGeneralizedAffinePhysicalSpaceRepresentation(3x3Matrix,Origin)
{
Extract D
Extract S
Extract O

this->SetDirection(D);
this->SetOrigin(O);
this->SetSpacing(S);
}

SetGeneralizedAffinePhysicalSpaceRepresentation(4x4Matrix)
{
Test that no projects are present, else fail.
this->SetGeneralizedAffinePhysicalSpaceRepresentation(4x4Matrix->GetRotatio
n(), 4x4Matrix->GetTranslation());
}


##2   Given that D may now NOT contain an orthonormal matrix, what do we
do when a storage request to DICOM is made (where only 6 parameters are
set and the other 3 are implied based on the orthonormal assumption).   I
would recommend loud failure here.

Hans
















________________________________
Notice: This UI Health Care e-mail (including attachments) is covered by the Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is confidential and may be legally privileged.  If you are not the intended recipient, you are hereby notified that any retention, dissemination, distribution, or copying of this communication is strictly prohibited.  Please reply to the sender that you have received the message in error, then delete it.  Thank you.
________________________________