[Insight-developers] Get/SetPixel() with taking an offsetparameter
Luis Ibanez
luis.ibanez at kitware.com
Thu Feb 19 14:37:43 EST 2009
Hi Gaetan,
I see your point, and agree with Jim's comments.
It seems that adding the method
itk::Image::SetPixel(offset, val);
is reasonable.
After all, as you pointed out, developers can
already do this via the GetBufferPointer(),
therefore, by adding the SetPixel(offset,value)
method we are not reducing the encapsulation of
the image anything more than what already is.
My suggestion will be to hold on this change
until the release of ITK 3.12. (just as you
suggested already).
This will be an addition to a very fundamental
class, and if we get it wrong just before a
release, we will be stuck with that API forever.
We should be cutting the release by the end of
next week (Feb 26~28).
Luis
------------------------
Gaëtan Lehmann wrote:
>
> Le 19 févr. 09 à 16:21, Luis Ibanez a écrit :
>
>>
>> Hi Gaetan,
>>
>> It seems that what you actually need is a new customized Iterator,
>> not a new pixel access method in the itkImage itself.
>>
>> You could take the code that selects the pixel to be visited,
>> and you could move it from the Filter into a new specialized
>> Iterator, (that may be reused by others).
>>
>> Inside the new Iterator you have direct access to the
>> internals of the itk::Image.
>>
>
> Hi Luis,
>
> I guess that's an objection :-)
>
> I must agree that iterators has many good points, and they can make
> developers life way easier and the code much efficient.
>
> However, they are sometime just not what we need. In the STL for
> example, all the containers are accessible through iterators. This
> gives us an easy way to visit their content. But the content is also
> available in other ways, and that's often what gives them there
> characteristics. Would you find nice to be forced to use an iterator to
> access to Nth element of a vector?
>
> For several algorithms, we need to access to the image data in a
> constant time and in non sequencial order, as we would do with a
> std::vector or std::deque. It's not a matter of personal choice, that's
> just that the algorithms are made that way, and that's often what is
> making them efficient.
>
> Developing a new iterator to be able to access the data the way the
> algorithm needs it is just overkill, and has some bad side effects,
> like adding an extra level of complexity, which may convince the
> developer to find some workarounds, like copying the input data in a
> data structure suitable for a simple linear access like std::vector, or
> use the itk::Indexes directly, to store them or to use them with
> Set/GetPixel(), which is both memory and computationally less efficient.
> That's not a supposition - Richard, an experimented developer, has done
> that in his code.
>
> I think that the iterators are simply not the best way to go for
> everything. Many published algorithm are using non sequential access to
> the image, and, at this time, can't be implemented in ITK as
> efficiently as they should because of the lack of linear access the
> that image.
> Also, as Jim said, exposing this code has no effect on the other memory
> organizations we could think about.
>
> Could you see any reason to not go that way?
>
> Gaëtan
>
>
>>
>> Luis
>>
>>
>> ----------------------
>> Gaëtan Lehmann wrote:
>>
>>> Hi,
>>> I recalled I've already asked for that, but I'm a bit surprised to
>>> see that it was more than two years ago.
>>> I'm working on improving a bit the memory usage in Richard Beare's
>>> area opening/closing filters, and what would help me a lot is to
>>> have the Get/SetPixel() method which are taking that offset parameter.
>>> Any objection to the addition of those methods (after the release)?
>>> Regards,
>>> Gaëtan
>>> Le 20 déc. 06 à 22:33, Miller, James V ((GE, Research)) a écrit :
>>>
>>>> Gaetan,
>>>>
>>>> You might also want to look at using std::deque instead of
>>>> std::vector. I probably overuse std::vector and should use
>>>> std::deque more. Deque doesn't use contiguous memory, so when it
>>>> needs to increase its size it just allocates another block (it is
>>>> basically a link list of fixed length vectors).
>>>>
>>>> Using a deque, you won't have to prescan so you can allocate
>>>> vectors of the exact length. You can be lazy and use push_back()
>>>> and have simpler code. Plus you don't pay the log2-type allocation
>>>> penalty of a vector as the deque grows. So you should get a
>>>> performance increase over not preallocating the vectors.
>>>>
>>>> So you might be able to use an std::map<PixelType,
>>>> std::deque<IndexType> > or an std::map<PixelType,
>>>> std::deque<OffsetValueType> >.
>>>>
>>>> If we do introduce an image type with non-continguous memory, I
>>>> will push that we maintain Image with an API that provides a one-
>>>> dimensional view of ComputeIndex()/ComputeOffset(). What I would
>>>> want to "remove" is the access to GetBufferPointer() and hence
>>>> prevent people writing code equivalent to
>>>>
>>>> x = img->GetBufferPointer() + img->ComputeOffset(index);
>>>> *x = val;
>>>>
>>>> But still allow a 1D access
>>>>
>>>> offset = img->ComputeOffset(index);
>>>> img->SetPixel(offset, val);
>>>>
>>>> Or
>>>>
>>>> for (i=0; i < numPixel; ++i)
>>>> img->SetPixel(I, val);
>>>>
>>>> So from a programming perspective, you can think of Image as
>>>> having a continguous block of memory. But internally, we'd have
>>>> the freedom to divide the image up differently.
>>>>
>>>> We can do this if we modify the iterators and the IO. We'd have
>>>> to pay a double dereference penalty to access a pixel. And we'd
>>>> upset anyone out in the world that is relying on
>>>> GetBufferPointer(). But we may be able to do this using new image
>>>> types and traits for creating iterators, etc....
>>>>
>>>> Those are my thoughts.
>>>>
>>>> Jim
>>>>
>>>>
>>>> -----Original Message-----
>>>> From: insight-developers-bounces+millerjv=crd.ge.com at itk.org
>>>> [mailto:insight-developers-bounces+millerjv=crd.ge.com at itk.org ] On
>>>> Behalf Of Gaetan Lehmann
>>>> Sent: Wednesday, December 20, 2006 4:37 AM
>>>> To: Miller, James V (GE, Research); Luis Ibanez
>>>> Cc: insight-developers at itk.org
>>>> Subject: Re: [Insight-developers] Get/SetPixel() with taking an
>>>> offsetparameter
>>>>
>>>>
>>>> Here is what I have tried about std::map< PixelType, std::vector<
>>>> IndexType > >:
>>>>
>>>> (a) - make a first pass to compute the number of pixels for all
>>>> the values
>>>> - create the map, and reserve the exact size which will be
>>>> required for all the vector
>>>> - fill the vectors
>>>>
>>>> (b) replace std::map< PixelType, std::vector< IndexType > > by
>>>> std::map< PixelType, std::vector< OffsetValueType > >, the
>>>> conversion being done by the ComputeOffset() and computeIndex()
>>>> methods.
>>>>
>>>> both (a) and (b) where done to decrease the memory usage - and it
>>>> works very well regarding that.
>>>> More surprising, both methods also give better performance. With
>>>> both, my test run in 8.2 sec instead of 10.7 before the changes.
>>>> I think that's because the vectors no more need to increase their
>>>> size, and that it avoid the cost of this increase, and because
>>>> manipulating 3 times less memory compensate the cost of the
>>>> conversion.
>>>>
>>>> The only problem is the much compex code, with additional
>>>> ComputeOffset() and ComputeIndex() methods in it because of (b),
>>>> and the 2 more steps done in (a).
>>>> And there is still the cost of the conversion index <-> offset
>>>> that I would like to avoid
>>>>
>>>> Gaetan
>>>>
>>>>
>>>> On Tue, 19 Dec 2006 22:48:52 +0100, Miller, James V (GE, Research)
>>>> <millerjv at crd.ge.com > wrote:
>>>>
>>>>> Wrt to the std::map< PixelType, std::vector< IndexType > >, I have
>>>>> thought about this a few times and always verred away because of the
>>>>> memory requirements.
>>>>>
>>>>> If we are adopting these high memory implementations, we may want to
>>>>> have an option to switch to a lower memory albeit slower version.
>>>>>
>>>>> I just had a thought that maybe some of these algorithms that use a
>>>>> map as above could be done in several iterations. Less than the
>>>>> number of pixel intensity values but more than one. Say you divide
>>>>> the intensity space into k bands, say you build the map with the top
>>>>> 10% of the pixels, then do whatever processing you can do, then build
>>>>> the map with the next 10%, etc.
>>>>>
>>>>>
>>>>> -----Original Message-----
>>>>> From: insight-developers-bounces+millerjv=crd.ge.com at itk.org
>>>>> [mailto:insight-developers-bounces+millerjv=crd.ge.com at itk.org] On
>>>>> Behalf Of Gaetan Lehmann
>>>>> Sent: Tuesday, December 12, 2006 6:13 AM
>>>>> To: Luis Ibanez
>>>>> Cc: insight-developers at itk.org
>>>>> Subject: Re: [Insight-developers] Get/SetPixel() with taking an
>>>>> offsetparameter
>>>>>
>>>>>
>>>>> 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=s
>>>>>>> ummary),
>>>>>>> ...) 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
>>>> _______________________________________________
>>>> Insight-developers mailing list
>>>> Insight-developers at itk.org
>>>> http://www.itk.org/mailman/listinfo/insight-developers
>
>
More information about the Insight-developers
mailing list