ITK  4.13.0
Insight Segmentation and Registration Toolkit
itkParabolicMorphUtils.h
Go to the documentation of this file.
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 itkParabolicMorphUtils_h
19 #define itkParabolicMorphUtils_h
20 
21 #include <itkArray.h>
22 
23 #include "itkProgressReporter.h"
24 namespace itk
25 {
26 // contact point algorithm
27 template< class LineBufferType, class RealType, bool doDilate >
28 void DoLineCP(LineBufferType & LineBuf, LineBufferType & tmpLineBuf,
29  const RealType magnitude, const RealType m_Extreme)
30 {
31  // contact point algorithm
32  long koffset = 0, newcontact = 0; // how far away the search starts.
33 
34  const long LineLength = LineBuf.size();
35 
36  // negative half of the parabola
37  for ( long pos = 0; pos < LineLength; pos++ )
38  {
39  RealType BaseVal = (RealType)m_Extreme; // the base value for
40  // comparison
41  for ( long krange = koffset; krange <= 0; krange++ )
42  {
43  // difference needs to be paramaterised
44  RealType T = LineBuf[pos + krange] - magnitude * krange * krange;
45  // switch on template parameter - hopefully gets optimized away.
46  if ( doDilate ? ( T >= BaseVal ) : ( T <= BaseVal ) )
47  {
48  BaseVal = T;
49  newcontact = krange;
50  }
51  }
52  tmpLineBuf[pos] = BaseVal;
53  koffset = newcontact - 1;
54  }
55  // positive half of parabola
56  koffset = newcontact = 0;
57  for ( long pos = LineLength - 1; pos >= 0; pos-- )
58  {
59  RealType BaseVal = (RealType)m_Extreme; // the base value for comparison
60  for ( long krange = koffset; krange >= 0; krange-- )
61  {
62  RealType T = tmpLineBuf[pos + krange] - magnitude * krange * krange;
63  if ( doDilate ? ( T >= BaseVal ) : ( T <= BaseVal ) )
64  {
65  BaseVal = T;
66  newcontact = krange;
67  }
68  }
69  LineBuf[pos] = BaseVal;
70  koffset = newcontact + 1;
71  }
72 }
73 
74 // intersection algorithm
75 // This algorithm has been described a couple of times. First by van
76 // den Boomgaard and more recently by Felzenszwalb and Huttenlocher,
77 // in the context of generalized distance transform
78 template< class LineBufferType, class IndexBufferType,
79  class EnvBufferType, class RealType, bool doDilate >
80 void DoLineIntAlg(LineBufferType & LineBuf, EnvBufferType & F,
81  IndexBufferType & v, EnvBufferType & z,
82  const RealType magnitude)
83 {
84  int k; /* Index of rightmost parabola in lower envelope */
85  /* Locations of parabolas in lower envelope */
86  /* these need to be int, rather than unsigned, to make boundary
87  conditions easy to test */
88 
89  // v stores locations of parabolas in the lower envelope
90  // z stores thr location of boundaries between parabolas
91 
92  // I've gone nuts with the static casts etc, because I seemed to
93  // have strange behaviour when I didn't do this. Also managed to get
94  // rid of all the warnings by sticking to size_t and equivalents.
95  RealType s;
96 
97  /* holds precomputed scale*f(q) + q^2 for speedup */
98 // LineBufferType F(LineBuf.size());
99 
100  // initialize
101  k = 0;
102  v[0] = 0;
104  z[1] = NumericTraits< int >::max();
105  F[0] = LineBuf[0] / magnitude;
106  const size_t N( LineBuf.size() );
107 
108  for ( size_t q = 1; q < N; q++ ) /* main loop */
109  {
110  if ( doDilate )
111  {
112  /* precompute f(q) + q^2 for speedup */
113  F[q] = ( LineBuf[q] / magnitude ) - ( static_cast< RealType >( q ) * static_cast< RealType >( q ) );
114  k++;
115  do
116  {
117  /* remove last parabola from surface */
118  k--;
119  /* compute intersection */
120  s = ( F[q] - F[v[k]] ) / ( 2.0 * ( v[k] - static_cast< RealType >( q ) ) );
121  }
122  while ( s <= z[k] );
123  /* bump k to add new parabola */
124  k++;
125  }
126  else
127  {
128  /* precompute f(q) + q^2 for speedup */
129  F[q] = ( LineBuf[q] / magnitude )
130  + ( static_cast< RealType >( q ) * static_cast< RealType >( q ) );
131  k++;
132  do
133  {
134  /* remove last parabola from surface */
135  k--;
136  /* compute intersection */
137  s = ( F[q] - F[v[k]] ) / ( 2.0 * ( static_cast< RealType >( q ) - v[k] ) );
138  }
139  while ( s <= z[k] );
140  /* bump k to add new parabola */
141  k++;
142  }
143  v[k] = q;
144  z[k] = s;
145  itkAssertInDebugAndIgnoreInReleaseMacro( (size_t)( k + 1 ) <= N );
146  z[k + 1] = NumericTraits< int >::max();
147  } /* for q */
148  /* now reconstruct output */
149  if ( doDilate )
150  {
151  k = 0;
152  for ( size_t q = 0; q < N; q++ )
153  {
154  while ( z[k + 1] < static_cast< typename IndexBufferType::ValueType >( q ) )
155  {
156  k++;
157  }
158  itkAssertInDebugAndIgnoreInReleaseMacro(static_cast< size_t >( v[k] ) < N);
159  itkAssertInDebugAndIgnoreInReleaseMacro(static_cast< size_t >( v[k] ) >= 0);
160  LineBuf[q] =
161  static_cast< RealType >( ( F[v[k]]
162  - ( static_cast< RealType >( q )
163  * ( static_cast< RealType >( q ) - 2 * v[k] ) ) ) * magnitude );
164  }
165  }
166  else
167  {
168  k = 0;
169  for ( size_t q = 0; q < N; q++ )
170  {
171  while ( z[k + 1] < static_cast< typename IndexBufferType::ValueType >( q ) )
172  {
173  k++;
174  }
175  itkAssertInDebugAndIgnoreInReleaseMacro(static_cast< size_t >( v[k] ) < N);
176  itkAssertInDebugAndIgnoreInReleaseMacro(static_cast< size_t >( v[k] ) >= 0);
177  LineBuf[q] =
178  ( ( static_cast< RealType >( q ) * ( static_cast< RealType >( q ) - 2 * v[k] ) + F[v[k]] ) * magnitude );
179  }
180  }
181 }
182 
183 template< class TInIter, class TOutIter, class RealType,
184  class OutputPixelType, bool doDilate >
185 void doOneDimension(TInIter & inputIterator, TOutIter & outputIterator,
186  ProgressReporter & progress,
187  const long LineLength,
188  const unsigned direction,
189  const int m_MagnitudeSign,
190  const bool m_UseImageSpacing,
191  const RealType m_Extreme,
192  const RealType image_scale,
193  const RealType Sigma,
194  int ParabolicAlgorithmChoice)
195 {
196  enum ParabolicAlgorithm {
197  NOCHOICE = 0, // decices based on scale - experimental
198  CONTACTPOINT = 1, // sometimes faster at low scale
199  INTERSECTION = 2 // default
200  };
201 
202 // typedef typename std::vector<RealType> LineBufferType;
203 
204  // message from M.Starring suggested performance gain using Array
205  // instead of std::vector.
206  typedef typename itk::Array< RealType > LineBufferType;
207  RealType iscale = 1.0;
208  if ( m_UseImageSpacing )
209  {
210  iscale = image_scale;
211  }
212  if ( ParabolicAlgorithmChoice == NOCHOICE )
213  {
214  // both set to true or false - use scale to figure it out
215  if ( ( 2.0 * Sigma ) < 0.2 )
216  {
217  ParabolicAlgorithmChoice = CONTACTPOINT;
218  }
219  else
220  {
221  ParabolicAlgorithmChoice = INTERSECTION;
222  }
223  }
224 
225  if ( ParabolicAlgorithmChoice == CONTACTPOINT )
226  {
227  // using the contact point algorithm
228 
229 // const RealType magnitude = m_MagnitudeSign * 1.0/(2.0 *
230 // Sigma/(iscale*iscale));
231  // restructure equation to reduce numerical error
232  const RealType magnitudeCP = ( m_MagnitudeSign * iscale * iscale ) / ( 2.0 * Sigma );
233 
234  LineBufferType LineBuf(LineLength);
235  LineBufferType tmpLineBuf(LineLength);
236  inputIterator.SetDirection(direction);
237  outputIterator.SetDirection(direction);
238  inputIterator.GoToBegin();
239  outputIterator.GoToBegin();
240 
241  unsigned count = 0;
242  while ( !inputIterator.IsAtEnd() && !outputIterator.IsAtEnd() )
243  {
244  // process this direction
245  // fetch the line into the buffer - this methodology is like
246  // the gaussian filters
247 
248  unsigned int i = 0;
249  while ( !inputIterator.IsAtEndOfLine() )
250  {
251  LineBuf[i++] = static_cast< RealType >( inputIterator.Get() );
252  ++inputIterator;
253  }
254 
255  DoLineCP< LineBufferType, RealType, doDilate >(LineBuf, tmpLineBuf, magnitudeCP, m_Extreme);
256  // copy the line back
257  unsigned int j = 0;
258  while ( !outputIterator.IsAtEndOfLine() )
259  {
260  outputIterator.Set( static_cast< OutputPixelType >( LineBuf[j++] ) );
261  ++outputIterator;
262  }
263 
264  ++count;
265  // now onto the next line
266  inputIterator.NextLine();
267  outputIterator.NextLine();
268  progress.CompletedPixel();
269  }
270  }
271  else
272  {
273  // using the Intersection algorithm
274  typedef typename itk::Array< int > IndexBufferType;
275 
276  const RealType magnitudeInt = ( iscale * iscale ) / ( 2.0 * Sigma );
277  LineBufferType LineBuf(LineLength);
278  LineBufferType Fbuf(LineLength);
279  IndexBufferType Vbuf(LineLength);
280  LineBufferType Zbuf(LineLength + 1);
281 
282  inputIterator.SetDirection(direction);
283  outputIterator.SetDirection(direction);
284  inputIterator.GoToBegin();
285  outputIterator.GoToBegin();
286 
287  unsigned count = 0;
288  while ( !inputIterator.IsAtEnd() && !outputIterator.IsAtEnd() )
289  {
290  // process this direction
291  // fetch the line into the buffer - this methodology is like
292  // the gaussian filters
293 
294  unsigned int i = 0;
295  while ( !inputIterator.IsAtEndOfLine() )
296  {
297  LineBuf[i++] = static_cast< RealType >( inputIterator.Get() );
298  ++inputIterator;
299  }
300  DoLineIntAlg< LineBufferType, IndexBufferType,
301  LineBufferType, RealType, doDilate >(LineBuf, Fbuf, Vbuf, Zbuf, magnitudeInt);
302  // copy the line back
303  unsigned int j = 0;
304  while ( !outputIterator.IsAtEndOfLine() )
305  {
306  outputIterator.Set( static_cast< OutputPixelType >( LineBuf[j++] ) );
307  ++outputIterator;
308  }
309 
310  ++count;
311  // now onto the next line
312  inputIterator.NextLine();
313  outputIterator.NextLine();
314  progress.CompletedPixel();
315  }
316  }
317 }
318 }
319 #endif
void DoLineCP(LineBufferType &LineBuf, LineBufferType &tmpLineBuf, const RealType magnitude, const RealType m_Extreme)
void doOneDimension(TInIter &inputIterator, TOutIter &outputIterator, ProgressReporter &progress, const long LineLength, const unsigned direction, const int m_MagnitudeSign, const bool m_UseImageSpacing, const RealType m_Extreme, const RealType image_scale, const RealType Sigma, int ParabolicAlgorithmChoice)
void DoLineIntAlg(LineBufferType &LineBuf, EnvBufferType &F, IndexBufferType &v, EnvBufferType &z, const RealType magnitude)
static ITK_CONSTEXPR_FUNC T max(const T &)
Implements progress tracking for a filter.
static ITK_CONSTEXPR_FUNC T NonpositiveMin()