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 itkLevelSetEquationChanAndVeseExternalTerm_hxx
20 #define itkLevelSetEquationChanAndVeseExternalTerm_hxx
21
22 #include "itkLevelSetEquationChanAndVeseExternalTerm.h"
23
24 namespace itk
25 {
26 template< typename TInput, typename TLevelSetContainer >
27 LevelSetEquationChanAndVeseExternalTerm< TInput, TLevelSetContainer >
LevelSetEquationChanAndVeseExternalTerm()28 ::LevelSetEquationChanAndVeseExternalTerm()
29 {
30 this->m_TermName = "External Chan And Vese term";
31 this->m_RequiredData.insert( "Value" );
32 this->m_DomainMapImageFilter = nullptr;
33 this->m_CacheImage = nullptr;
34 }
35
36 template< typename TInput, typename TLevelSetContainer >
37 void LevelSetEquationChanAndVeseExternalTerm< TInput, TLevelSetContainer >
ComputeProduct(const LevelSetInputIndexType & iP,LevelSetOutputRealType & prod)38 ::ComputeProduct( const LevelSetInputIndexType& iP, LevelSetOutputRealType& prod )
39 {
40 this->ComputeProductTerm( iP, prod );
41 LevelSetPointer levelSet = this->m_LevelSetContainer->GetLevelSet( this->m_CurrentLevelSetId );
42 LevelSetOutputRealType value = levelSet->Evaluate( iP );
43 prod *= -(1 - this->m_Heaviside->Evaluate( -value ) );
44 }
45
46 template< typename TInput, typename TLevelSetContainer >
47 void LevelSetEquationChanAndVeseExternalTerm< TInput, TLevelSetContainer >
ComputeProductTerm(const LevelSetInputIndexType & iP,LevelSetOutputRealType & prod)48 ::ComputeProductTerm( const LevelSetInputIndexType& iP, LevelSetOutputRealType& prod )
49 {
50 prod = -1 * NumericTraits< LevelSetOutputRealType >::OneValue();
51
52 if( this->m_LevelSetContainer->HasDomainMap() )
53 {
54 if(this->m_DomainMapImageFilter == nullptr)
55 {
56 this->m_DomainMapImageFilter = this->m_LevelSetContainer->GetModifiableDomainMapFilter();
57 this->m_CacheImage = this->m_DomainMapImageFilter->GetOutput();
58 }
59 const LevelSetIdentifierType id = this->m_CacheImage->GetPixel( iP );
60
61 using DomainMapType = typename DomainMapImageFilterType::DomainMapType;
62 const DomainMapType domainMap = this->m_DomainMapImageFilter->GetDomainMap();
63 auto levelSetMapItr = domainMap.find(id);
64
65 if( levelSetMapItr != domainMap.end() )
66 {
67 const IdListType * idList = levelSetMapItr->second.GetIdList();
68
69 LevelSetIdentifierType kk;
70 LevelSetPointer levelSet;
71 LevelSetOutputRealType value;
72
73 auto idListIt = idList->begin();
74 while( idListIt != idList->end() )
75 {
76 //! \todo Fix me for string identifiers
77 kk = *idListIt - 1;
78 if( kk != this->m_CurrentLevelSetId )
79 {
80 levelSet = this->m_LevelSetContainer->GetLevelSet( kk );
81 value = levelSet->Evaluate( iP );
82 prod *= ( NumericTraits< LevelSetOutputRealType >::OneValue() - this->m_Heaviside->Evaluate( -value ) );
83 }
84 ++idListIt;
85 }
86 }
87 }
88 else
89 {
90 LevelSetIdentifierType kk;
91 LevelSetPointer levelSet;
92 LevelSetOutputRealType value;
93
94 typename LevelSetContainerType::Iterator lsIt = this->m_LevelSetContainer->Begin();
95
96 while( lsIt != this->m_LevelSetContainer->End() )
97 {
98 kk = lsIt->GetIdentifier();
99 if( kk != this->m_CurrentLevelSetId )
100 {
101 levelSet = this->m_LevelSetContainer->GetLevelSet( kk );
102 value = levelSet->Evaluate( iP );
103 prod *= ( NumericTraits< LevelSetOutputRealType >::OneValue() - this->m_Heaviside->Evaluate( -value ) );
104 }
105 ++lsIt;
106 }
107 }
108 }
109
110 template< typename TInput, typename TLevelSetContainer >
111 void LevelSetEquationChanAndVeseExternalTerm< TInput, TLevelSetContainer >
UpdatePixel(const LevelSetInputIndexType & iP,const LevelSetOutputRealType & oldValue,const LevelSetOutputRealType & newValue)112 ::UpdatePixel( const LevelSetInputIndexType& iP, const LevelSetOutputRealType & oldValue,
113 const LevelSetOutputRealType & newValue )
114 {
115 // Compute the product factor
116 LevelSetOutputRealType prod;
117
118 this->ComputeProductTerm( iP, prod );
119
120 // For each affected h val: h val = new hval (this will dirty some cvals)
121 InputPixelType input = this->m_Input->GetPixel( iP );
122
123 const LevelSetOutputRealType oldH = this->m_Heaviside->Evaluate( -oldValue );
124 const LevelSetOutputRealType newH = this->m_Heaviside->Evaluate( -newValue );
125 const LevelSetOutputRealType change = oldH - newH;//(1 - newH) - (1 - oldH);
126
127 // Determine the change in the product factor
128 const LevelSetOutputRealType productChange = -( prod * change );
129
130 this->m_TotalH += change;
131 this->m_TotalValue += input * productChange;
132 }
133
134 }
135
136 #endif
137