OTB  9.0.0
Orfeo Toolbox
otbCBAMI.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 #ifndef otbCBAMI_h
22 #define otbCBAMI_h
23 
24 #include <vector>
25 #include "itkMath.h"
26 
27 namespace otb
28 {
29 
30 namespace Functor
31 {
32 
33 template <class TInput1, class TInput2, class TOutput>
34 class CBAMI
35 {
36 public:
37  typedef typename std::vector<TOutput> VectorType;
38  typedef typename VectorType::iterator IteratorType;
39  typedef typename std::vector<VectorType> VectorOfVectorType;
40  typedef typename VectorOfVectorType::iterator VecOfVecIteratorType;
41 
43  {
44  }
45  virtual ~CBAMI()
46  {
47  }
48  inline TOutput operator()(const TInput1& itA, const TInput2& itB) const
49  {
50  double epsilon = 0.01;
51  VectorType vecA;
52  VectorType vecB;
53 
54  for (unsigned long pos = 0; pos < itA.Size(); ++pos)
55  {
56 
57  vecA.push_back(static_cast<double>(itA.GetPixel(pos)));
58  vecB.push_back(static_cast<double>(itB.GetPixel(pos)));
59  }
60 
61  normalizeInPlace(vecA);
62  normalizeInPlace(vecB);
63 
64  return static_cast<TOutput>(-std::log(static_cast<double>(PhiMI(vecA, vecB) + epsilon)));
65  }
66 
67 protected:
68  inline void normalizeInPlace(VectorType vx) const
69  {
70 
71  TOutput Ex = 0.0;
72 
73  IteratorType itx;
74 
75  for (itx = vx.begin(); itx < vx.end(); ++itx)
76  {
77  Ex += static_cast<TOutput>(*itx);
78  }
79 
80  Ex /= (vx.size());
81 
82  TOutput Vx = 0.0;
83 
84  for (itx = vx.begin(); itx < vx.end(); ++itx)
85  {
86  Vx += static_cast<TOutput>(std::pow(static_cast<double>((*itx) - Ex), 2));
87  }
88 
89  Vx /= (vx.size());
90 
91  for (itx = vx.begin(); itx < vx.end(); ++itx)
92  {
93  (*itx) = ((*itx) - Ex) / static_cast<TOutput>(std::sqrt(static_cast<double>(Vx)));
94  }
95  }
96  inline TOutput Exyc(VectorType vx, VectorType vy) const
97  {
98 
99  TOutput Exy = 0.0;
100  TOutput Ex = 0.0;
101  TOutput Ey = 0.0;
102 
103  IteratorType itx;
104  IteratorType ity;
105 
106  for (itx = vx.begin(), ity = vy.begin(); itx < vx.end(); ++itx, ++ity)
107  {
108  // Ex += (*itx);
109  // Ey += (*ity);
110  Exy += (*itx) * (*ity);
111  }
112 
113  // Ex /= (vx.size());
114  // Ey /= (vy.size());
115  Exy /= (vx.size());
116 
117  return Exy - Ex * Ey;
118  }
119 
120  inline TOutput Exyztc(VectorType vx, VectorType vy, VectorType vz, VectorType vt) const
121  {
122 
123  TOutput Exyzt = 0.0;
124 
125  TOutput Exyz = 0.0;
126  TOutput Exyt = 0.0;
127  TOutput Exzt = 0.0;
128  TOutput Eyzt = 0.0;
129 
130  TOutput Exy = 0.0;
131  TOutput Exz = 0.0;
132  TOutput Ext = 0.0;
133  TOutput Eyz = 0.0;
134  TOutput Eyt = 0.0;
135  TOutput Ezt = 0.0;
136 
137  TOutput Ex = 0.0;
138  TOutput Ey = 0.0;
139  TOutput Ez = 0.0;
140  TOutput Et = 0.0;
141 
142  IteratorType itx;
143  IteratorType ity;
144  IteratorType itz;
145  IteratorType itt;
146 
147  for (itx = vx.begin(), ity = vy.begin(), itz = vz.begin(), itt = vt.begin(); itx < vx.end(); ++itx, ++ity, itz++, itt++)
148  {
149  // Ex += (*itx);
150  // Ey += (*ity);
151  // Ez += (*itz);
152  // Et += (*itt);
153 
154  Exy += (*itx) * (*ity);
155  Exz += (*itx) * (*itz);
156  Ext += (*itx) * (*itt);
157  Eyz += (*ity) * (*itz);
158  Eyt += (*ity) * (*itt);
159  Ezt += (*itz) * (*itt);
160 
161  Exyz += (*itx) * (*ity) * (*itz);
162  Exyt += (*itx) * (*ity) * (*itt);
163  Exzt += (*itx) * (*itz) * (*itt);
164  Eyzt += (*ity) * (*itz) * (*itt);
165 
166  Exyzt += (*itx) * (*ity) * (*itz) * (*itt);
167  }
168 
169  /*Ex /= (vx.size());
170  Ey /= (vx.size());
171  Ez /= (vx.size());
172  Et /= (vx.size()); */
173 
174  Exy /= (vx.size());
175  Exz /= (vx.size());
176  Ext /= (vx.size());
177  Eyz /= (vx.size());
178  Eyt /= (vx.size());
179  Ezt /= (vx.size());
180 
181  Exyz /= (vx.size());
182  Exyt /= (vx.size());
183  Exzt /= (vx.size());
184  Eyzt /= (vx.size());
185 
186  TOutput result = Exyzt - Exyz * Et - Exyt * Ez - Exzt * Ey - Eyzt * Ex + Exy * Ez * Et + Exz * Et * Ey + Ext * Ey * Ez + Eyz * Et * Ex + Eyt * Ex * Ez +
187  Ezt * Ex * Ey - 3 * Ex * Ey * Ez * Et;
188 
189  return result;
190  }
191 
192  inline TOutput Rxy(const VectorType &va, const VectorType &vb) const
193  {
194 
195  return Exyc(va, vb);
196  }
197 
198  inline TOutput Qxijkl(const VectorType &va, const VectorType &vb, const VectorType &vc, const VectorType &vd) const
199  {
200  // IteratorType ita;
201  // IteratorType itb;
202  // IteratorType itc;
203  // IteratorType itd;
204 
205  TOutput Eabcd_c = Exyztc(va, vb, vc, vd);
206 
207  TOutput Eab_c = Exyc(va, vb);
208  TOutput Eac_c = Exyc(va, vc);
209  TOutput Ead_c = Exyc(va, vd);
210  TOutput Ecd_c = Exyc(vc, vd);
211  TOutput Ebd_c = Exyc(vb, vd);
212  TOutput Ebc_c = Exyc(vb, vc);
213 
214  return Eabcd_c - Eab_c * Ecd_c - Eac_c * Ebd_c - Ead_c * Ebc_c;
215  }
216 
217  inline TOutput PhiMI(VectorType v1, VectorType v2) const
218  {
219 
220  VectorOfVectorType donnees;
221  donnees.push_back(v1);
222  donnees.push_back(v2);
223 
228 
229  TOutput termeR = 0.0;
230  TOutput termeQ = 0.0;
231 
232  for (iti = donnees.begin(); iti < donnees.end(); ++iti)
233  for (itj = donnees.begin(); itj < donnees.end(); ++itj)
234  {
235  if (iti != itj)
236  termeR += static_cast<TOutput>(std::pow(static_cast<double>(Rxy((*iti), (*itj))), 2));
237 
238  for (itk = donnees.begin(); itk < donnees.end(); ++itk)
239  for (itl = donnees.begin(); itl < donnees.end(); itl++)
240  {
241  if ((iti != itj) || (iti != itk) || (iti != itl))
242  termeQ += static_cast<TOutput>(std::pow(static_cast<double>(Qxijkl((*iti), (*itj), (*itk), (*itl))), 2));
243  }
244  }
245 
246  return 1.0 / 4.0 * termeR + 1.0 / 48.0 * termeQ;
247  }
248 };
249 }
250 }
251 
252 #endif
otb::Functor::CBAMI::VectorType
std::vector< TOutput > VectorType
Definition: otbCBAMI.h:37
otb::Functor::CBAMI::IteratorType
VectorType::iterator IteratorType
Definition: otbCBAMI.h:38
otb::Functor::CBAMI::PhiMI
TOutput PhiMI(VectorType v1, VectorType v2) const
Definition: otbCBAMI.h:217
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::Functor::CBAMI::Exyc
TOutput Exyc(VectorType vx, VectorType vy) const
Definition: otbCBAMI.h:96
otb::Functor::CBAMI::Qxijkl
TOutput Qxijkl(const VectorType &va, const VectorType &vb, const VectorType &vc, const VectorType &vd) const
Definition: otbCBAMI.h:198
otb::Functor::CBAMI::CBAMI
CBAMI()
Definition: otbCBAMI.h:42
otb::Functor::CBAMI::Exyztc
TOutput Exyztc(VectorType vx, VectorType vy, VectorType vz, VectorType vt) const
Definition: otbCBAMI.h:120
otb::Functor::CBAMI::normalizeInPlace
void normalizeInPlace(VectorType vx) const
Definition: otbCBAMI.h:68
otb::Functor::CBAMI::Rxy
TOutput Rxy(const VectorType &va, const VectorType &vb) const
Definition: otbCBAMI.h:192
otb::Functor::CBAMI
Definition: otbCBAMI.h:34
otb::Functor::CBAMI::operator()
TOutput operator()(const TInput1 &itA, const TInput2 &itB) const
Definition: otbCBAMI.h:48
itk
Definition: otbNoDataHelper.h:31
otb::Functor::CBAMI::VectorOfVectorType
std::vector< VectorType > VectorOfVectorType
Definition: otbCBAMI.h:39
otb::Functor::CBAMI::VecOfVecIteratorType
VectorOfVectorType::iterator VecOfVecIteratorType
Definition: otbCBAMI.h:40
otb::Functor::CBAMI::~CBAMI
virtual ~CBAMI()
Definition: otbCBAMI.h:45