Explanation of principal axis and principal moments

marabuuu · April 20, 2023, 10:41am

Dear all,

I am writing a blog post for biologists about different measurements including the ones provided by simpleITK and have problems to understand the measurements for principal axis and principal moments. I know that they are somehow used to compute a matrix, so are they maybe used to compute another measurement?

If so, could you please point me to these measurements and the formulas needed to compute them?

My aim is to explain these measurements in a way that it easy to understand.

Thank you in advance!

Best,
Mara

dzenanz · April 20, 2023, 1:23pm

Principal moments and principal axis code seems to be here:

github.com

InsightSoftwareConsortium/ITK/blob/815d0f0531eaff6a9ee18e3f1a2bbb2e251eff12/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx#L819-L891


      
          // Project the physical points onto principal axes
          VNLMatrixType transformedPixelLocations = principalAxesBasisMatrix * pixelLocations;
          
          
// find the bounds in the projected domain
          assert(transformedPixelLocations.columns() != 0);
          VNLVectorType minimumPrincipalAxis = transformedPixelLocations.get_column(0);
          VNLVectorType maximumPrincipalAxis = transformedPixelLocations.get_column(0);
          
          
for (unsigned int column = 1; column < transformedPixelLocations.columns(); ++column)
          {
            for (unsigned int i = 0; i < ImageDimension; ++i)
            {
              const double value = transformedPixelLocations(i, column);
              minimumPrincipalAxis[i] = std::min(minimumPrincipalAxis[i], value);
              maximumPrincipalAxis[i] = std::max(maximumPrincipalAxis[i], value);
            }
          }
          
          
// The minimumPrincipalAxis/maximumPrincipalAxis is from center of pixel to center of pixel
          // in the principal axis basis. The full extent of the pixels needs

This file has been truncated. show original

Some principal moments calculation is also here.

marabuuu · April 25, 2023, 12:37pm

Hi @dzenanz,

Thanks for the hint. It would be super cool if you or somebody else could explain this code to me as I have a hard time understanding it. I am interested in:

What are principal moments and axis?
When are these measurements of interest?
Are they maybe used to compute other measurement?

Best,
Mara

phcerdan · April 25, 2023, 12:45pm

Hi @marabuuu, who is interested in writing the blog post? I suggest you ask ChatGPT those questions, I am sure it can guide you!

Best,
Pablo

blowekamp · April 25, 2023, 4:57pm

There are two or three papers in the Insight Journal which describe the original implementation here is one of them:
https://www.insight-journal.org/browse/publication/176

I can’t seem to find the others at the moment.

Here is some additional code documentation about the computation:

github.com

InsightSoftwareConsortium/ITK/blob/815d0f0531eaff6a9ee18e3f1a2bbb2e251eff12/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx#L212-L231


      
          // moments computation
          //
          //  This computation has changed from what is documented in the
          //  original publication ( see class level documentation for
          //  reference).  It has been re derived to properly support the
          //  direction cosine matrix in the image.
          //
          // Using the same optimization of the computation and substitutions,
          // the new computation is derived with the following:
          //
          //   p_i = o_i + s_0*d_i_0 * x, where s_0 is spacing the line run, and
          //   d_i_0, it the components of the first column of the direction cosine
          //        matrix, and x is the index offset from o_i.
          //
          // Then the elements of the central moments _all_ are:
          //   S_i_j = sum_L_in_O( sum_p_in_L( p_i dot p_j ) )
          //
          // Then we follow the paper, in substituting p, expanding and
          // substituting for known summations over x. This is very similar to
          // equation 9 in the paper but with p_i dot p_j and NOT p_i dot p_i.

There are a couple aspects of what is implemented in ITK that makes it more complicated than basic index based moment calculations. First the computation take into considerations the image’s geometry which includes the pixel spaces and the direction cosine matrix. Second, it is optimized for the run-length encoding that the LabelObject used to store each label.

If you want to explain the computation is an easy way, the above optimization and spacial awareness should be omitted. And the wikipedia use for the standard index based computations as a reference.