[Insight-developers] Coordinate system semantics of ImageIOBase?

[Luis Ibanez]

Hi Steve,

Please find answers below.

   Regards,


       Luis


--------------------
Steve M. Robbins wrote:
> On Mon, Dec 29, 2008 at 12:10:22PM -0500, Luis Ibanez wrote:
> 
> 
>>Since you list is long and very detailed, please find my
>>comments below, interleaved with your text.
> 
> 
> Thanks for your detailed responses.  You have answered
> most of my questions, but a few questions were not
> asked quite right so I have some follow-on clarifying
> questions, below.
> 
> 
> 
>>>To do so successfully, I need to understand the semantics of
>>>its base class, ImageIOBase.  I'm getting tangled up in
>>>coordinate system issues.
>>>
>>>For example, "GetDimensions" is documented as follows
>>>
>>>  /** Set/Get the image dimensions in the x, y, z, etc. directions.
>>>   * GetDimensions() is typically used after reading the data; the
>>>   * SetDimensions() is used prior to writing the data. */
>>>
>>>... but I think the intent is that the dimension are NOT
>>>world-space x, y, and z, but rather column, row, and slice.
>>>Is that correct?
>>>
>>
>>
>>  GetDimensions() and SetDimensions() relate to the number of Pixel
>>  along each of the directions of space.
> 
> 
> What I meant to ask is: "what space?".  In your response, you are
> calling the axes X, Y, and Z.  But this is an "image space" distinct
> from the "world space", right?  This is the space that some people
> would say uses IJK coordinates?
> 
> Furthermore, in this "image space" isn't it true that the pixel
> spacings are precisely 1.  The spacing values in ImageIOBase
> are for the world space only, right?
> 
> 


      * They are in the IJK coordinates.
      * The values correspond to "number of pixels".
      * Spacing is irrelevant here.


> 
>>>What should one assume about the world space coordinates?  I believe
>>>it is safe to assume a right-handed coordinate system.  Do we assume
>>>they should be LPS (as in DICOM)?  If so, do we convert incoming
>>>coordinates (MINC uses the RAS convention)?  I would assume not as
>>>I've read about folks using ITK for non-medical images like in
>>>astronomy.  But if no conversion is done, how does an application
>>>handle loading from multiple formats; e.g. DICOM and MINC?
>>>
>>
>>
>>      This is a delicate point.
>>
>>      The standard in ITK is to use LPS.
> 
> 
> So what happens if your image is not of a person?  Is the stance of
> ITK that:
> 
>   1. ITK is only used for medical data, so there is always a
>      patient-centric world reference frame.
> or
> 
>   2. ITK specifically uses LPS for medical data.  Other image sources
>      are free to choose their own world reference frame.
> 
> Sorry to belabour the point, but I am serious about augmenting the
> documentation of this class and I want to be both specific and
> accurate.  If the case is 2, are there any other reference systems
> that have been settled on?  For example, I thought someone was using
> ITK for astronomy.
> 


      When the image is not a person,
      we worry much less about this.

      In that case,
      your application is free to choose
      any interpretation for the Direction matrix.

      A point to keep in mind is that coordinate systems
      are always *relative to* something.

      What to use as the reference is a choice that must
      be made at the application level.

      ITK simply provides the placeholders for carrying
      the information of position and orientation relative
      to the "reference" of your choice.


> 
> 
>>>The documentation for Read() is as follows
>>>
>>>  /** Reads the data from disk into the memory buffer provided. */
>>>  virtual void Read(void* buffer) = 0;
>>>
>>>Do we dump the bits into buffer exactly as they appear on disk with no
>>>type conversion?  We're supposed to pay attention to the
>>>ImageIORegion, right?  So the buffer is guaranteed to be large enough
>>>to hold the region last set with SetIORegion() ?
>>>
>>
>>       Yes, you must pay attention to (e.g. Respect) the ImageIORegion.
>>
>>       In the Read() method, you should simply read the pixel
>>       data the way it is stored in the image file.
>>
>>       Pixel type conversion will be managed for you by the
>>       ImageFileReader. Please see lines: 401-417 of
>>       Insight/Code/IO/itkImageFileReader.txx
> 
> 
> Understood.  But this gets in the way of pixel intensity scaling ...
> 


          Yes, you must consider both issues.
          This is not a trivial matter, since the order
          of operations here is important.

          You may have to apply rescalings first, using
          a pixel type 'float' or 'double' as intermediary,
          and then cast the resulting value to the final
          pixel type.


> 
>>>MINC files can have either per-slice or per-volume intensity
>>>transformation (i.e. slope and intercept); where should this be done?
>>>
>>
>>
>>        Potential intensity transformations must be done in the
>>        Read() method. See for example: itkGDCMImageIO.txx:
>>        lines 575-738.
>>
>>        Note that you must pay attention here to the pixel type,
>>        to make sure that the resulting intensities can be stored
>>        without losses.
> 
> 
> This is the sticking point.  MR images often have "real" values
> ranging up to 10^5 or 10^6 quantized down into 16 bits.  MINC may even
> store these using *signed* shorts.  If the Read() method has to stuff
> the signed shorts into the buffer provided AND rescale them back to
> the "real" values of 10^5 or 10^6, we're sure to lose precision.
> 


      Yes, you may end up lossing precision.
      But that happened when you quantized the values down to 16 bits.
      That is, the problem happened at *writing* time.
      There isn't much you can do at read time to correct for that
      lost information.  At read time... it is too late.

      If you really care about precision, (and given the low price
      of memory storage) you probably should consider storing the
      values as float and in this way, solve the real problem.



> DICOM must suffer from the same issue; how is it handled there?
> 
> ... steve peruses itkGDCMImageIO.cxx ...
> 
> So the rescaling blindly happens.  Does this not cause anyone
> problems?  Oh, wait ... here's a subtlety that I didn't pick up on
> previously: itkGDCMImageIO changes the IOComponentType to suit the
> data in the file.
> 
> So am I correct to summarize the semantics of the IOComponentType as
> follows?
> 
>   1. It is set according to the file data, NOT by the user of the
>      ImageIOBase (subclass)?



     The IOComponentType is intended to represent the Type used
     to store the pixel data in the file.



>   2. It represents the required component type to hold the data
>      *after* applying any intensity transformation to the pixel
>      values?
> 

      This is a matter of choice...

      you can certainly make your ImageIO class use an intermediate type.

      For example:

         a) File may contain pixel data encoded in 16 bits
            (but in need of rescaling)

         b) At read time, you can rescale the 16bits (short) data
            into float pixel type, and tell the ImageFileReader that
            this was the native pixel type.

         c) The user may have decided to instantiate the ImageFileReader
            using an image of pixel type "int".  Therefore the pixel
            data will be C-style casted from "float" (b) to "int".


       In other words, the ImageFileReader relies on the ImageIO to
       report what the "image file" pixel type is.  You can decide
       to report an intermediate pixel type *AS LONG AS* you make
       sure that for all purposes, this type is what the ImageFileReader
       will actually gets in the buffer filled up by the ImageIO.


> The intensity transformation in point 2 does include LUTs and DICOM's
> Palette Colour transformations (meaning that the output would be RGB)?
> What do you do when there are two LUTs for the image?
> 


        I'm not sure how you arrived to this conclusion.

        My reading of the code indicates that a palette colour
        LUT can (conditionally) be applied in lines 694-701,
        and that an intensity rescaling can (conditionally) be
        applied in lines 719-720.

        You seem to be assuming that it would be possible to
        receive DICOM images that have simultaneouly set a
        photometric interpretation to palette colour, *and*
        also set the Slope an Intercept to non-trivial values.

        We would need a DICOM guru to advice us here, on whether
        such combination can actually arise in practice.

        --

        From the seat of a software maintainer, this is a textbook
        example of why patching source code with "if()" is a road
        to manheim.

        The number of 'if' statements per number of lines of code,
        should be used as a metric of how much quick-and-dirty
        "patching" work has been made in a piece of source code.


        I'll go now to put my offerings in the State Machine shrine...




> Thanks,
> -Steve