Orfeo Toolbox  3.16
otbContinuousMinimumMaximumImageCalculator.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 
19 #ifndef __otbContinuousMinimumMaximumImageCalculator_txx
20 #define __otbContinuousMinimumMaximumImageCalculator_txx
21 
22 #include "otbContinuousMinimumMaximumImageCalculator.h"
23 #include "itkImageRegionConstIteratorWithIndex.h"
24 #include "itkNumericTraits.h"
25 
26 namespace otb
27 {
28 
29 /*
30  * Constructor
31  */
32 template<class TInputImage>
33 ContinuousMinimumMaximumImageCalculator<TInputImage>
34 ::ContinuousMinimumMaximumImageCalculator()
35 {
36  m_Image = TInputImage::New();
37  m_Maximum = itk::NumericTraits<PixelType>::NonpositiveMin();
38  m_Minimum = itk::NumericTraits<PixelType>::max();
39  m_IndexOfMinimum.Fill(0);
40  m_IndexOfMaximum.Fill(0);
41  m_RegionSetByUser = false;
42 }
43 
44 /*
45  * Compute Min and Max of m_Image
46  */
47 template<class TInputImage>
48 void
49 ContinuousMinimumMaximumImageCalculator<TInputImage>
50 ::Compute(void)
51 {
52  if (!m_RegionSetByUser)
53  {
54  m_Region = m_Image->GetRequestedRegion();
55  }
56 
57  itk::ImageRegionConstIteratorWithIndex<TInputImage> it(m_Image, m_Region);
58  m_Maximum = itk::NumericTraits<PixelType>::NonpositiveMin();
59  m_Minimum = itk::NumericTraits<PixelType>::max();
60 
61  while (!it.IsAtEnd())
62  {
63  const RealPixelType value = it.Get();
64  if (value > static_cast<RealPixelType>(m_Maximum))
65  {
66  m_Maximum = static_cast<PixelType>(value);
67  m_IndexOfMaximum = it.GetIndex();
68  }
69  if (value < static_cast<RealPixelType>(m_Minimum))
70  {
71  m_Minimum = static_cast<PixelType>(value);
72  m_IndexOfMinimum = it.GetIndex();
73  }
74  ++it;
75  }
76 
77  IndexType indexNeighbor;
78 
79  { //Continuous Minimum calculation
80  //Compute horizontal offset
81  indexNeighbor[0] = m_IndexOfMinimum[0] - 1;
82  indexNeighbor[1] = m_IndexOfMinimum[1];
83  it.SetIndex(indexNeighbor);
84  const RealPixelType leftValue = it.Get();
85  indexNeighbor[0] = m_IndexOfMinimum[0] + 1;
86  indexNeighbor[1] = m_IndexOfMinimum[1];
87  it.SetIndex(indexNeighbor);
88  const RealPixelType rightValue = it.Get();
89 
90  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Minimum));
91 
92  //Compute vertical offset
93  indexNeighbor[0] = m_IndexOfMinimum[0];
94  indexNeighbor[1] = m_IndexOfMinimum[1] - 1;
95  it.SetIndex(indexNeighbor);
96  const RealPixelType topValue = it.Get();
97  indexNeighbor[0] = m_IndexOfMinimum[0];
98  indexNeighbor[1] = m_IndexOfMinimum[1] + 1;
99  it.SetIndex(indexNeighbor);
100  const RealPixelType bottomValue = it.Get();
101 
102  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Minimum));
103 
104  m_ContinuousIndexOfMinimum[0] = m_IndexOfMinimum[0] + hOffset;
105  m_ContinuousIndexOfMinimum[1] = m_IndexOfMinimum[1] + vOffset;
106  }
107 
108  { //Continuous Maximum calculation
109  //Compute horizontal offset
110  indexNeighbor[0] = m_IndexOfMaximum[0] - 1;
111  indexNeighbor[1] = m_IndexOfMaximum[1];
112  it.SetIndex(indexNeighbor);
113  const RealPixelType leftValue = it.Get();
114  indexNeighbor[0] = m_IndexOfMaximum[0] + 1;
115  indexNeighbor[1] = m_IndexOfMaximum[1];
116  it.SetIndex(indexNeighbor);
117  const RealPixelType rightValue = it.Get();
118 
119  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Maximum));
120 
121  //Compute vertical offset
122  indexNeighbor[0] = m_IndexOfMaximum[0];
123  indexNeighbor[1] = m_IndexOfMaximum[1] - 1;
124  it.SetIndex(indexNeighbor);
125  const RealPixelType topValue = it.Get();
126  indexNeighbor[0] = m_IndexOfMaximum[0];
127  indexNeighbor[1] = m_IndexOfMaximum[1] + 1;
128  it.SetIndex(indexNeighbor);
129  const RealPixelType bottomValue = it.Get();
130 
131  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Maximum));
132 
133  m_ContinuousIndexOfMaximum[0] = m_IndexOfMaximum[0] + hOffset;
134  m_ContinuousIndexOfMaximum[1] = m_IndexOfMaximum[1] + vOffset;
135  }
136 
137 }
138 
139 /*
140  * Compute the minimum intensity value of the image
141  */
142 template<class TInputImage>
143 void
144 ContinuousMinimumMaximumImageCalculator<TInputImage>
145 ::ComputeMinimum(void)
146 {
147  if (!m_RegionSetByUser)
148  {
149  m_Region = m_Image->GetRequestedRegion();
150  }
151  itk::ImageRegionConstIteratorWithIndex<TInputImage> it(m_Image, m_Region);
152  m_Minimum = itk::NumericTraits<PixelType>::max();
153 
154  while (!it.IsAtEnd())
155  {
156  const RealPixelType value = it.Get();
157  if (value < static_cast<RealPixelType>(m_Minimum))
158  {
159  m_Minimum = static_cast<PixelType>(value);
160  m_IndexOfMinimum = it.GetIndex();
161  }
162  ++it;
163  }
164 
165  IndexType indexNeighbor;
166 
167  //Compute horizontal offset
168  indexNeighbor[0] = m_IndexOfMinimum[0] - 1;
169  indexNeighbor[1] = m_IndexOfMinimum[1];
170  it.SetIndex(indexNeighbor);
171  const RealPixelType leftValue = it.Get();
172  indexNeighbor[0] = m_IndexOfMinimum[0] + 1;
173  indexNeighbor[1] = m_IndexOfMinimum[1];
174  it.SetIndex(indexNeighbor);
175  const RealPixelType rightValue = it.Get();
176 
177  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Minimum));
178 
179  //Compute vertical offset
180  indexNeighbor[0] = m_IndexOfMinimum[0];
181  indexNeighbor[1] = m_IndexOfMinimum[1] - 1;
182  it.SetIndex(indexNeighbor);
183  const RealPixelType topValue = it.Get();
184  indexNeighbor[0] = m_IndexOfMinimum[0];
185  indexNeighbor[1] = m_IndexOfMinimum[1] + 1;
186  it.SetIndex(indexNeighbor);
187  const RealPixelType bottomValue = it.Get();
188 
189  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Minimum));
190 
191  m_ContinuousIndexOfMinimum[0] = m_IndexOfMinimum[0] + hOffset;
192  m_ContinuousIndexOfMinimum[1] = m_IndexOfMinimum[1] + vOffset;
193 
194 }
195 
196 /*
197  * Compute the maximum intensity value of the image
198  */
199 template<class TInputImage>
200 void
201 ContinuousMinimumMaximumImageCalculator<TInputImage>
202 ::ComputeMaximum(void)
203 {
204  if (!m_RegionSetByUser)
205  {
206  m_Region = m_Image->GetRequestedRegion();
207  }
208  itk::ImageRegionConstIteratorWithIndex<TInputImage> it(m_Image, m_Region);
209  m_Maximum = itk::NumericTraits<PixelType>::NonpositiveMin();
210 
211  while (!it.IsAtEnd())
212  {
213  const RealPixelType value = it.Get();
214  if (value > static_cast<RealPixelType>(m_Maximum))
215  {
216  m_Maximum = static_cast<PixelType>(value);
217  m_IndexOfMaximum = it.GetIndex();
218  }
219  ++it;
220  }
221 
222  //Solve equations a, b, c
223 // y0 = a*x0^2 + b*x0 +c
224 // y1 = a*x1^2 + b*x1 +c
225 // y2 = a*x2^2 + b*x2 +c
226 //
227 // y0 = a - b +c
228 // y1 = c
229 // y2 = a + b +c
230 //
231 //
232 // Max is at -b/2a
233 // -(y2-y0)/(2*(y0+y2-2y1))
234  IndexType indexNeighbor;
235 
236  //Compute horizontal offset
237  indexNeighbor[0] = m_IndexOfMaximum[0] - 1;
238  indexNeighbor[1] = m_IndexOfMaximum[1];
239  it.SetIndex(indexNeighbor);
240  const RealPixelType leftValue = it.Get();
241  indexNeighbor[0] = m_IndexOfMaximum[0] + 1;
242  indexNeighbor[1] = m_IndexOfMaximum[1];
243  it.SetIndex(indexNeighbor);
244  const RealPixelType rightValue = it.Get();
245 
246  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Maximum));
247 
248  //Compute vertical offset
249  indexNeighbor[0] = m_IndexOfMaximum[0];
250  indexNeighbor[1] = m_IndexOfMaximum[1] - 1;
251  it.SetIndex(indexNeighbor);
252  const RealPixelType topValue = it.Get();
253  indexNeighbor[0] = m_IndexOfMaximum[0];
254  indexNeighbor[1] = m_IndexOfMaximum[1] + 1;
255  it.SetIndex(indexNeighbor);
256  const RealPixelType bottomValue = it.Get();
257 
258  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Maximum));
259 
260  m_ContinuousIndexOfMaximum[0] = m_IndexOfMaximum[0] + hOffset;
261  m_ContinuousIndexOfMaximum[1] = m_IndexOfMaximum[1] + vOffset;
262 
263 }
264 
265 template<class TInputImage>
266 void
267 ContinuousMinimumMaximumImageCalculator<TInputImage>
268 ::SetRegion(const RegionType& region)
269 {
270  m_Region = region;
271  m_RegionSetByUser = true;
272 }
273 
274 template<class TInputImage>
275 void
276 ContinuousMinimumMaximumImageCalculator<TInputImage>
277 ::PrintSelf(std::ostream& os, itk::Indent indent) const
278 {
279  Superclass::PrintSelf(os, indent);
280 
281  os << indent << "Minimum: "
282  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(m_Minimum)
283  << std::endl;
284  os << indent << "Maximum: "
285  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(m_Maximum)
286  << std::endl;
287  os << indent << "Index of Minimum: " << m_IndexOfMinimum << std::endl;
288  os << indent << "Index of Maximum: " << m_IndexOfMaximum << std::endl;
289  os << indent << "Continuous Index of Minimum: " << m_ContinuousIndexOfMinimum << std::endl;
290  os << indent << "Continuous Index of Maximum: " << m_ContinuousIndexOfMaximum << std::endl;
291 
292  os << indent << "Image: " << std::endl;
293  m_Image->Print(os, indent.GetNextIndent());
294  os << indent << "Region: " << std::endl;
295  m_Region.Print(os, indent.GetNextIndent());
296  os << indent << "Region set by User: " << m_RegionSetByUser << std::endl;
297 }
298 
299 } // end namespace otb
300 
301 #endif
Generated at Sun Feb 3 2013 00:19:41 for Orfeo Toolbox with doxygen 1.8.1.1