Compare View

switch
from
...
to
 
Commits (3)
include/iota2FeatureExtraction.h
... ... @@ -103,12 +103,12 @@ public:
103 103 m_CopyInputBands{pars.CopyInputBands}
104 104 {
105 105 m_NumberOfDates = m_NumberOfInputComponents/m_ComponentsPerDate;
106   - auto numberOfOutputFeatures = m_NumberOfFeatures-
  106 + const auto numberOfOutputFeatures = m_NumberOfFeatures-
107 107 ((m_RelativeReflectances&&m_RemoveDuplicates)?1:0);
108 108 m_NumberOfOutputComponents = ( numberOfOutputFeatures +
109 109 (m_CopyInputBands?
110 110 m_ComponentsPerDate:0))*m_NumberOfDates;
111   - auto max_index_band = std::max({m_RedIndex, m_NIRIndex, m_SWIRIndex});
  111 + const auto max_index_band = std::max({m_RedIndex, m_NIRIndex, m_SWIRIndex});
112 112 if(max_index_band > m_ComponentsPerDate)
113 113 throw std::domain_error("Band indices and components per date are not coherent.");
114 114 };
... ... @@ -119,8 +119,8 @@ public:
119 119 throw std::domain_error("Pixel size incoherent with number of components per date.");
120 120 PixelType result(m_NumberOfOutputComponents);
121 121 //use std vectors instead of pixels
122   - auto inVec = VectorType(p.GetDataPointer(),
123   - p.GetDataPointer()+p.GetSize());
  122 + const auto inVec = VectorType(p.GetDataPointer(),
  123 + p.GetDataPointer()+p.GetSize());
124 124 //copy the spectral bands
125 125 auto outVec = VectorType(m_NumberOfOutputComponents);
126 126 //copy the input reflectances
... ... @@ -148,6 +148,7 @@ public:
148 148 }
149 149  
150 150 protected:
  151 + inline
151 152 void AddReflectances(const VectorType& inVec, VectorType& outVec)
152 153 {
153 154 if(!m_RelativeReflectances)
... ... @@ -156,17 +157,35 @@ protected:
156 157 }
157 158 else
158 159 {
159   - auto inIt = inVec.cbegin();
160   - auto outIt = outVec.begin();
161   - while(inIt != inVec.cend())
  160 + AddRelativeReflectances(inVec, outVec);
  161 + }
  162 + }
  163 +
  164 + inline
  165 + void AddRelativeReflectances(const VectorType& inVec, VectorType& outVec)
  166 + {
  167 + for(size_t d=0; d<m_NumberOfDates; ++d)
  168 + {
  169 + for(size_t c=0; c<m_ComponentsPerDate; ++c)
162 170 {
163   - (*outIt) = normalized_index(*inIt, *(inIt+m_ReferenceIndex-1));
164   - ++inIt;
165   - ++outIt;
  171 + const size_t date_offset = m_ComponentsPerDate*d;
  172 + const size_t position = c+date_offset;
  173 + const size_t refrefl_position = m_ReferenceIndex-1+date_offset;
  174 + if(position != refrefl_position)
  175 + {
  176 + outVec[position] = normalized_index(inVec[position],
  177 + inVec[refrefl_position]) *
  178 + m_NormalizedIndexFactor;
  179 + }
  180 + else
  181 + {
  182 + outVec[position] = inVec[position];
  183 + }
166 184 }
167 185 }
168 186 }
169   -
  187 +
  188 + inline
170 189 void ComputeFeatures(const VectorType& inVec, VectorType& outVec)
171 190 {
172 191 size_t copyOffset = (m_CopyInputBands?m_NumberOfInputComponents:0);
... ... @@ -188,19 +207,20 @@ protected:
188 207 else
189 208 {
190 209 //compute the features
191   - auto red = *(inIt+m_RedIndex-1);
192   - auto nir = *(inIt+m_NIRIndex-1);
193   - auto swir = *(inIt+m_SWIRIndex-1);
  210 + const auto red = *(inIt+m_RedIndex-1);
  211 + const auto nir = *(inIt+m_NIRIndex-1);
  212 + const auto swir = *(inIt+m_SWIRIndex-1);
194 213 VectorType tmpVec(m_ComponentsPerDate);
195 214 std::transform(inIt, inIt+m_ComponentsPerDate,tmpVec.begin(),
196 215 [](decltype(*inIt)x){ return x*x;});
197   - auto brightness = std::sqrt(std::accumulate(tmpVec.begin(), tmpVec.end(),
198   - ValueType{0}));
  216 + const auto brightness = std::sqrt(std::accumulate(tmpVec.begin(), tmpVec.end(),
  217 + ValueType{0}));
199 218 //append the features
200 219 size_t featureOffset{0};
201   -
  220 + //ndvi
202 221 AddNormalizedIndexMaybe(nir, red, m_RedIndex, featureOffset,
203 222 copyOffset, outVec, date_counter);
  223 + //ndwi
204 224 AddNormalizedIndexMaybe(swir, nir, m_NIRIndex, featureOffset,
205 225 copyOffset, outVec, date_counter);
206 226 outVec[copyOffset+m_NumberOfDates*featureOffset+date_counter] = brightness;
... ... @@ -211,10 +231,11 @@ protected:
211 231 }
212 232 }
213 233  
  234 + inline
214 235 void AddNormalizedIndexMaybe(ValueType refl, ValueType refrefl,
215   - size_t refindex, size_t& featureOffset,
216   - size_t copyOffset, VectorType& outVec,
217   - size_t date_counter)
  236 + size_t refindex, size_t& featureOffset,
  237 + size_t copyOffset, VectorType& outVec,
  238 + size_t date_counter)
218 239 {
219 240 if(!m_RemoveDuplicates || m_ReferenceIndex != refindex)
220 241 {
... ... @@ -249,3 +270,5 @@ protected:
249 270  
250 271  
251 272  
  273 +
  274 +
... ...
test/iota2FeatureExtractionFunctor.cxx
... ... @@ -94,6 +94,7 @@ int fexFunctor(int argc, char * argv[])
94 94 }
95 95  
96 96  
  97 + auto tmpp0 = p[0];
97 98 p[0] = ndv;
98 99 res = func(p);
99 100 ndvi1_res = res[outOffset*nbd];
... ... @@ -102,6 +103,41 @@ int fexFunctor(int argc, char * argv[])
102 103 std::cout << p[0] << "\t" << res[1] << " --> should be no data\n";
103 104 return EXIT_FAILURE;
104 105 }
  106 + p[0] = tmpp0;
105 107  
106   - return EXIT_SUCCESS;
107   -}
  108 +
  109 + // test the relative reflectances
  110 + copyInputBands = true;
  111 + pars.CopyInputBands = copyInputBands;
  112 + pars.RelativeReflectances = true;
  113 + pars.ReferenceIndex = pars.RedIndex;
  114 + pars.RemoveDuplicates = false;
  115 +
  116 + func = iota2::FeatureExtractionFunctor<PixelType>(pars);
  117 + res = func(p);
  118 + ndvi1_res = res[cpd*nbd];
  119 + auto relref_res = res[pars.NIRIndex-1];
  120 + if(ndvi1_res != relref_res)
  121 + {
  122 + std::cout << "Relative reflectance for NIR != ndvi \n";
  123 + return EXIT_FAILURE;
  124 + }
  125 + if(p[pars.ReferenceIndex-1]!=res[pars.ReferenceIndex-1])
  126 + {
  127 + std::cout << "The relative reflectance of the reference index should be the original one\n";
  128 + return EXIT_FAILURE;
  129 + }
  130 +
  131 + //test remove duplicates
  132 + pars.RemoveDuplicates = true;
  133 + func = iota2::FeatureExtractionFunctor<PixelType>(pars);
  134 + res = func(p);
  135 + if(res.GetSize() != (cpd+2)*nbd)
  136 + {
  137 + std::cout << "Remove duplicates yields " << res.GetSize() << " components\n";
  138 + return EXIT_FAILURE;
  139 + }
  140 +
  141 +
  142 + return EXIT_SUCCESS;
  143 + }
... ...