[Insight-users] Comparison of Sethians equation and itkLevelSetFunction implementation

[Karl Krissian]

Hi,

I think the original code is correct,
vnl_math_sqr is the square and not the square root,

then there is the following operation which takes the square root:

propagation_term *= vcl_sqrt(propagation_gradient)

Regards,

Karl


Christoph Palm wrote:

>Dear all,
>
>to understand the levelSet implementation within ITK
>I went very deep into the code and find an inconsistancy
>between the cited bood of Sethian and the implementation
>of itkLevelSetFunction:
>
>The following lines code refer obviously to equation (6.15)
>on page 65 chapter 6 of Sethians famous bool about LevelSets:
>      
>for(i = 0; i< ImageDimension; i++)
>{
>   propagation_gradient += vnl_math_sqr(    
>                           vnl_math_max(gd->m_dx_backward[i], ZERO) )
>          + vnl_math_sqr( vnl_math_min(gd->m_dx_forward[i],  ZERO) );
>}
>
>To my mind, there is the square of each derivative missing and the
>square root should not be used only for one term.
>Therefore, the code should be like the following:
>
>for(i = 0; i< ImageDimension; i++)
>{
>   propagation_gradient +=    vnl_math_max(gd->m_dx_backward[i], ZERO) 
>                            * vnl_math_max(gd->m_dx_backward[i], ZERO)
>                            + vnl_math_min(gd->m_dx_forward[i],  ZERO) 
>                            * vnl_math_min(gd->m_dx_forward[i],  ZERO);
>}
>propagation_gradient = vnl_math_sqr(propagation_gradient);
>
>(the else part in the code should be changed accordingly) 
>
>Am I right or should I have to look closer into the code?
>Comments on that are highly appreciated.
>
>Thanks
>
>Christoph
>
>  
>