1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #ifndef itkDiscreteLevelSetImage_hxx
20 #define itkDiscreteLevelSetImage_hxx
21
22 #include "itkDiscreteLevelSetImage.h"
23
24 namespace itk
25 {
26
27 // ----------------------------------------------------------------------------
28 template< typename TOutput, unsigned int VDimension >
29 typename DiscreteLevelSetImage< TOutput, VDimension >::GradientType
EvaluateGradient(const InputType & inputIndex) const30 DiscreteLevelSetImage< TOutput, VDimension >::EvaluateGradient( const InputType& inputIndex ) const
31 {
32 InputType inputIndexA = inputIndex;
33 InputType inputIndexB = inputIndex;
34
35 GradientType dx;
36
37 for( unsigned int dim = 0; dim < Dimension; dim++ )
38 {
39 inputIndexA[dim] += 1;
40 inputIndexB[dim] -= 1;
41
42 if( !this->IsInsideDomain( inputIndexA ) )
43 {
44 inputIndexA[dim] = inputIndex[dim];
45 }
46
47 if( !this->IsInsideDomain( inputIndexB ) )
48 {
49 inputIndexB[dim] = inputIndex[dim];
50 }
51
52 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
53 const auto valueB = static_cast< OutputRealType >( this->Evaluate( inputIndexB ) );
54
55 // division by 0 only if image is a single pixel
56 const OutputRealType scale = this->m_NeighborhoodScales[dim] / (inputIndexA[dim] - inputIndexB[dim]);
57
58 dx[dim] = ( valueA - valueB ) * scale;
59
60 inputIndexA[dim] = inputIndexB[dim] = inputIndex[dim];
61
62 }
63
64 return dx;
65 }
66
67 // ----------------------------------------------------------------------------
68 template< typename TOutput, unsigned int VDimension >
69 typename DiscreteLevelSetImage< TOutput, VDimension >::GradientType
70 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateForwardGradient(const InputType & inputIndex) const71 ::EvaluateForwardGradient( const InputType& inputIndex ) const
72 {
73 const auto centerValue = static_cast< OutputRealType >( this->Evaluate( inputIndex ) );
74
75 InputType inputIndexA = inputIndex;
76
77 GradientType dx;
78
79 for( unsigned int dim = 0; dim < Dimension; dim++ )
80 {
81 inputIndexA[dim] += 1;
82
83 if( !this->IsInsideDomain( inputIndexA ) )
84 {
85 inputIndexA[dim] = inputIndex[dim];
86 }
87
88 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
89 const OutputRealType scale = this->m_NeighborhoodScales[dim];
90
91 dx[dim] = ( valueA - centerValue ) * scale;
92
93 inputIndexA[dim] = inputIndex[dim];
94 }
95
96 return dx;
97 }
98
99 // ----------------------------------------------------------------------------
100 template< typename TOutput, unsigned int VDimension >
101 typename DiscreteLevelSetImage< TOutput, VDimension >::GradientType
102 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateBackwardGradient(const InputType & inputIndex) const103 ::EvaluateBackwardGradient( const InputType& inputIndex ) const
104 {
105 const auto centerValue = static_cast< OutputRealType >( this->Evaluate( inputIndex ) );
106
107 InputType inputIndexA = inputIndex;
108
109 GradientType dx;
110
111 for( unsigned int dim = 0; dim < Dimension; dim++ )
112 {
113 inputIndexA[dim] -= 1;
114
115 if( !this->IsInsideDomain( inputIndexA ) )
116 {
117 inputIndexA[dim] = inputIndex[dim];
118 }
119
120 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
121 const OutputRealType scale = this->m_NeighborhoodScales[dim];
122
123 dx[dim] = ( centerValue - valueA ) * scale;
124
125 inputIndexA[dim] = inputIndex[dim];
126 }
127 return dx;
128 }
129
130 // ----------------------------------------------------------------------------
131 template< typename TOutput, unsigned int VDimension >
132 typename DiscreteLevelSetImage< TOutput, VDimension >::HessianType
133 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateHessian(const InputType & inputIndex) const134 ::EvaluateHessian( const InputType& inputIndex ) const
135 {
136 HessianType oHessian;
137
138 const auto centerValue = static_cast< OutputRealType >( this->Evaluate( inputIndex ) );
139
140 InputType inputIndexA = inputIndex;
141 InputType inputIndexB = inputIndex;
142
143 InputType inputIndexAa;
144 InputType inputIndexBa;
145 InputType inputIndexCa;
146 InputType inputIndexDa;
147
148 for( unsigned int dim1 = 0; dim1 < Dimension; dim1++ )
149 {
150 inputIndexA[dim1] += 1;
151 inputIndexB[dim1] -= 1;
152
153 if( !this->IsInsideDomain( inputIndexA ) )
154 {
155 inputIndexA[dim1] = inputIndex[dim1];
156 }
157
158 if( !this->IsInsideDomain( inputIndexB ) )
159 {
160 inputIndexB[dim1] = inputIndex[dim1];
161 }
162
163 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
164 const auto valueB = static_cast< OutputRealType >( this->Evaluate( inputIndexB ) );
165
166 oHessian[dim1][dim1] = ( valueA + valueB - 2.0 * centerValue )
167 * itk::Math::sqr( this->m_NeighborhoodScales[dim1] );
168
169 inputIndexAa = inputIndexB;
170 inputIndexBa = inputIndexB;
171
172 inputIndexCa = inputIndexA;
173 inputIndexDa = inputIndexA;
174
175 for( unsigned int dim2 = dim1 + 1; dim2 < Dimension; dim2++ )
176 {
177 inputIndexAa[dim2] -= 1;
178 inputIndexBa[dim2] += 1;
179
180 inputIndexCa[dim2] -= 1;
181 inputIndexDa[dim2] += 1;
182
183 if( !this->IsInsideDomain( inputIndexAa ) )
184 {
185 inputIndexAa[dim2] = inputIndexB[dim2];
186 }
187
188 if( !this->IsInsideDomain( inputIndexBa ) )
189 {
190 inputIndexBa[dim2] = inputIndexB[dim2];
191 }
192
193 if( !this->IsInsideDomain( inputIndexCa ) )
194 {
195 inputIndexCa[dim2] = inputIndexA[dim2];
196 }
197
198 if( !this->IsInsideDomain( inputIndexDa ) )
199 {
200 inputIndexDa[dim2] = inputIndexA[dim2];
201 }
202
203 const auto valueAa = static_cast< OutputRealType >( this->Evaluate( inputIndexAa ) );
204 const auto valueBa = static_cast< OutputRealType >( this->Evaluate( inputIndexBa ) );
205 const auto valueCa = static_cast< OutputRealType >( this->Evaluate( inputIndexCa ) );
206 const auto valueDa = static_cast< OutputRealType >( this->Evaluate( inputIndexDa ) );
207
208 oHessian[dim1][dim2] = oHessian[dim2][dim1] =
209 0.25 * ( valueAa - valueBa - valueCa + valueDa )
210 * this->m_NeighborhoodScales[dim1] * this->m_NeighborhoodScales[dim2];
211
212 inputIndexAa[dim2] = inputIndexB[dim2];
213 inputIndexBa[dim2] = inputIndexB[dim2];
214
215 inputIndexCa[dim2] = inputIndexA[dim2];
216 inputIndexDa[dim2] = inputIndexA[dim2];
217 }
218
219 inputIndexA[dim1] = inputIndex[dim1];
220 inputIndexB[dim1] = inputIndex[dim1];
221 }
222
223 return oHessian;
224 }
225
226 // ----------------------------------------------------------------------------
227 template< typename TOutput, unsigned int VDimension >
228 typename DiscreteLevelSetImage< TOutput, VDimension >::OutputRealType
229 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateLaplacian(const InputType & inputIndex) const230 ::EvaluateLaplacian( const InputType& inputIndex ) const
231 {
232 OutputRealType oLaplacian = NumericTraits< OutputRealType >::ZeroValue();
233
234 const auto centerValue = static_cast< OutputRealType >( this->Evaluate( inputIndex ) );
235
236 InputType inputIndexA = inputIndex;
237 InputType inputIndexB = inputIndex;
238
239 for( unsigned int dim1 = 0; dim1 < Dimension; dim1++ )
240 {
241 inputIndexA[dim1] += 1;
242 inputIndexB[dim1] -= 1;
243
244 if( !this->IsInsideDomain( inputIndexA ) )
245 {
246 inputIndexA[dim1] = inputIndex[dim1];
247 }
248
249 if( !this->IsInsideDomain( inputIndexB ) )
250 {
251 inputIndexB[dim1] = inputIndex[dim1];
252 }
253
254 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
255 const auto valueB = static_cast< OutputRealType >( this->Evaluate( inputIndexB ) );
256
257 oLaplacian += ( valueA + valueB - 2.0 * centerValue )
258 * itk::Math::sqr(this->m_NeighborhoodScales[dim1]);
259
260 inputIndexA[dim1] = inputIndex[dim1];
261 inputIndexB[dim1] = inputIndex[dim1];
262 }
263
264 return oLaplacian;
265 }
266
267 // ----------------------------------------------------------------------------
268 template< typename TOutput, unsigned int VDimension >
269 void
270 DiscreteLevelSetImage< TOutput, VDimension >
Evaluate(const InputType & inputIndex,LevelSetDataType & data) const271 ::Evaluate( const InputType& inputIndex, LevelSetDataType& data ) const
272 {
273 // If it has not already been computed before
274 if( data.Value.m_Computed )
275 {
276 return;
277 }
278
279 data.Value.m_Value = this->Evaluate( inputIndex );
280 data.Value.m_Computed = true;
281 }
282
283 // ----------------------------------------------------------------------------
284 template< typename TOutput, unsigned int VDimension >
285 void
286 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateGradient(const InputType & inputIndex,LevelSetDataType & data) const287 ::EvaluateGradient( const InputType& inputIndex, LevelSetDataType& data ) const
288 {
289 if( data.Gradient.m_Computed )
290 {
291 return;
292 }
293
294 // If it has not already been computed before
295
296 // compute the gradient
297
298 InputType inputIndexA = inputIndex;
299 InputType inputIndexB = inputIndex;
300
301 for( unsigned int dim = 0; dim < Dimension; dim++ )
302 {
303 inputIndexA[dim] += 1;
304 inputIndexB[dim] -= 1;
305
306 if( !this->IsInsideDomain( inputIndexA ) )
307 {
308 inputIndexA[dim] = inputIndex[dim];
309 }
310
311 if( !this->IsInsideDomain( inputIndexB ) )
312 {
313 inputIndexB[dim] = inputIndex[dim];
314 }
315
316 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
317 const auto valueB = static_cast< OutputRealType >( this->Evaluate( inputIndexB ) );
318 const OutputRealType scale = this->m_NeighborhoodScales[dim] / (inputIndexA[dim] - inputIndexB[dim]);
319
320 data.Gradient.m_Value[dim] = ( valueA - valueB ) * scale;
321
322 inputIndexA[dim] = inputIndexB[dim] = inputIndex[dim];
323 }
324
325 data.Gradient.m_Computed = true;
326 }
327
328 // ----------------------------------------------------------------------------
329 template< typename TOutput, unsigned int VDimension >
330 void
331 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateHessian(const InputType & inputIndex,LevelSetDataType & data) const332 ::EvaluateHessian( const InputType& inputIndex, LevelSetDataType& data ) const
333 {
334 if( data.Hessian.m_Computed )
335 {
336 return;
337 }
338
339 if( !data.Value.m_Computed )
340 {
341 data.Value.m_Computed = true;
342 data.Value.m_Value = this->Evaluate( inputIndex );
343 }
344
345 // compute the hessian
346 auto centerValue = static_cast< OutputRealType >( data.Value.m_Value );
347
348 InputType inputIndexA = inputIndex;
349 InputType inputIndexB = inputIndex;
350
351 InputType inputIndexAa;
352 InputType inputIndexBa;
353 InputType inputIndexCa;
354 InputType inputIndexDa;
355
356 bool backward = data.BackwardGradient.m_Computed;
357 bool forward = data.ForwardGradient.m_Computed;
358
359 for( unsigned int dim1 = 0; dim1 < Dimension; dim1++ )
360 {
361 inputIndexA[dim1] += 1;
362 inputIndexB[dim1] -= 1;
363
364 if( !this->IsInsideDomain( inputIndexA ) )
365 {
366 inputIndexA[dim1] = inputIndex[dim1];
367 }
368
369 if( !this->IsInsideDomain( inputIndexB ) )
370 {
371 inputIndexB[dim1] = inputIndex[dim1];
372 }
373
374 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
375 const auto valueB = static_cast< OutputRealType >( this->Evaluate( inputIndexB ) );
376
377 data.Hessian.m_Value[dim1][dim1] =
378 ( valueA + valueB - 2.0 * centerValue ) * itk::Math::sqr( this->m_NeighborhoodScales[dim1] );
379
380 if( !backward )
381 {
382 data.BackwardGradient.m_Computed = true;
383 data.BackwardGradient.m_Value[dim1] =
384 ( centerValue - valueB ) * this->m_NeighborhoodScales[dim1];
385 }
386
387 if( !forward )
388 {
389 data.ForwardGradient.m_Computed = true;
390 data.ForwardGradient.m_Value[dim1] =
391 ( valueA - centerValue ) * this->m_NeighborhoodScales[dim1];
392 }
393
394 inputIndexAa = inputIndexB;
395 inputIndexBa = inputIndexB;
396
397 inputIndexCa = inputIndexA;
398 inputIndexDa = inputIndexA;
399
400 for( unsigned int dim2 = dim1 + 1; dim2 < Dimension; dim2++ )
401 {
402 inputIndexAa[dim2] -= 1;
403 inputIndexBa[dim2] += 1;
404
405 inputIndexCa[dim2] -= 1;
406 inputIndexDa[dim2] += 1;
407
408 if( !this->IsInsideDomain( inputIndexAa ) )
409 {
410 inputIndexAa[dim2] = inputIndexB[dim2];
411 }
412
413 if( !this->IsInsideDomain( inputIndexBa ) )
414 {
415 inputIndexBa[dim2] = inputIndexB[dim2];
416 }
417
418 if( !this->IsInsideDomain( inputIndexCa ) )
419 {
420 inputIndexCa[dim2] = inputIndexA[dim2];
421 }
422
423 if( !this->IsInsideDomain( inputIndexDa ) )
424 {
425 inputIndexDa[dim2] = inputIndexA[dim2];
426 }
427
428 const auto valueAa = static_cast< OutputRealType >( this->Evaluate( inputIndexAa ) );
429 const auto valueBa = static_cast< OutputRealType >( this->Evaluate( inputIndexBa ) );
430 const auto valueCa = static_cast< OutputRealType >( this->Evaluate( inputIndexCa ) );
431 const auto valueDa = static_cast< OutputRealType >( this->Evaluate( inputIndexDa ) );
432
433 data.Hessian.m_Value[dim1][dim2] =
434 data.Hessian.m_Value[dim2][dim1] =
435 0.25 * ( valueAa - valueBa - valueCa + valueDa )
436 * this->m_NeighborhoodScales[dim1] * this->m_NeighborhoodScales[dim2];
437
438 inputIndexAa[dim2] = inputIndexB[dim2];
439 inputIndexBa[dim2] = inputIndexB[dim2];
440
441 inputIndexCa[dim2] = inputIndexA[dim2];
442 inputIndexDa[dim2] = inputIndexA[dim2];
443 }
444
445 inputIndexA[dim1] = inputIndex[dim1];
446 inputIndexB[dim1] = inputIndex[dim1];
447 }
448
449 data.Hessian.m_Computed = true;
450 }
451
452 // ----------------------------------------------------------------------------
453 template< typename TOutput, unsigned int VDimension >
454 typename DiscreteLevelSetImage< TOutput, VDimension >::OutputRealType
455 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateMeanCurvature(const InputType & inputIndex) const456 ::EvaluateMeanCurvature( const InputType& inputIndex ) const
457 {
458 OutputRealType oValue = NumericTraits< OutputRealType >::ZeroValue();
459
460 HessianType hessian = this->EvaluateHessian( inputIndex );
461 GradientType grad = this->EvaluateGradient( inputIndex );
462
463 for( unsigned int i = 0; i < Dimension; i++ )
464 {
465 for( unsigned int j = 0; j < Dimension; j++ )
466 {
467 if( j != i )
468 {
469 oValue -= grad[i] * grad[j] * hessian[i][j];
470 oValue += hessian[j][j] * grad[i] * grad[i];
471 }
472 }
473 }
474
475 OutputRealType gradNorm = grad.GetNorm();
476
477 if( gradNorm > itk::Math::eps )
478 {
479 oValue /= ( gradNorm * gradNorm * gradNorm );
480 }
481 else
482 {
483 oValue /= ( NumericTraits< OutputRealType >::OneValue() + gradNorm );
484 }
485
486 return oValue;
487 }
488
489 // ----------------------------------------------------------------------------
490 template< typename TOutput, unsigned int VDimension >
491 void
492 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateLaplacian(const InputType & inputIndex,LevelSetDataType & data) const493 ::EvaluateLaplacian( const InputType& inputIndex, LevelSetDataType& data ) const
494 {
495 if( data.Laplacian.m_Computed )
496 {
497 return;
498 }
499
500 if( !data.Value.m_Computed )
501 {
502 data.Value.m_Value = this->Evaluate( inputIndex );
503 data.Value.m_Computed = true;
504 }
505
506 const auto centerValue = static_cast< OutputRealType >( data.Value.m_Value );
507
508 InputType inputIndexA =inputIndex;
509 InputType inputIndexB = inputIndex;
510
511 for( unsigned int dim1 = 0; dim1 < Dimension; dim1++ )
512 {
513 inputIndexA[dim1] += 1;
514 inputIndexB[dim1] -= 1;
515
516 if( !this->IsInsideDomain( inputIndexA ) )
517 {
518 inputIndexA[dim1] = inputIndex[dim1];
519 }
520
521 if( !this->IsInsideDomain( inputIndexB ) )
522 {
523 inputIndexB[dim1] = inputIndex[dim1];
524 }
525
526 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
527 const auto valueB = static_cast< OutputRealType >( this->Evaluate( inputIndexB ) );
528
529 data.Laplacian.m_Value +=
530 ( valueA + valueB - 2.0 * centerValue ) * itk::Math::sqr( this->m_NeighborhoodScales[dim1] );
531
532 inputIndexA[dim1] = inputIndex[dim1];
533 inputIndexB[dim1] = inputIndex[dim1];
534 }
535
536 data.Laplacian.m_Computed = true;
537 }
538
539 // ----------------------------------------------------------------------------
540 template< typename TOutput, unsigned int VDimension >
541 void
542 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateMeanCurvature(const InputType & inputIndex,LevelSetDataType & data) const543 ::EvaluateMeanCurvature( const InputType& inputIndex, LevelSetDataType& data ) const
544 {
545 if( !data.MeanCurvature.m_Computed )
546 {
547 if( !data.Hessian.m_Computed )
548 {
549 this->EvaluateHessian( inputIndex, data );
550 }
551
552 if( !data.Gradient.m_Computed )
553 {
554 this->EvaluateGradient( inputIndex, data );
555 }
556
557 if( !data.GradientNorm.m_Computed )
558 {
559 this->EvaluateGradientNorm( inputIndex, data );
560 }
561
562 data.MeanCurvature.m_Computed = true;
563 data.MeanCurvature.m_Value = NumericTraits< OutputRealType >::ZeroValue();
564
565 for( unsigned int i = 0; i < Dimension; i++ )
566 {
567 for( unsigned int j = 0; j < Dimension; j++ )
568 {
569 if( j != i )
570 {
571 data.MeanCurvature.m_Value -= data.Gradient.m_Value[i]
572 * data.Gradient.m_Value[j] * data.Hessian.m_Value[i][j];
573 data.MeanCurvature.m_Value += data.Hessian.m_Value[j][j]
574 * data.Gradient.m_Value[i] * data.Gradient.m_Value[i];
575 }
576 }
577 }
578
579 OutputRealType temp = data.GradientNorm.m_Value;
580
581 if( temp > itk::Math::eps )
582 {
583 data.MeanCurvature.m_Value /= ( temp * temp * temp );
584 }
585 else
586 {
587 data.MeanCurvature.m_Value /= ( NumericTraits< OutputRealType >::OneValue() + temp );
588 }
589 }
590 }
591
592 // ----------------------------------------------------------------------------
593 template< typename TOutput, unsigned int VDimension >
594 void
595 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateForwardGradient(const InputType & inputIndex,LevelSetDataType & data) const596 ::EvaluateForwardGradient( const InputType& inputIndex, LevelSetDataType& data ) const
597 {
598 if( data.ForwardGradient.m_Computed )
599 {
600 return;
601 }
602
603 // compute the gradient
604 if( !data.Value.m_Computed )
605 {
606 data.Value.m_Computed = true;
607 data.Value.m_Value = this->Evaluate( inputIndex );
608 }
609
610 const auto centerValue = static_cast< OutputRealType >( data.Value.m_Value );
611
612 InputType inputIndexA = inputIndex;
613
614 GradientType dx;
615
616 for( unsigned int dim = 0; dim < Dimension; dim++ )
617 {
618 inputIndexA[dim] += 1;
619
620 if( !this->IsInsideDomain( inputIndexA ) )
621 {
622 inputIndexA[dim] = inputIndex[dim];
623 }
624
625 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
626 const OutputRealType scale = this->m_NeighborhoodScales[dim];
627
628 dx[dim] = ( valueA - centerValue ) * scale;
629
630 inputIndexA[dim] = inputIndex[dim];
631 }
632 data.ForwardGradient.m_Value = dx;
633
634 data.ForwardGradient.m_Computed = true;
635 }
636
637 // ----------------------------------------------------------------------------
638 template< typename TOutput, unsigned int VDimension >
639 void
640 DiscreteLevelSetImage< TOutput, VDimension >
EvaluateBackwardGradient(const InputType & inputIndex,LevelSetDataType & data) const641 ::EvaluateBackwardGradient( const InputType& inputIndex, LevelSetDataType& data ) const
642 {
643 if( data.BackwardGradient.m_Computed )
644 {
645 return;
646 }
647
648 // compute the gradient
649 if( !data.Value.m_Computed )
650 {
651 data.Value.m_Value = this->Evaluate( inputIndex );
652 data.Value.m_Computed = true;
653 }
654
655 const auto centerValue = static_cast< OutputRealType >( data.Value.m_Value );
656
657 InputType inputIndexA = inputIndex;
658
659 GradientType dx;
660
661 for( unsigned int dim = 0; dim < Dimension; dim++ )
662 {
663 inputIndexA[dim] -= 1;
664
665 if( !this->IsInsideDomain( inputIndexA ) )
666 {
667 inputIndexA[dim] = inputIndex[dim];
668 }
669
670 const auto valueA = static_cast< OutputRealType >( this->Evaluate( inputIndexA ) );
671 const OutputRealType scale = this->m_NeighborhoodScales[dim];
672
673 dx[dim] = ( centerValue - valueA ) * scale;
674
675 inputIndexA[dim] = inputIndex[dim];
676 }
677 data.BackwardGradient.m_Value = dx;
678
679 data.BackwardGradient.m_Computed = true;
680 }
681
682 // ----------------------------------------------------------------------------
683 template< typename TOutput, unsigned int VDimension >
684 void
685 DiscreteLevelSetImage< TOutput, VDimension >
Initialize()686 ::Initialize()
687 {
688 Superclass::Initialize();
689 }
690
691 // ----------------------------------------------------------------------------
692 template< typename TOutput, unsigned int VDimension >
693 void
694 DiscreteLevelSetImage< TOutput, VDimension >
CopyInformation(const DataObject * data)695 ::CopyInformation(const DataObject *data)
696 {
697 Superclass::CopyInformation( data );
698
699 const auto * LevelSet = dynamic_cast< const Self * >( data );
700 if ( !LevelSet )
701 {
702 // pointer could not be cast back down
703 itkExceptionMacro( << "itk::DiscreteLevelSetImage::CopyInformation() cannot cast "
704 << typeid( data ).name() << " to "
705 << typeid( Self * ).name() );
706 }
707 }
708
709 // ----------------------------------------------------------------------------
710 template< typename TOutput, unsigned int VDimension >
711 void
712 DiscreteLevelSetImage< TOutput, VDimension >
Graft(const DataObject * data)713 ::Graft( const DataObject* data )
714 {
715 Superclass::Graft( data );
716 const auto * LevelSet = dynamic_cast< const Self* >( data );
717
718 if ( !LevelSet )
719 {
720 // pointer could not be cast back down
721 itkExceptionMacro( << "itk::DiscreteLevelSetImage::CopyInformation() cannot cast "
722 << typeid( data ).name() << " to "
723 << typeid( Self * ).name() );
724 }
725
726 this->m_NeighborhoodScales = LevelSet->m_NeighborhoodScales;
727 }
728
729 }
730 #endif // itkDiscreteLevelSetImage_hxx
731