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 #ifndef itkOrthogonalSwath2DPathFilter_hxx
19 #define itkOrthogonalSwath2DPathFilter_hxx
20
21 #include "itkOrthogonalSwath2DPathFilter.h"
22 #include "itkMath.h"
23 #include "itkNumericTraits.h"
24
25 namespace itk
26 {
27 /**
28 * Constructor
29 */
30 template< typename TParametricPath, typename TSwathMeritImage >
31 OrthogonalSwath2DPathFilter< TParametricPath, TSwathMeritImage >
OrthogonalSwath2DPathFilter()32 ::OrthogonalSwath2DPathFilter()
33 {
34 SizeType size;
35
36 // Initialize the member variables
37 size[0] = 0;
38 size[1] = 0;
39 m_SwathSize = size;
40 m_StepValues = nullptr;
41 m_MeritValues = nullptr;
42 m_OptimumStepsValues = nullptr;
43 m_FinalOffsetValues = OrthogonalCorrectionTableType::New();
44 }
45
46 /**
47 * Destructor
48 */
49 template< typename TParametricPath, typename TSwathMeritImage >
50 OrthogonalSwath2DPathFilter< TParametricPath, TSwathMeritImage >
~OrthogonalSwath2DPathFilter()51 ::~OrthogonalSwath2DPathFilter()
52 {
53 delete[] m_StepValues;
54 delete[] m_MeritValues;
55 delete[] m_OptimumStepsValues;
56 }
57
58 /**
59 * GenerateData Performs the reflection
60 */
61 template< typename TParametricPath, typename TSwathMeritImage >
62 void
63 OrthogonalSwath2DPathFilter< TParametricPath, TSwathMeritImage >
GenerateData()64 ::GenerateData()
65 {
66 // Get a convenience pointer
67 ImageConstPointer swathMeritImage = this->GetImageInput();
68
69 // Re-initialize the member variables
70 m_SwathSize = swathMeritImage->GetLargestPossibleRegion().GetSize();
71 delete[] m_StepValues;
72 delete[] m_MeritValues;
73 delete[] m_OptimumStepsValues;
74 m_StepValues = new int[m_SwathSize[0] * m_SwathSize[1] * m_SwathSize[1]];
75 m_MeritValues = new double[m_SwathSize[0] * m_SwathSize[1] * m_SwathSize[1]];
76 m_OptimumStepsValues = new int[m_SwathSize[0]];
77 m_FinalOffsetValues->Initialize();
78
79 // Perform the remaining calculations; use dynamic programming
80
81 // current swath column (all previous columns have been fully processed)
82 unsigned int x;
83 // current first row and last row of the swath.
84 unsigned int F, L;
85 // index used to access the processed swath image; filled in with x, F, & L
86 IndexType index;
87
88 // CalcFirstStep (x=0)
89 // Enter the initial merit values
90 index[0] = 0;
91 for ( F = 0; F < m_SwathSize[1]; F++ )
92 {
93 for ( L = 0; L < m_SwathSize[1]; L++ )
94 {
95 if ( F == L )
96 {
97 index[1] = F;
98 MeritValue(F, L, 0) = (double)swathMeritImage->GetPixel(index);
99 StepValue(F, L, 0) = F;
100 }
101 else
102 {
103 MeritValue(F, L, 0) = NumericTraits< double >::NonpositiveMin();
104 StepValue(F, L, 0) = F;
105 }
106 }
107 }
108 // end of double for-loop covering F & L
109
110 // PrepForRemainingSteps
111 for ( F = 0; F < m_SwathSize[1]; F++ )
112 {
113 for ( L = 0; L < m_SwathSize[1]; L++ )
114 {
115 // find merit for x=1
116 if ( itk::Math::abs(F - L) <= 1 )
117 {
118 IndexType index2; // we need a second index here
119 index[0] = 0;
120 index[1] = F;
121 index2[0] = 1;
122 index2[1] = L;
123 // Here we know in advance that Pixel(0,F) =
124 // Max(l=L-1..L+1){Merit(F,l,0)}
125 MeritValue(F, L, 1) = double( swathMeritImage->GetPixel(index)
126 + swathMeritImage->GetPixel(index2) );
127 }
128 else
129 {
130 MeritValue(F, L, 1) = NumericTraits< double >::NonpositiveMin();
131 }
132 // Enter the final step values (x=SWATH_COLUMNS-1)
133 StepValue(F, L, m_SwathSize[0] - 1) = L;
134 }
135 }
136 // end of double for-loop covering F & L
137
138 // CalcRestPath
139 for ( x = 1; x < m_SwathSize[0] - 1; x++ )
140 {
141 for ( F = 0; F < m_SwathSize[1]; F++ )
142 {
143 for ( L = 0; L < m_SwathSize[1]; L++ )
144 {
145 int bestL = FindAndStoreBestErrorStep(x, F, L);
146 index[0] = x + 1;
147 index[1] = L;
148 MeritValue(F, L, x + 1) = MeritValue(F, bestL, x)
149 + double( swathMeritImage->GetPixel(index) );
150 }
151 }
152 }
153 // end of tripple for-loop covering x & F & L
154
155 // Find the best starting and ending points (F & L) for the path
156 int bestF = 0, bestL = 0;
157 double meritTemp, meritMax = NumericTraits< double >::NonpositiveMin();
158 for ( F = 0; F < m_SwathSize[1]; F++ )
159 {
160 for ( L = 0; L < m_SwathSize[1]; L++ )
161 {
162 if ( itk::Math::abs(F - L) <= 1 ) // only accept closed paths
163 {
164 meritTemp = MeritValue(F, L, m_SwathSize[0] - 1);
165 if ( meritTemp > meritMax )
166 {
167 meritMax = meritTemp;
168 bestF = F;
169 bestL = L;
170 }
171 }
172 }
173 }
174 // end of double for-loop covering F & L
175
176 // Fill in the optimum path error-step (orthogonal correction) values
177 m_OptimumStepsValues[m_SwathSize[0] - 1] = bestL;
178 for ( x = m_SwathSize[0] - 2;; x-- )
179 {
180 m_OptimumStepsValues[x] = StepValue(bestF, m_OptimumStepsValues[x + 1], x);
181 if ( 0 == x ) { break; }
182 }
183
184 // Convert from absolute indices to +/- orthogonal offset values
185 for ( x = 0; x < m_SwathSize[0]; x++ )
186 {
187 m_FinalOffsetValues->InsertElement(
188 x, double( m_OptimumStepsValues[x] - int(m_SwathSize[1] / 2) ) );
189 }
190
191 // setup the output path
192 OutputPathPointer outputPtr = this->GetOutput(0);
193 outputPtr->SetOriginalPath( this->GetPathInput() );
194 outputPtr->SetOrthogonalCorrectionTable(m_FinalOffsetValues);
195 }
196
197 template< typename TParametricPath, typename TSwathMeritImage >
198 void
199 OrthogonalSwath2DPathFilter< TParametricPath, TSwathMeritImage >
PrintSelf(std::ostream & os,Indent indent) const200 ::PrintSelf(std::ostream & os, Indent indent) const
201 {
202 Superclass::PrintSelf(os, indent);
203 os << indent << "StepValues: " << m_StepValues << std::endl;
204 os << indent << "MeritValues: " << m_MeritValues << std::endl;
205 os << indent << "OptimumStepsValues: " << m_OptimumStepsValues << std::endl;
206 os << indent << "FinalOffsetValues: " << m_FinalOffsetValues << std::endl;
207 }
208
209 // The next three functions are private helper functions
210 template< typename TParametricPath, typename TSwathMeritImage >
211 unsigned int
212 OrthogonalSwath2DPathFilter< TParametricPath, TSwathMeritImage >
FindAndStoreBestErrorStep(unsigned int x,unsigned int F,unsigned int L)213 ::FindAndStoreBestErrorStep(unsigned int x, unsigned int F, unsigned int L)
214 {
215 unsigned int bestL; // L with largest merit of L and its 2 neighbors L-1 & L+1
216
217 // Handle perimeter boundaries of the vert. gradient image
218 if ( L == 0 )
219 {
220 if ( MeritValue(F, L + 1, x) > MeritValue(F, L, x) )
221 {
222 bestL = L + 1;
223 }
224 else
225 {
226 bestL = L;
227 }
228 }
229 else if ( L == m_SwathSize[1] - 1 )
230 {
231 if ( MeritValue(F, L - 1, x) > MeritValue(F, L, x) )
232 {
233 bestL = L - 1;
234 }
235 else
236 {
237 bestL = L;
238 }
239 }
240 else
241 {
242 // We are now free to consider all 3 cases for bestL
243 if ( MeritValue(F, L + 1, x) > MeritValue(F, L, x)
244 && MeritValue(F, L + 1, x) > MeritValue(F, L - 1, x) )
245 {
246 bestL = L + 1;
247 }
248 else if ( MeritValue(F, L - 1, x) > MeritValue(F, L, x)
249 && MeritValue(F, L - 1, x) > MeritValue(F, L + 1, x) )
250 {
251 bestL = L - 1;
252 }
253 else
254 {
255 bestL = L;
256 }
257 }
258 StepValue(F, L, x) = bestL;
259 return bestL;
260 }
261 } // end namespace itk
262
263 #endif
264