[Insight-users] SparseFieldLevelSetImageFilter Boundary Conditions

[Brady McCary]

insight-users,

I think I have discovered the issue.

1. The current implementation of SparseFieldLevelSetImageFilter uses
an Image to keep track of the sparse-field-algorithm-status for each
pixel. In addition, there is a special status which indicates that a
pixel is a boundary pixel.
2. Neighborhood iterators have an optional optimization, in which you
can tell it not to check for being on the boundary.
3. The current algorithm checks if it ever hits the boundary status in
the status image and then turns on the neighborhood iterator's
checking for boundary conditions.
4. So, if the sparse field never hits a boundary status pixel, you
don't have to use the neighborhood iterator boundary conditions, a
convenient optimization.

That would be a nice optimization to have, but the implementation is
incorrect for (at least) the following reasons.

1. The boundary-status is set before the initialization of the sparse
field, and being on the boundary does not depend on time. However, if
the sparse field intersects the boundary, then the boundary-status is
lost when the sparse field status values are set. The sparse field
intersecting the boundary is a reasonable and natural thing.
2. The boundary-status pixels which are not overwritten during
initialization do not receive updates in the same fashion as the
non-boundary pixels. This occasionally results in non-updated level
set pixels on the boundary.

I am attaching a patch for the .h and .txx files, named patch.h and
patch.txx, which undo this optimization. The results are now what I
expect. Additionally, these patches reduce the memory allocation by
caching the difference function's global data pointer in a similar way
to ParallelSparseFieldLevelSetImageFilter.

As far as I am concerned these patches are trivial and in the public domain.

Brady

On Fri, Feb 6, 2009 at 2:09 PM, Brady McCary <brady.mccary+ITK at gmail.com> wrote:
> insight-users,
>
> To add to the discussion, when I disable all code that computes and
> checks the ``on-the-boundary'' status, the level set evolution is
> exactly what I expect it to be.
>
> Brady
>
> On Fri, Feb 6, 2009 at 2:01 PM, Brady McCary <brady.mccary+ITK at gmail.com> wrote:
>> insight-users,
>>
>> I have found evidence that makes me strongly believe that this is in
>> fact a bug in ITK.
>>
>> In order to understand this discussion, note that the Sparse Field
>> algorithm has to keep track of the STATUS of every pixel (whether it
>> is moving between layers, etc.) so that it can properly track and
>> construct the layers used in the algorithm. In particular, there is a
>> status value which signifies that a particular pixel is on the
>> boundary of the image. This status should (logically) be constant
>> throughout the entire iteration. However, if the initialization of
>> level set is such that the zero level set intersects the boundary of
>> the image, then the ``on-the-boundary'' status is lost.
>>
>> I.e., the status of a pixel being on a boundary is assumed to be
>> constant, but in fact this status depends on the initialization of the
>> level set, because if the zero level set intersects the boundary of
>> the image, the status values associated with the layers close to the
>> zero level set overwrite the ``on-the-boundary'' status value.
>>
>> I have been viewing the the source code and I am not clear why the
>> ``on-the-image-boundary'' status is significant, but I am looking. If
>> the author of this class could comment, it would be greatly
>> appreciated.
>>
>> I have posted two sequences of images to illustrate my findings. In
>> particular, the images named sta_NNNN.nii are the status of each pixel
>> at iteration NNNN. In the sequence which has an initialization where
>> image boundary pixels are INSIDE, the
>>
>> http://www.utdallas.edu/~bcm052000/example.tar.bz2
>>
>> Brady
>>
>> On Fri, Feb 6, 2009 at 1:00 PM, Koen Vermeer <koen at vermeer.tv> wrote:
>>> On Fri, 2009-02-06 at 10:53 -0600, Brady McCary wrote:
>>>> As an update to this thread, I have observed the following. To make
>>> [...]
>>>> dynamics in my problem, it is clear that this is not appropriate.
>>>
>>> Not very useful, but: I experienced the same issues. In my case, I'm
>>> looking at objects that are larger than the image, meaning that the
>>> level set at the boundary should move - but often it doesn't. I simply
>>> ignored this problem for now, as I'm just exploring possible solutions.
>>>
>>> So, basically, this is just a 'me too'.
>>>
>>> Best,
>>> Koen
>>>
>>>
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>>>
>>> Please keep messages on-topic and check the ITK FAQ at: http://www.itk.org/Wiki/ITK_FAQ
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: patch.h
Type: text/x-chdr
Size: 78 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090216/f23693c7/attachment.h>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: patch.txx
Type: application/octet-stream
Size: 1572 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090216/f23693c7/attachment.obj>