[Insight-developers] question about ComputeGlobalTimeStep in LevelSetFunction

[Arnaud Gelas]

Hi guys,

I was looking at this method (LevelSetFunction:: ComputeGlobalTimeStep  
I copy the code below), and I have several questions...

I was wondering what's the theoretical and numerical difference  
between m_DT and m_WaveDT?

what happens when d->m_MaxCurvatureChange is negative?
I thought the global time step had to be positive?
Is d->m_MaxCurvatureChange always positive? If so why using a  
vnl_math_abs( d->m_MaxCurvatureChange ) in the first condition?
Am I missing something?

Why dt = m_DT / max( ( m_MaxAdvectionChange + m_PropagationChange ),  
m_MaxCurvatureChange )?
and not dt = m_DT / ( abs( m_MaxAdvectionChange ) +  
abs( m_PropagationChange ) + abs( m_MaxCurvatureChange ) )?

Thanks in advance for your answers

Arnaud

224 template< class TImageType >
225 typename LevelSetFunction< TImageType >::TimeStepType
226 LevelSetFunction<TImageType>
227 ::ComputeGlobalTimeStep(void *GlobalData) const
228 {
229   TimeStepType dt;
230
231   GlobalDataStruct *d = (GlobalDataStruct *)GlobalData;
232
233   d->m_MaxAdvectionChange += d->m_MaxPropagationChange;
234
235   if (vnl_math_abs(d->m_MaxCurvatureChange) > 0.0)
236     {
237     if (d->m_MaxAdvectionChange > 0.0)
238       {
239       dt = vnl_math_min((m_WaveDT / d->m_MaxAdvectionChange),
240                         (    m_DT / d->m_MaxCurvatureChange ));
241       }
242     else
243       {
244       dt = m_DT / d->m_MaxCurvatureChange;
245       }
246     }
247   else
248     {
249     if (d->m_MaxAdvectionChange > 0.0)
250       {
251       dt = m_WaveDT / d->m_MaxAdvectionChange;
252       }
253     else
254       {
255       dt = 0.0;
256       }
257     }
258
259   double maxScaleCoefficient = 0.0;
260   for (unsigned int i=0; i<ImageDimension; i++)
261     {
262     maxScaleCoefficient = vnl_math_max(this- 
 >m_ScaleCoefficients[i],maxScaleCoefficient);
263     }
264   dt /= maxScaleCoefficient;
265
266   // reset the values
267   d->m_MaxAdvectionChange   = NumericTraits<ScalarValueType>::Zero;
268   d->m_MaxPropagationChange = NumericTraits<ScalarValueType>::Zero;
269   d->m_MaxCurvatureChange   = NumericTraits<ScalarValueType>::Zero;
270
271   return dt;
272 }