Commit 19f19137cc3b44cc05c8213814319f9edcfbac1d

Authored by Jordi Inglada
1 parent 5ec1482e

ENH: implement atmospherically corrected normalized indices

app/iota2FeatureExtraction.cxx
... ... @@ -91,6 +91,9 @@ private:
91 91 AddParameter(ParameterType_Empty, "keepduplicates", "Keep duplicate relative reflectances (true/false). Default value is false");
92 92 MandatoryOff("keepduplicates");
93 93  
  94 + AddParameter(ParameterType_Empty, "acorfeat", "Apply atmospherically corrected features http://www.cesbio.ups-tlse.fr/multitemp/?p=12746 (true/false). Default value is false");
  95 + MandatoryOff("acorfeat");
  96 +
94 97  
95 98  
96 99 AddRAMParameter();
... ... @@ -144,6 +147,11 @@ private:
144 147 std::cout << " relative index " << pars.ReferenceIndex << " \n";
145 148 }
146 149  
  150 + if(IsParameterEnabled("acorfeat"))
  151 + {
  152 + std::cout << " Atmospherically corrected features \n";
  153 + pars.ACorFeat = true;
  154 + }
147 155  
148 156 auto fef = FeatureExtractionFunctorType(pars);
149 157 m_FeatureExtractionFilter = FeatureExtractionFilterType::New();
... ...
include/iota2FeatureExtraction.h
... ... @@ -64,217 +64,231 @@ private:
64 64 };
65 65  
66 66 template <typename T>
67   -constexpr T normalized_index(T refl, T refrefl, T epsilon=10e-6)
  67 +constexpr T normalized_index(const T refl, const T refrefl,
  68 + const bool acorfeat=false,
  69 + const T acorepsilon=0.05,
  70 + const T epsilon=10e-6)
68 71 {
69   - return std::fabs(refl+refrefl)<epsilon?
  72 + if(!acorfeat)
  73 + {
  74 + return std::fabs(refl+refrefl)<epsilon?
70 75 T{0}:(refl-refrefl)/(refl+refrefl);
  76 + }
  77 + else
  78 + {
  79 + return (refl-refrefl-acorepsilon)/(refl+refrefl+acorepsilon);
  80 + }
71 81 }
72 82  
73   -template <typename PixelType>
74   -class FeatureExtractionFunctor
75   -{
76   -public:
77   - using ValueType = typename PixelType::ValueType;
78   - using VectorType = std::vector<ValueType>;
79   -
80   - struct Parameters {
81   - size_t ComponentsPerDate{1};
82   - size_t RedIndex{3};
83   - size_t NIRIndex{4};
84   - size_t SWIRIndex{5};
85   - bool RelativeReflectances{false};
86   - size_t ReferenceIndex{3};
87   - bool RemoveDuplicates{true};
88   - ValueType NormalizedIndexFactor{1000};
89   - ValueType NoDataValue{-10000};
90   - size_t NumberOfInputComponents{1};
91   - bool CopyInputBands{false};
92   - };
93   - FeatureExtractionFunctor() = default;
94   - FeatureExtractionFunctor(Parameters pars)
95   - : m_ComponentsPerDate{pars.ComponentsPerDate}, m_RedIndex{pars.RedIndex},
96   - m_NIRIndex{pars.NIRIndex}, m_SWIRIndex{pars.SWIRIndex},
97   - m_RelativeReflectances{pars.RelativeReflectances},
98   - m_ReferenceIndex{pars.ReferenceIndex},
99   - m_RemoveDuplicates{pars.RemoveDuplicates},
100   - m_NormalizedIndexFactor{pars.NormalizedIndexFactor},
101   - m_NoDataValue{pars.NoDataValue},
102   - m_NumberOfInputComponents{pars.NumberOfInputComponents},
103   - m_CopyInputBands{pars.CopyInputBands}
  83 + template <typename PixelType>
  84 + class FeatureExtractionFunctor
104 85 {
105   - m_NumberOfDates = m_NumberOfInputComponents/m_ComponentsPerDate;
106   - UpdateNumberOfFeatures();
107   - m_NumberOfOutputComponents = ( m_NumberOfFeatures +
108   - (m_CopyInputBands?
109   - m_ComponentsPerDate:0))*m_NumberOfDates;
110   - const auto max_index_band = std::max({m_RedIndex, m_NIRIndex, m_SWIRIndex});
111   - if(max_index_band > m_ComponentsPerDate)
112   - throw std::domain_error("Band indices and components per date are not coherent.");
113   - };
  86 + public:
  87 + using ValueType = typename PixelType::ValueType;
  88 + using VectorType = std::vector<ValueType>;
114 89  
115   - PixelType operator()(const PixelType& p)
116   - {
117   - if(p.GetSize()%m_ComponentsPerDate != 0)
118   - throw std::domain_error("Pixel size incoherent with number of components per date.");
119   - PixelType result(m_NumberOfOutputComponents);
120   - //use std vectors instead of pixels
121   - const auto inVec = VectorType(p.GetDataPointer(),
122   - p.GetDataPointer()+p.GetSize());
123   - //copy the spectral bands
124   - auto outVec = VectorType(m_NumberOfOutputComponents);
125   - //copy the input reflectances
126   - if(m_CopyInputBands)
127   - {
128   - AddReflectances(inVec, outVec);
129   - }
  90 + struct Parameters {
  91 + size_t ComponentsPerDate{1};
  92 + size_t RedIndex{3};
  93 + size_t NIRIndex{4};
  94 + size_t SWIRIndex{5};
  95 + bool RelativeReflectances{false};
  96 + size_t ReferenceIndex{3};
  97 + bool RemoveDuplicates{true};
  98 + ValueType NormalizedIndexFactor{1000};
  99 + ValueType NoDataValue{-10000};
  100 + size_t NumberOfInputComponents{1};
  101 + bool CopyInputBands{false};
  102 + bool ACorFeat{false};
  103 + };
  104 + FeatureExtractionFunctor() = default;
  105 + FeatureExtractionFunctor(Parameters pars)
  106 + : m_ComponentsPerDate{pars.ComponentsPerDate}, m_RedIndex{pars.RedIndex},
  107 + m_NIRIndex{pars.NIRIndex}, m_SWIRIndex{pars.SWIRIndex},
  108 + m_RelativeReflectances{pars.RelativeReflectances},
  109 + m_ReferenceIndex{pars.ReferenceIndex},
  110 + m_RemoveDuplicates{pars.RemoveDuplicates},
  111 + m_NormalizedIndexFactor{pars.NormalizedIndexFactor},
  112 + m_NoDataValue{pars.NoDataValue},
  113 + m_NumberOfInputComponents{pars.NumberOfInputComponents},
  114 + m_CopyInputBands{pars.CopyInputBands},
  115 + m_ACorFeat{pars.ACorFeat}
  116 + {
  117 + m_NumberOfDates = m_NumberOfInputComponents/m_ComponentsPerDate;
  118 + UpdateNumberOfFeatures();
  119 + m_NumberOfOutputComponents = ( m_NumberOfFeatures +
  120 + (m_CopyInputBands?
  121 + m_ComponentsPerDate:0))*m_NumberOfDates;
  122 + const auto max_index_band = std::max({m_RedIndex, m_NIRIndex, m_SWIRIndex});
  123 + if(max_index_band > m_ComponentsPerDate)
  124 + throw std::domain_error("Band indices and components per date are not coherent.");
  125 + };
130 126  
131   - ComputeFeatures(inVec, outVec);
  127 + PixelType operator()(const PixelType& p)
  128 + {
  129 + if(p.GetSize()%m_ComponentsPerDate != 0)
  130 + throw std::domain_error("Pixel size incoherent with number of components per date.");
  131 + PixelType result(m_NumberOfOutputComponents);
  132 + //use std vectors instead of pixels
  133 + const auto inVec = VectorType(p.GetDataPointer(),
  134 + p.GetDataPointer()+p.GetSize());
  135 + //copy the spectral bands
  136 + auto outVec = VectorType(m_NumberOfOutputComponents);
  137 + //copy the input reflectances
  138 + if(m_CopyInputBands)
  139 + {
  140 + AddReflectances(inVec, outVec);
  141 + }
  142 +
  143 + ComputeFeatures(inVec, outVec);
132 144  
133   - //convert the result to a pixel
134   - for(size_t i=0; i<m_NumberOfOutputComponents; i++)
135   - result[i] = outVec[i];
136   - return result;
137   - }
  145 + //convert the result to a pixel
  146 + for(size_t i=0; i<m_NumberOfOutputComponents; i++)
  147 + result[i] = outVec[i];
  148 + return result;
  149 + }
138 150  
139   - bool operator!=(FeatureExtractionFunctor<PixelType> f)
140   - {
141   - return m_ComponentsPerDate != f.m_ComponentsPerDate;
142   - }
  151 + bool operator!=(FeatureExtractionFunctor<PixelType> f)
  152 + {
  153 + return m_ComponentsPerDate != f.m_ComponentsPerDate;
  154 + }
143 155  
144   - size_t GetNumberOfOutputComponents() const
145   - {
146   - return m_NumberOfOutputComponents;
147   - }
  156 + size_t GetNumberOfOutputComponents() const
  157 + {
  158 + return m_NumberOfOutputComponents;
  159 + }
148 160  
149   -protected:
150   - inline void UpdateNumberOfFeatures()
151   - {
152   - if(m_SWIRIndex==0) --m_NumberOfFeatures;
153   - if((m_RelativeReflectances && m_RemoveDuplicates &&
154   - (m_ReferenceIndex==m_RedIndex || m_ReferenceIndex==m_NIRIndex)))
155   - --m_NumberOfFeatures;
156   - }
  161 + protected:
  162 + inline void UpdateNumberOfFeatures()
  163 + {
  164 + if(m_SWIRIndex==0) --m_NumberOfFeatures;
  165 + if((m_RelativeReflectances && m_RemoveDuplicates &&
  166 + (m_ReferenceIndex==m_RedIndex || m_ReferenceIndex==m_NIRIndex)))
  167 + --m_NumberOfFeatures;
  168 + }
157 169  
158   - inline
159   - void AddReflectances(const VectorType& inVec, VectorType& outVec)
160   - {
161   - if(!m_RelativeReflectances)
162   - {
163   - std::copy(inVec.cbegin(), inVec.cend(), outVec.begin());
164   - }
165   - else
166   - {
167   - AddRelativeReflectances(inVec, outVec);
168   - }
169   - }
  170 + inline
  171 + void AddReflectances(const VectorType& inVec, VectorType& outVec)
  172 + {
  173 + if(!m_RelativeReflectances)
  174 + {
  175 + std::copy(inVec.cbegin(), inVec.cend(), outVec.begin());
  176 + }
  177 + else
  178 + {
  179 + AddRelativeReflectances(inVec, outVec);
  180 + }
  181 + }
170 182  
171   - inline
172   - void AddRelativeReflectances(const VectorType& inVec, VectorType& outVec)
173   - {
174   - for(size_t d=0; d<m_NumberOfDates; ++d)
175   - {
176   - for(size_t c=0; c<m_ComponentsPerDate; ++c)
  183 + inline
  184 + void AddRelativeReflectances(const VectorType& inVec, VectorType& outVec)
  185 + {
  186 + for(size_t d=0; d<m_NumberOfDates; ++d)
177 187 {
178   - const size_t date_offset = m_ComponentsPerDate*d;
179   - const size_t position = c+date_offset;
180   - const size_t refrefl_position = m_ReferenceIndex-1+date_offset;
181   - if(position != refrefl_position)
182   - {
183   - outVec[position] = normalized_index(inVec[position],
184   - inVec[refrefl_position]) *
185   - m_NormalizedIndexFactor;
186   - }
187   - else
  188 + for(size_t c=0; c<m_ComponentsPerDate; ++c)
188 189 {
189   - outVec[position] = inVec[position];
  190 + const size_t date_offset = m_ComponentsPerDate*d;
  191 + const size_t position = c+date_offset;
  192 + const size_t refrefl_position = m_ReferenceIndex-1+date_offset;
  193 + if(position != refrefl_position)
  194 + {
  195 + outVec[position] = normalized_index(inVec[position],
  196 + inVec[refrefl_position],
  197 + m_ACorFeat) *
  198 + m_NormalizedIndexFactor;
  199 + }
  200 + else
  201 + {
  202 + outVec[position] = inVec[position];
  203 + }
190 204 }
191 205 }
192   - }
193   - }
  206 + }
194 207  
195   - inline
196   - void ComputeFeatures(const VectorType& inVec, VectorType& outVec)
197   - {
198   - size_t copyOffset = (m_CopyInputBands?m_NumberOfInputComponents:0);
199   - size_t date_counter{0};
200   - auto inIt = inVec.cbegin();
201   - while(inIt != inVec.cend())
202   - {
203   - //check for invalid values
204   - if(std::any_of(inIt, inIt+m_ComponentsPerDate,
205   - [&](ValueType x)
206   - {
207   - return std::fabs(x - m_NoDataValue)<0.1;
208   - }))
  208 + inline
  209 + void ComputeFeatures(const VectorType& inVec, VectorType& outVec)
  210 + {
  211 + size_t copyOffset = (m_CopyInputBands?m_NumberOfInputComponents:0);
  212 + size_t date_counter{0};
  213 + auto inIt = inVec.cbegin();
  214 + while(inIt != inVec.cend())
209 215 {
210   - for(size_t feat=0; feat<m_NumberOfFeatures; ++feat)
  216 + //check for invalid values
  217 + if(std::any_of(inIt, inIt+m_ComponentsPerDate,
  218 + [&](ValueType x)
  219 + {
  220 + return std::fabs(x - m_NoDataValue)<0.1;
  221 + }))
211 222 {
212   - outVec[copyOffset+m_NumberOfDates*feat+date_counter] = m_NoDataValue;
  223 + for(size_t feat=0; feat<m_NumberOfFeatures; ++feat)
  224 + {
  225 + outVec[copyOffset+m_NumberOfDates*feat+date_counter] = m_NoDataValue;
  226 + }
213 227 }
214   - }
215   - else
216   - {
217   - //compute the features
218   - const auto red = *(inIt+m_RedIndex-1);
219   - const auto nir = *(inIt+m_NIRIndex-1);
220   - const auto swir = *(inIt+(m_SWIRIndex>0?m_SWIRIndex:1)-1);
221   - VectorType tmpVec(m_ComponentsPerDate);
222   - std::transform(inIt, inIt+m_ComponentsPerDate,tmpVec.begin(),
223   - [](decltype(*inIt)x){ return x*x;});
224   - const auto brightness = std::sqrt(std::accumulate(tmpVec.begin(), tmpVec.end(),
225   - ValueType{0}));
226   - //append the features
227   - size_t featureOffset{0};
228   - //ndvi
229   - featureOffset = AddNormalizedIndexMaybe(nir, red, m_RedIndex, featureOffset,
230   - copyOffset, outVec, date_counter);
231   - //ndwi
232   - if(m_SWIRIndex!=0)
  228 + else
233 229 {
234   - featureOffset = AddNormalizedIndexMaybe(swir, nir, m_NIRIndex, featureOffset,
  230 + //compute the features
  231 + const auto red = *(inIt+m_RedIndex-1);
  232 + const auto nir = *(inIt+m_NIRIndex-1);
  233 + const auto swir = *(inIt+(m_SWIRIndex>0?m_SWIRIndex:1)-1);
  234 + VectorType tmpVec(m_ComponentsPerDate);
  235 + std::transform(inIt, inIt+m_ComponentsPerDate,tmpVec.begin(),
  236 + [](decltype(*inIt)x){ return x*x;});
  237 + const auto brightness = std::sqrt(std::accumulate(tmpVec.begin(), tmpVec.end(),
  238 + ValueType{0}));
  239 + //append the features
  240 + size_t featureOffset{0};
  241 + //ndvi
  242 + featureOffset = AddNormalizedIndexMaybe(nir, red, m_RedIndex, featureOffset,
235 243 copyOffset, outVec, date_counter);
  244 + //ndwi
  245 + if(m_SWIRIndex!=0)
  246 + {
  247 + featureOffset = AddNormalizedIndexMaybe(swir, nir, m_NIRIndex, featureOffset,
  248 + copyOffset, outVec, date_counter);
  249 + }
  250 + outVec[copyOffset+m_NumberOfDates*featureOffset+date_counter] = brightness;
236 251 }
237   - outVec[copyOffset+m_NumberOfDates*featureOffset+date_counter] = brightness;
  252 + //move to the next date
  253 + std::advance(inIt, m_ComponentsPerDate);
  254 + ++date_counter;
238 255 }
239   - //move to the next date
240   - std::advance(inIt, m_ComponentsPerDate);
241   - ++date_counter;
242   - }
243   - }
  256 + }
244 257  
245   - inline
246   - size_t AddNormalizedIndexMaybe(ValueType refl, ValueType refrefl,
247   - size_t refindex, size_t featureOffset,
248   - size_t copyOffset, VectorType& outVec,
249   - size_t date_counter)
250   - {
251   - auto result = featureOffset;
252   - if(!(m_RelativeReflectances && m_RemoveDuplicates && m_ReferenceIndex == refindex))
253   - {
254   - outVec[copyOffset+m_NumberOfDates*featureOffset+date_counter] =
255   - normalized_index(refl, refrefl) * m_NormalizedIndexFactor;
256   - ++result;
257   - }
258   - return result;
259   - }
  258 + inline
  259 + size_t AddNormalizedIndexMaybe(ValueType refl, ValueType refrefl,
  260 + size_t refindex, size_t featureOffset,
  261 + size_t copyOffset, VectorType& outVec,
  262 + size_t date_counter)
  263 + {
  264 + auto result = featureOffset;
  265 + if(!(m_RelativeReflectances && m_RemoveDuplicates && m_ReferenceIndex == refindex))
  266 + {
  267 + outVec[copyOffset+m_NumberOfDates*featureOffset+date_counter] =
  268 + normalized_index(refl, refrefl, m_ACorFeat) * m_NormalizedIndexFactor;
  269 + ++result;
  270 + }
  271 + return result;
  272 + }
260 273  
261   - size_t m_ComponentsPerDate;
262   - size_t m_RedIndex;
263   - size_t m_NIRIndex;
264   - size_t m_SWIRIndex;
265   - bool m_RelativeReflectances;
266   - size_t m_ReferenceIndex; //which reflectance is used as reference if
267   - //relative reflectances are used
268   - bool m_RemoveDuplicates; //If relative reflectances, NDVI or NDWI
269   - //may be redundant
270   - ValueType m_NormalizedIndexFactor;
271   - ValueType m_NoDataValue;
272   - size_t m_NumberOfInputComponents;
273   - bool m_CopyInputBands;
274   - size_t m_NumberOfOutputComponents;
275   - size_t m_NumberOfDates;
276   - size_t m_NumberOfFeatures = 3;
277   -};
  274 + size_t m_ComponentsPerDate;
  275 + size_t m_RedIndex;
  276 + size_t m_NIRIndex;
  277 + size_t m_SWIRIndex;
  278 + bool m_RelativeReflectances;
  279 + size_t m_ReferenceIndex; //which reflectance is used as reference if
  280 + //relative reflectances are used
  281 + bool m_RemoveDuplicates; //If relative reflectances, NDVI or NDWI
  282 + //may be redundant
  283 + ValueType m_NormalizedIndexFactor;
  284 + ValueType m_NoDataValue;
  285 + size_t m_NumberOfInputComponents;
  286 + bool m_CopyInputBands;
  287 + bool m_ACorFeat;
  288 + size_t m_NumberOfOutputComponents;
  289 + size_t m_NumberOfDates;
  290 + size_t m_NumberOfFeatures = 3;
  291 + };
278 292 } // end namespace iota2
279 293  
280 294  
... ...