OTB  6.1.0
Orfeo Toolbox
otbContinuousMinimumMaximumImageCalculator.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 
22 #ifndef otbContinuousMinimumMaximumImageCalculator_txx
23 #define otbContinuousMinimumMaximumImageCalculator_txx
24 
27 
28 namespace otb
29 {
30 
31 /*
32  * Constructor
33  */
34 template<class TInputImage>
37 {
38  m_Image = TInputImage::New();
41  m_IndexOfMinimum.Fill(0);
42  m_IndexOfMaximum.Fill(0);
43  m_RegionSetByUser = false;
44 }
45 
46 /*
47  * Compute Min and Max of m_Image
48  */
49 template<class TInputImage>
50 void
52 ::Compute(void)
53 {
54  if (!m_RegionSetByUser)
55  {
56  m_Region = m_Image->GetRequestedRegion();
57  }
58 
62 
63  while (!it.IsAtEnd())
64  {
65  const RealPixelType value = it.Get();
66  if (value > static_cast<RealPixelType>(m_Maximum))
67  {
68  m_Maximum = static_cast<PixelType>(value);
69  m_IndexOfMaximum = it.GetIndex();
70  }
71  if (value < static_cast<RealPixelType>(m_Minimum))
72  {
73  m_Minimum = static_cast<PixelType>(value);
74  m_IndexOfMinimum = it.GetIndex();
75  }
76  ++it;
77  }
78 
79  IndexType indexNeighbor;
80 
81  { //Continuous Minimum calculation
82  //Compute horizontal offset
83  indexNeighbor[0] = m_IndexOfMinimum[0] - 1;
84  indexNeighbor[1] = m_IndexOfMinimum[1];
85  it.SetIndex(indexNeighbor);
86  const RealPixelType leftValue = it.Get();
87  indexNeighbor[0] = m_IndexOfMinimum[0] + 1;
88  indexNeighbor[1] = m_IndexOfMinimum[1];
89  it.SetIndex(indexNeighbor);
90  const RealPixelType rightValue = it.Get();
91 
92  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Minimum));
93 
94  //Compute vertical offset
95  indexNeighbor[0] = m_IndexOfMinimum[0];
96  indexNeighbor[1] = m_IndexOfMinimum[1] - 1;
97  it.SetIndex(indexNeighbor);
98  const RealPixelType topValue = it.Get();
99  indexNeighbor[0] = m_IndexOfMinimum[0];
100  indexNeighbor[1] = m_IndexOfMinimum[1] + 1;
101  it.SetIndex(indexNeighbor);
102  const RealPixelType bottomValue = it.Get();
103 
104  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Minimum));
105 
106  m_ContinuousIndexOfMinimum[0] = m_IndexOfMinimum[0] + hOffset;
107  m_ContinuousIndexOfMinimum[1] = m_IndexOfMinimum[1] + vOffset;
108  }
109 
110  { //Continuous Maximum calculation
111  //Compute horizontal offset
112  indexNeighbor[0] = m_IndexOfMaximum[0] - 1;
113  indexNeighbor[1] = m_IndexOfMaximum[1];
114  it.SetIndex(indexNeighbor);
115  const RealPixelType leftValue = it.Get();
116  indexNeighbor[0] = m_IndexOfMaximum[0] + 1;
117  indexNeighbor[1] = m_IndexOfMaximum[1];
118  it.SetIndex(indexNeighbor);
119  const RealPixelType rightValue = it.Get();
120 
121  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Maximum));
122 
123  //Compute vertical offset
124  indexNeighbor[0] = m_IndexOfMaximum[0];
125  indexNeighbor[1] = m_IndexOfMaximum[1] - 1;
126  it.SetIndex(indexNeighbor);
127  const RealPixelType topValue = it.Get();
128  indexNeighbor[0] = m_IndexOfMaximum[0];
129  indexNeighbor[1] = m_IndexOfMaximum[1] + 1;
130  it.SetIndex(indexNeighbor);
131  const RealPixelType bottomValue = it.Get();
132 
133  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Maximum));
134 
135  m_ContinuousIndexOfMaximum[0] = m_IndexOfMaximum[0] + hOffset;
136  m_ContinuousIndexOfMaximum[1] = m_IndexOfMaximum[1] + vOffset;
137  }
138 
139 }
140 
141 /*
142  * Compute the minimum intensity value of the image
143  */
144 template<class TInputImage>
145 void
148 {
149  if (!m_RegionSetByUser)
150  {
151  m_Region = m_Image->GetRequestedRegion();
152  }
155 
156  while (!it.IsAtEnd())
157  {
158  const RealPixelType value = it.Get();
159  if (value < static_cast<RealPixelType>(m_Minimum))
160  {
161  m_Minimum = static_cast<PixelType>(value);
162  m_IndexOfMinimum = it.GetIndex();
163  }
164  ++it;
165  }
166 
167  IndexType indexNeighbor;
168 
169  //Compute horizontal offset
170  indexNeighbor[0] = m_IndexOfMinimum[0] - 1;
171  indexNeighbor[1] = m_IndexOfMinimum[1];
172  it.SetIndex(indexNeighbor);
173  const RealPixelType leftValue = it.Get();
174  indexNeighbor[0] = m_IndexOfMinimum[0] + 1;
175  indexNeighbor[1] = m_IndexOfMinimum[1];
176  it.SetIndex(indexNeighbor);
177  const RealPixelType rightValue = it.Get();
178 
179  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Minimum));
180 
181  //Compute vertical offset
182  indexNeighbor[0] = m_IndexOfMinimum[0];
183  indexNeighbor[1] = m_IndexOfMinimum[1] - 1;
184  it.SetIndex(indexNeighbor);
185  const RealPixelType topValue = it.Get();
186  indexNeighbor[0] = m_IndexOfMinimum[0];
187  indexNeighbor[1] = m_IndexOfMinimum[1] + 1;
188  it.SetIndex(indexNeighbor);
189  const RealPixelType bottomValue = it.Get();
190 
191  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Minimum));
192 
193  m_ContinuousIndexOfMinimum[0] = m_IndexOfMinimum[0] + hOffset;
194  m_ContinuousIndexOfMinimum[1] = m_IndexOfMinimum[1] + vOffset;
195 
196 }
197 
198 /*
199  * Compute the maximum intensity value of the image
200  */
201 template<class TInputImage>
202 void
205 {
206  if (!m_RegionSetByUser)
207  {
208  m_Region = m_Image->GetRequestedRegion();
209  }
212 
213  while (!it.IsAtEnd())
214  {
215  const RealPixelType value = it.Get();
216  if (value > static_cast<RealPixelType>(m_Maximum))
217  {
218  m_Maximum = static_cast<PixelType>(value);
219  m_IndexOfMaximum = it.GetIndex();
220  }
221  ++it;
222  }
223 
224  //Solve equations a, b, c
225 // y0 = a*x0^2 + b*x0 +c
226 // y1 = a*x1^2 + b*x1 +c
227 // y2 = a*x2^2 + b*x2 +c
228 //
229 // y0 = a - b +c
230 // y1 = c
231 // y2 = a + b +c
232 //
233 //
234 // Max is at -b/2a
235 // -(y2-y0)/(2*(y0+y2-2y1))
236  IndexType indexNeighbor;
237 
238  //Compute horizontal offset
239  indexNeighbor[0] = m_IndexOfMaximum[0] - 1;
240  indexNeighbor[1] = m_IndexOfMaximum[1];
241  it.SetIndex(indexNeighbor);
242  const RealPixelType leftValue = it.Get();
243  indexNeighbor[0] = m_IndexOfMaximum[0] + 1;
244  indexNeighbor[1] = m_IndexOfMaximum[1];
245  it.SetIndex(indexNeighbor);
246  const RealPixelType rightValue = it.Get();
247 
248  double hOffset = -(rightValue - leftValue) / (2 * (rightValue + leftValue - 2 * m_Maximum));
249 
250  //Compute vertical offset
251  indexNeighbor[0] = m_IndexOfMaximum[0];
252  indexNeighbor[1] = m_IndexOfMaximum[1] - 1;
253  it.SetIndex(indexNeighbor);
254  const RealPixelType topValue = it.Get();
255  indexNeighbor[0] = m_IndexOfMaximum[0];
256  indexNeighbor[1] = m_IndexOfMaximum[1] + 1;
257  it.SetIndex(indexNeighbor);
258  const RealPixelType bottomValue = it.Get();
259 
260  double vOffset = -(bottomValue - topValue) / (2 * (bottomValue + topValue - 2 * m_Maximum));
261 
262  m_ContinuousIndexOfMaximum[0] = m_IndexOfMaximum[0] + hOffset;
263  m_ContinuousIndexOfMaximum[1] = m_IndexOfMaximum[1] + vOffset;
264 
265 }
266 
267 template<class TInputImage>
268 void
270 ::SetRegion(const RegionType& region)
271 {
272  m_Region = region;
273  m_RegionSetByUser = true;
274 }
275 
276 template<class TInputImage>
277 void
279 ::PrintSelf(std::ostream& os, itk::Indent indent) const
280 {
281  Superclass::PrintSelf(os, indent);
282 
283  os << indent << "Minimum: "
284  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(m_Minimum)
285  << std::endl;
286  os << indent << "Maximum: "
287  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(m_Maximum)
288  << std::endl;
289  os << indent << "Index of Minimum: " << m_IndexOfMinimum << std::endl;
290  os << indent << "Index of Maximum: " << m_IndexOfMaximum << std::endl;
291  os << indent << "Continuous Index of Minimum: " << m_ContinuousIndexOfMinimum << std::endl;
292  os << indent << "Continuous Index of Maximum: " << m_ContinuousIndexOfMaximum << std::endl;
293 
294  os << indent << "Image: " << std::endl;
295  m_Image->Print(os, indent.GetNextIndent());
296  os << indent << "Region: " << std::endl;
297  m_Region.Print(os, indent.GetNextIndent());
298  os << indent << "Region set by User: " << m_RegionSetByUser << std::endl;
299 }
300 
301 } // end namespace otb
302 
303 #endif
PixelType Get(void) const
const IndexType & GetIndex() const
void SetIndex(const IndexType &ind)
static ITK_CONSTEXPR_FUNC T max(const T &)
static ITK_CONSTEXPR_FUNC T NonpositiveMin()
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
Indent GetNextIndent()