Skip to content

Add notes on the order of matrix elements of Direction and GridDirection parameters#1208

Merged
N-Dekker merged 2 commits intomainfrom
DOC-order-Direction-matrix-elements
Aug 16, 2024
Merged

Add notes on the order of matrix elements of Direction and GridDirection parameters#1208
N-Dekker merged 2 commits intomainfrom
DOC-order-Direction-matrix-elements

Conversation

@N-Dekker
Copy link
Copy Markdown
Member

@N-Dekker N-Dekker commented Jul 24, 2024

Triggered by the discussion "Elastix idiosyncracy?", started by Soren Christensen, July 11, 2024, at https://discourse.itk.org/t/elastix-idiosyncracy/7049

@N-Dekker
Copy link
Copy Markdown
Member Author

For the record, elastix::ResamplerBase::ReadFromFile() did already read the Direction matrix in column-major order with commit 695f922, "ENH: write result images with correct direction cosines (as in fixed image). Also store direction cosines in TransformParameters file", 19 June 2008:

for (unsigned int j = 0; j < ImageDimension; j++ )
{
this->m_Configuration->ReadParameter( direction(j,i), "Direction", i*ImageDimension+j );
}

@N-Dekker
Copy link
Copy Markdown
Member Author

Just double-checked: a direction matrix for 3D images specified by elastix parameter (Direction v0 v1 v2 v3 v4 v5 v6 v7 v8) will internally be stored in the data buffer of a 3x3 itk::Matrix as { v0, v3, v6, v1, v4, v7, v2, v5, v8 }.

* \transformparameter Direction: The direction cosines matrix of the fixed image
* that was used during registration, and which is used for resampling the deformed moving image
* if the UseDirectionCosines parameter is set to "true".\n
* Please note: The matrix elements should be specified in column-major order!\n
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would describe more explicitly that this differs from ITK conventions:

Suggested change
* Please note: The matrix elements should be specified in column-major order!\n
* Please note: The matrix elements are specified in column-major order (Fortran-style).
* This is not the same as ITK's convention of using row-major order (C-style) for flattening
* direction matrix of an image into a list.\n

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think @lassoan's comments reduces the room for error even further!

Copy link
Copy Markdown
Member Author

@N-Dekker N-Dekker Aug 8, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lassoan Thank you for your suggestion 👍 I'm fine with your suggestion, at least if it's the only place where it should the documented. (If this specification should be documented on multiple places, I'd rather elaborate on the meaning of "column-major order" just once, rather than on each of those places. But this appears the one and only place where it needs to be documented, right?)

One nit-pick question, why do you suggest replacing

The matrix elements should be specified in column-major order!

with

The matrix elements are specified in column-major order!
?

I wanted to express that the user should specify the matrix elements in that order.

Copy link
Copy Markdown
Contributor

@lassoan lassoan Aug 8, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, this detailed note should only appear once in the Elastix documentation. So, if this needs to be mentioned elsewhere in the future than that place should refer to this note or this note could be moved out elsewhere.

To me, "should" has some degree of uncertainty. If "should" is needed to make it clear that it is up to the user has to specify the values then we could use "must" or "shall" to make it clear that column-major order is not just a recommendation but a requirement.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, thanks, then shall we make it "must"?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lassoan @sorenchr2011 Stefan (@stefanklein) just pointer out to me that MetaIO also uses column-major order for its TransformMatrix representation. See https://discourse.itk.org/t/metaio-image-rotation-conventions/4525 and https://itk.org/Wiki/ITK/MetaIO/Documentation#Associated_transformations

So then, is ITK's convention really using row-major order? Shouldn't the comment say something like:

Please note: The matrix elements are specified in column-major order (Fortran-style, like MetaIO's matrix representation).
ITK often uses row-major order (C-style) for flattening direction matrix of an image into a list.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lassoan @sorenchr2011 Stefan (@stefanklein) just pointer out to me that MetaIO also uses column-major order for its TransformMatrix representation. See https://discourse.itk.org/t/metaio-image-rotation-conventions/4525 and https://itk.org/Wiki/ITK/MetaIO/Documentation#Associated_transformations

So then, is ITK's convention really using row-major order?

ITK internally uses row-major order. Various file formats don't have to use the same conventions as ITK. ITK-based libraries (such as Elastix) don't have to use the same conventions as ITK either, but when a user gets a vector from an ITK method then he expects to be able to use it as input in another method in an ITK-based library without any conversion, so it is worth trying to harmonize conventions. Bringing in MetaIO into this discussion is unnecessary and may be distracting.

The matrix elements are specified in column-major order (Fortran-style, like MetaIO's matrix representation).
ITK often uses row-major order (C-style) for flattening direction matrix of an image into a list.

I don't think this is accurate, because ITK is consistent in how it stores matrices in memory and how they may be exposed as flat arrays via Python wrapping.

Anyway, any documentation would improve things and I don't really mind the details.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Anyway, any documentation would improve things and I don't really mind the details.

@lassoan OK, thanks! We obviously agree that the documentation should specify the order of the matrix elements.

Stefan (@stefanklein) further clarified their choice for column-major order (back in 2008) by saying that as each column represents one vector of direction cosines, it may be more readable to have the direction cosines of a vector together, in the flattered representation.

@N-Dekker N-Dekker force-pushed the DOC-order-Direction-matrix-elements branch from 448a43f to 16added Compare August 9, 2024 10:27
@N-Dekker
Copy link
Copy Markdown
Member Author

N-Dekker commented Aug 9, 2024

@lassoan I added you as co-author, please check!

@N-Dekker N-Dekker force-pushed the DOC-order-Direction-matrix-elements branch from 16added to 1771017 Compare August 9, 2024 10:29
@N-Dekker
Copy link
Copy Markdown
Member Author

N-Dekker commented Aug 9, 2024

For the record, my little test to observe the ordering of matrix elements:

static constexpr unsigned int                       ImageDimension = 3;
itk::Matrix<double, ImageDimension, ImageDimension> direction;

// Create a "flat" array of matrix elements, {1, 2, ..., N}
std::array<double, ImageDimension * ImageDimension> flatArray;
std::iota(flatArray.begin(), flatArray.end(), 1.0);

// Retrieve the elements from the flat array by a nested for-loop, just like elastix does at:
// https://github.com/SuperElastix/elastix/blob/995f471aeb17929f8c0557135997e78f54326390/Core/Kernel/elxElastixTemplate.hxx#L1081-L1090
for (unsigned int i = 0; i < ImageDimension; ++i)
{
  for (unsigned int j = 0; j < ImageDimension; ++j)
  {
    direction(j, i) = flatArray[i * ImageDimension + j];
  }
}

for (const double x : direction.GetVnlMatrix())
{
  std::cout << " v" << x;
}

std::cout << std::endl;

Based on

for (unsigned int i = 0; i < FixedDimension; ++i)
{
for (unsigned int j = 0; j < FixedDimension; ++j)
{
if (!configuration.ReadParameter(directionRead(j, i), "Direction", i * FixedDimension + j, false))
{
return false;
}
}
}

Output:

v1 v4 v7 v2 v5 v8 v3 v6 v9

Similar nested for loops are in ResamplerBase::ReadFromFile() and SimilarityTransformElastix::ReadCenterOfRotationIndex(InputPointType &).

@N-Dekker
Copy link
Copy Markdown
Member Author

N-Dekker commented Aug 9, 2024

I still wonder, shouldn't "Direction" be documented in "elxResamplerBase.h", if it's actually being retrieved by ResamplerBase::ReadFromFile()...? (It's being used as OutputDirection of its ResampleImageFilter.) But then again, it's also being used for the dummyImage in SimilarityTransformElastix::ReadCenterOfRotationIndex, and for ElastixTemplate::GetOriginalFixedImageDirection.

@N-Dekker

This comment was marked as resolved.

@N-Dekker N-Dekker force-pushed the DOC-order-Direction-matrix-elements branch from 1771017 to 329558e Compare August 9, 2024 13:41
@N-Dekker N-Dekker marked this pull request as ready for review August 9, 2024 14:38
@N-Dekker N-Dekker force-pushed the DOC-order-Direction-matrix-elements branch 2 times, most recently from be28a43 to 6c58973 Compare August 12, 2024 14:25
@lassoan
Copy link
Copy Markdown
Contributor

lassoan commented Aug 12, 2024

I still wonder, shouldn't "Direction" be documented in "elxResamplerBase.h"

Maybe the best could be describe the convention in every method where the user can get or set a direction array, but keep it very short and simple, with a link to more details. Something like "direction: ..., matrix elements in column-major order (conversion may be needed, see details at URL)" .

* This is not the same as ITK's convention of using row-major order (C-style) for flattening
* direction matrix of an image into a list.\n
* example: <tt>(Direction -1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.1)</tt>\n
* Default: identity matrix. Elements are sorted as follows: [ d11 d21 d31 d12 d22 d32 d13 d23 d33] (in 3D).
Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see now, here it already says:

Elements are sorted as follows: [ d11 d21 d31 d12 d22 d32 d13 d23 d33] (in 3D)

@N-Dekker N-Dekker force-pushed the DOC-order-Direction-matrix-elements branch from 6c58973 to c2a6a96 Compare August 15, 2024 09:18
N-Dekker and others added 2 commits August 15, 2024 11:39
Triggered by the discussion "Elastix idiosyncracy?", started by Soren Christensen, July 11, 2024, at https://discourse.itk.org/t/elastix-idiosyncracy/7049

Co-authored-by: Andras Lasso <lasso@queensu.ca>
Triggered by the discussion "Elastix idiosyncracy?", started by Soren Christensen, July 11, 2024, at https://discourse.itk.org/t/elastix-idiosyncracy/7049
@N-Dekker N-Dekker force-pushed the DOC-order-Direction-matrix-elements branch from c2a6a96 to c304fb0 Compare August 15, 2024 09:41
@N-Dekker N-Dekker changed the title DOC: Add a note on the order of matrix elements of Direction parameter Add notes on the order of matrix elements of Direction and GridDirection parameters Aug 15, 2024
@N-Dekker
Copy link
Copy Markdown
Member Author

Going to merge now, but maybe the documentation on Direction and GridDirection should be re-organized somehow, eventually. (Direction is now documented in elastix::TransformBase while it's not specifically used by that base class. On the other hand, GridDirection is documented multiple times, redundantly, for four b-spline transforms.)

@N-Dekker N-Dekker merged commit 7a7d9be into main Aug 16, 2024
@N-Dekker N-Dekker deleted the DOC-order-Direction-matrix-elements branch August 16, 2024 10:11
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants