[Insight-developers] Get/SetPixel() with taking an offset parameter

[Gaetan Lehmann]

Hi,

>
> Hi Gaetan,
>
> That's an interesting option.
>
> However, please note that a better way of approaching this problem
> is to create an Iterator that is specific to the way you want to
> walk in your image.
>
> It seems that in your algorithm you are already defining a clever
> way of finding the offsets that will lead you to the next pixel
> of interest.  What we probably want to explore is the possibility
> of creating a new Image Iterator that will do the equivalent job.
>

What is done in ReconstructionByDilationImageFilter currently work only  
with shape iterators, so having the same operators on the other iterators  
would be great. I'm not sure that a dedicated iterator is needed.

> In general you want to avoid having direct access to the image buffer.
> One of the reasons for it is that soon we will be introducing in ITK
> an alternative memory layout where the pixel data is not stored in a
> contiguous block of memory. Instead, the pixel data may be stored as
> a collection of disconnected slices.
>
> Having a GetPixel()/SetPixel() methods with direct access to the
> buffer data will not work in this memory layout. We already will have
> to manage the modification of the existing ComputOffset() method.
>

I'm not sure that modifying the Get/SetPixel() method is the best way to  
go, but what is proposed has no extra cost, because the ComputeIndex() and  
ComputeOffset() are already exposed, and there is nothing more done that  
using these methods.

> You will also have a harder time making sure that the algorithm works
> in N-Dimensional images. It is in general, easier, to delegate to
> Iterators the complex details of dealing with the N-Dimensionality
> of the image.

It will not be harder, because an iterator is used to get the offsets :-)
Offsets are only a cheaper representation of Index, and the conversion can  
be easily done with the ComputeIndex() and ComputeOffset() methods.

>
> What is the specific way in which you need to walk through the image
> in the morphological reconstruction filters that you are preparing ?
>

A common step - at least common enough to make me implement it several  
times - is to sort the pixels by there value to process the highest value  
first for example. The data structure used to do that is:

   std::map< PixelType, std::vector< IndexType > >

after have filled this data structure with all the pixels in the image,  
it's easy to process all the pixels in the image in the right order.
If we take care of the memory allocation of the vectors (which require an  
extra iteration step), for a big image (1024 x 1024 x 80), this structure  
take about 1024*1024*80*4*3/1024^2 = 960 MB on a 32 bits system - twice  
more a 64bits one. When using offsets, the cost is reduced by the image  
dimension, 3 in that case = 320 MB. The difference is quite significant :-)

This example is "only" an optimization problem, as the pixels can be  
processed in the right order without that, by iterating over the image as  
many time as we have pixel values in the image (and one more time to find  
the maximum).

The next case is more difficult.
The component tree representation (also called min tree or max tree) don't  
store the image in a array, but in a tree of connected components in the  
image. All the nodes of that tree are associated to a list of pixels in  
the image - a list of itk::Index. That's no more an algorithm  
optimization, but a data representation problem, and the itk way of  
designing a pixel location make it highly inefficient compared to what can  
be done in other image analysis libs.

It is already possible to use offsets instead of itk::Index. That's what  
I'm doing in the component trees, for memory efficiency (please, don't  
make the ComputeOffset() and ComputeIndex() method deprecated !). However,  
the permanent conversion itk::Index <-> offset seem to be an important  
extra cost which should be avoided.
I have not measured the extra cost of the conversion, so perhaps I'm  
wrong, and the conversion cost is not significant compared to the rest of  
the algorithms.
Given the comments about performance in the code, I think that the  
measures of the execution time of those conversion have already been done.
Can you give a pointer to those results ?

Thanks,

Gaetan


>
>    Please let us know,
>
>
>       Thanks
>
>
>          Luis
>
>
>
> ----------------------
> Gaetan Lehmann wrote:
>>  Hi,
>>  Currently, the GetPixel() and SetPixel() methods of the itk::Image  
>> class  are taking a itk::Index parameter to know which pixel is  
>> concerned in the  image.
>> The itk::Image also have the ComputeIndex() and ComputeOffset() methods  
>> to  manipulate directly the offset in the buffer, and to convert this  
>> offset  from/to an index.
>>  Manipulating indexes in an algorithm (like in the current  
>> morphological  reconstruction filters, in the component tree  
>> contribution I'm preparing   
>> (http://voxel.jouy.inra.fr/darcsweb/darcsweb.cgi?r=componentTree;a=summary),   
>> ...) would be a lot more efficient, both in memory and in execution  
>> time,  by using the offsets instead of the indexes: store an itk::Index  
>> cost 3  long (with dim=3), while an offset is only 1, and there is no  
>> need to  convert the offset to an index and the index to an offset.
>> Unfortunately, there is no Get/SetPixel() method able to work with an   
>> offset.
>>  I think that it would be a good idea to implement them as follow (not   
>> tested yet, that's an example):
>>    in itkImage.h:
>>  void SetPixel(const OffsetValueType &offset, const TPixel& value)
>>   {
>>   (*m_Buffer)[offset] = value;
>>   }
>>  const TPixel& GetPixel(const OffsetValueType &offset) const
>>   {
>>   return ( (*m_Buffer)[offset] );
>>   }
>>  TPixel& GetPixel(const OffsetValueType &offset)
>>   {
>>   return ( (*m_Buffer)[offset] );
>>   }
>>  TPixel & operator[](const OffsetValueType &offset)
>>   { return this->GetPixel(offset); }
>>  const TPixel& operator[](const OffsetValueType &offset) const
>>   { return this->GetPixel(offset); }
>>   void SetPixel(const IndexType &index, const TPixel& value)
>>   {
>>   typename Superclass::OffsetValueType offset =  
>> this->ComputeOffset(index);
>>   this->SetPixel( offset, value );
>>   }
>>  const TPixel& GetPixel(const IndexType &index) const
>>   {
>>   typename Superclass::OffsetValueType offset =  
>> this->ComputeOffset(index);
>>   return this->GetPixel( offset );
>>   }
>>  TPixel& GetPixel(const IndexType &index)
>>   {
>>   typename Superclass::OffsetValueType offset =  
>> this->ComputeOffset(index);
>>   return this->GetPixel( offset );
>>   }
>>  TPixel & operator[](const IndexType &index)
>>   { return this->GetPixel(index); }
>>  const TPixel& operator[](const IndexType &index) const
>>   { return this->GetPixel(index); }
>>    in itkImageAdaptor.h:
>>  void SetPixel(const IndexType &index, const PixelType & value)
>>   { m_PixelAccessor.Set( m_Image->GetPixel(index), value ); }
>>  PixelType GetPixel(const IndexType &index) const
>>   { return m_PixelAccessor.Get( m_Image->GetPixel(index) ); }
>>  PixelType operator[](const IndexType &index) const
>>   { return m_PixelAccessor.Get( m_Image->GetPixel(index) ); }
>>   void SetPixel(const OffsetValueType &offset, const PixelType & value)
>>   { m_PixelAccessor.Set( m_Image->GetPixel(offset), value ); }
>>  PixelType GetPixel(const OffsetValueType &offset) const
>>   { return m_PixelAccessor.Get( m_Image->GetPixel(offset) ); }
>>  PixelType operator[](const OffsetValueType &offset) const
>>   { return m_PixelAccessor.Get( m_Image->GetPixel(offset) ); }
>>    Does it seem reasonable to do that ?
>>  Thanks,
>>  Gaetan
>>



-- 
Gaëtan Lehmann
Biologie du Développement et de la Reproduction
INRA de Jouy-en-Josas (France)
tel: +33 1 34 65 29 66    fax: 01 34 65 29 09
http://voxel.jouy.inra.fr