OTB  5.0.0
Orfeo Toolbox
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 
24 
25 namespace otb
26 {
27 
28 /*
29  * Constructor
30  */
31 template<class TInputImage>
34 {
35  m_Image = TInputImage::New();
38  m_IndexOfMinimum.Fill(0);
39  m_IndexOfMaximum.Fill(0);
40  m_RegionSetByUser = false;
41 }
42 
43 /*
44  * Compute Min and Max of m_Image
45  */
46 template<class TInputImage>
47 void
49 ::Compute(void)
50 {
51  if (!m_RegionSetByUser)
52  {
53  m_Region = m_Image->GetRequestedRegion();
54  }
55 
59 
60  while (!it.IsAtEnd())
61  {
62  const RealPixelType value = it.Get();
63  if (value > static_cast<RealPixelType>(m_Maximum))
64  {
65  m_Maximum = static_cast<PixelType>(value);
66  m_IndexOfMaximum = it.GetIndex();
67  }
68  if (value < static_cast<RealPixelType>(m_Minimum))
69  {
70  m_Minimum = static_cast<PixelType>(value);
71  m_IndexOfMinimum = it.GetIndex();
72  }
73  ++it;
74  }
75 
76  IndexType indexNeighbor;
77 
78  { //Continuous Minimum calculation
79  //Compute horizontal offset
80  indexNeighbor[0] = m_IndexOfMinimum[0] - 1;
81  indexNeighbor[1] = m_IndexOfMinimum[1];
82  it.SetIndex(indexNeighbor);
83  const RealPixelType leftValue = it.Get();
84  indexNeighbor[0] = m_IndexOfMinimum[0] + 1;
85  indexNeighbor[1] = m_IndexOfMinimum[1];
86  it.SetIndex(indexNeighbor);
87  const RealPixelType rightValue = it.Get();
88 
89  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Minimum));
90 
91  //Compute vertical offset
92  indexNeighbor[0] = m_IndexOfMinimum[0];
93  indexNeighbor[1] = m_IndexOfMinimum[1] - 1;
94  it.SetIndex(indexNeighbor);
95  const RealPixelType topValue = it.Get();
96  indexNeighbor[0] = m_IndexOfMinimum[0];
97  indexNeighbor[1] = m_IndexOfMinimum[1] + 1;
98  it.SetIndex(indexNeighbor);
99  const RealPixelType bottomValue = it.Get();
100 
101  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Minimum));
102 
103  m_ContinuousIndexOfMinimum[0] = m_IndexOfMinimum[0] + hOffset;
104  m_ContinuousIndexOfMinimum[1] = m_IndexOfMinimum[1] + vOffset;
105  }
106 
107  { //Continuous Maximum calculation
108  //Compute horizontal offset
109  indexNeighbor[0] = m_IndexOfMaximum[0] - 1;
110  indexNeighbor[1] = m_IndexOfMaximum[1];
111  it.SetIndex(indexNeighbor);
112  const RealPixelType leftValue = it.Get();
113  indexNeighbor[0] = m_IndexOfMaximum[0] + 1;
114  indexNeighbor[1] = m_IndexOfMaximum[1];
115  it.SetIndex(indexNeighbor);
116  const RealPixelType rightValue = it.Get();
117 
118  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Maximum));
119 
120  //Compute vertical offset
121  indexNeighbor[0] = m_IndexOfMaximum[0];
122  indexNeighbor[1] = m_IndexOfMaximum[1] - 1;
123  it.SetIndex(indexNeighbor);
124  const RealPixelType topValue = it.Get();
125  indexNeighbor[0] = m_IndexOfMaximum[0];
126  indexNeighbor[1] = m_IndexOfMaximum[1] + 1;
127  it.SetIndex(indexNeighbor);
128  const RealPixelType bottomValue = it.Get();
129 
130  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Maximum));
131 
132  m_ContinuousIndexOfMaximum[0] = m_IndexOfMaximum[0] + hOffset;
133  m_ContinuousIndexOfMaximum[1] = m_IndexOfMaximum[1] + vOffset;
134  }
135 
136 }
137 
138 /*
139  * Compute the minimum intensity value of the image
140  */
141 template<class TInputImage>
142 void
145 {
146  if (!m_RegionSetByUser)
147  {
148  m_Region = m_Image->GetRequestedRegion();
149  }
152 
153  while (!it.IsAtEnd())
154  {
155  const RealPixelType value = it.Get();
156  if (value < static_cast<RealPixelType>(m_Minimum))
157  {
158  m_Minimum = static_cast<PixelType>(value);
159  m_IndexOfMinimum = it.GetIndex();
160  }
161  ++it;
162  }
163 
164  IndexType indexNeighbor;
165 
166  //Compute horizontal offset
167  indexNeighbor[0] = m_IndexOfMinimum[0] - 1;
168  indexNeighbor[1] = m_IndexOfMinimum[1];
169  it.SetIndex(indexNeighbor);
170  const RealPixelType leftValue = it.Get();
171  indexNeighbor[0] = m_IndexOfMinimum[0] + 1;
172  indexNeighbor[1] = m_IndexOfMinimum[1];
173  it.SetIndex(indexNeighbor);
174  const RealPixelType rightValue = it.Get();
175 
176  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Minimum));
177 
178  //Compute vertical offset
179  indexNeighbor[0] = m_IndexOfMinimum[0];
180  indexNeighbor[1] = m_IndexOfMinimum[1] - 1;
181  it.SetIndex(indexNeighbor);
182  const RealPixelType topValue = it.Get();
183  indexNeighbor[0] = m_IndexOfMinimum[0];
184  indexNeighbor[1] = m_IndexOfMinimum[1] + 1;
185  it.SetIndex(indexNeighbor);
186  const RealPixelType bottomValue = it.Get();
187 
188  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Minimum));
189 
190  m_ContinuousIndexOfMinimum[0] = m_IndexOfMinimum[0] + hOffset;
191  m_ContinuousIndexOfMinimum[1] = m_IndexOfMinimum[1] + vOffset;
192 
193 }
194 
195 /*
196  * Compute the maximum intensity value of the image
197  */
198 template<class TInputImage>
199 void
202 {
203  if (!m_RegionSetByUser)
204  {
205  m_Region = m_Image->GetRequestedRegion();
206  }
209 
210  while (!it.IsAtEnd())
211  {
212  const RealPixelType value = it.Get();
213  if (value > static_cast<RealPixelType>(m_Maximum))
214  {
215  m_Maximum = static_cast<PixelType>(value);
216  m_IndexOfMaximum = it.GetIndex();
217  }
218  ++it;
219  }
220 
221  //Solve equations a, b, c
222 // y0 = a*x0^2 + b*x0 +c
223 // y1 = a*x1^2 + b*x1 +c
224 // y2 = a*x2^2 + b*x2 +c
225 //
226 // y0 = a - b +c
227 // y1 = c
228 // y2 = a + b +c
229 //
230 //
231 // Max is at -b/2a
232 // -(y2-y0)/(2*(y0+y2-2y1))
233  IndexType indexNeighbor;
234 
235  //Compute horizontal offset
236  indexNeighbor[0] = m_IndexOfMaximum[0] - 1;
237  indexNeighbor[1] = m_IndexOfMaximum[1];
238  it.SetIndex(indexNeighbor);
239  const RealPixelType leftValue = it.Get();
240  indexNeighbor[0] = m_IndexOfMaximum[0] + 1;
241  indexNeighbor[1] = m_IndexOfMaximum[1];
242  it.SetIndex(indexNeighbor);
243  const RealPixelType rightValue = it.Get();
244 
245  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Maximum));
246 
247  //Compute vertical offset
248  indexNeighbor[0] = m_IndexOfMaximum[0];
249  indexNeighbor[1] = m_IndexOfMaximum[1] - 1;
250  it.SetIndex(indexNeighbor);
251  const RealPixelType topValue = it.Get();
252  indexNeighbor[0] = m_IndexOfMaximum[0];
253  indexNeighbor[1] = m_IndexOfMaximum[1] + 1;
254  it.SetIndex(indexNeighbor);
255  const RealPixelType bottomValue = it.Get();
256 
257  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Maximum));
258 
259  m_ContinuousIndexOfMaximum[0] = m_IndexOfMaximum[0] + hOffset;
260  m_ContinuousIndexOfMaximum[1] = m_IndexOfMaximum[1] + vOffset;
261 
262 }
263 
264 template<class TInputImage>
265 void
267 ::SetRegion(const RegionType& region)
268 {
269  m_Region = region;
270  m_RegionSetByUser = true;
271 }
272 
273 template<class TInputImage>
274 void
276 ::PrintSelf(std::ostream& os, itk::Indent indent) const
277 {
278  Superclass::PrintSelf(os, indent);
279 
280  os << indent << "Minimum: "
281  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(m_Minimum)
282  << std::endl;
283  os << indent << "Maximum: "
284  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(m_Maximum)
285  << std::endl;
286  os << indent << "Index of Minimum: " << m_IndexOfMinimum << std::endl;
287  os << indent << "Index of Maximum: " << m_IndexOfMaximum << std::endl;
288  os << indent << "Continuous Index of Minimum: " << m_ContinuousIndexOfMinimum << std::endl;
289  os << indent << "Continuous Index of Maximum: " << m_ContinuousIndexOfMaximum << std::endl;
290 
291  os << indent << "Image: " << std::endl;
292  m_Image->Print(os, indent.GetNextIndent());
293  os << indent << "Region: " << std::endl;
294  m_Region.Print(os, indent.GetNextIndent());
295  os << indent << "Region set by User: " << m_RegionSetByUser << std::endl;
296 }
297 
298 } // end namespace otb
299 
300 #endif
PixelType Get(void) const
const IndexType & GetIndex() const
void SetIndex(const IndexType &ind)
static T max(const T &)
void PrintSelf(std::ostream &os, itk::Indent indent) const
static T NonpositiveMin()
Indent GetNextIndent()