ENH: CalaculateResponse method in ReduceSpectralResponse computation speedup.

Code/Simulation/otbReduceSpectralResponse.txx

@@ -63,6 +63,8 @@ ReduceSpectralResponse<TSpectralResponse , TRSR>
     PrecisionType lambda1;
     PrecisionType lambda2;
 
+
+
     typename VectorPairType::const_iterator it;
     VectorPairType pairs = (m_InputSatRSR->GetRSR())[numBand]->GetResponse();
     it = pairs.begin();
@@ -88,6 +90,7 @@ ReduceSpectralResponse<TSpectralResponse , TRSR>
         inputRSR1 = (*m_InputSpectralResponse)(lambda1);
         inputRSR2 = (*m_InputSpectralResponse)(lambda2);
 
+
         // lambda1 need to be resampled
         /*
 
@@ -126,11 +129,13 @@ ReduceSpectralResponse<TSpectralResponse , TRSR>
 
           PrecisionType lambdaDist = lambdaRSRmin - lambda1;
           PrecisionType ratio = lambdaDist / (lambda2 - lambda1);
-
+          std::cout<<"modif lambda 1 !!! "<< lambda1<<" ";
           lambda1 = lambdaRSRmin;
-
+          std::cout<<" "<<lambda1<<std::endl;
           inputSatRSR1 = ratio * inputSatRSR1 + (1 - ratio) * inputSatRSR2;
 
+          inputRSR1=(*m_InputSpectralResponse)(lambda1);
+
           }
 
         // lambda2 need to be resampled
@@ -170,14 +175,17 @@ ReduceSpectralResponse<TSpectralResponse , TRSR>
             }
           PrecisionType lambdaDist = lambdaRSRmax - lambda1;
           PrecisionType ratio = lambdaDist / (lambda2 - lambda1);
-
+          std::cout<<"modif lambda 2 !!! "<< lambda2<<" ";
           lambda2 = lambdaRSRmax;
-
+          std::cout<<" "<<lambda2<<std::endl;
           inputSatRSR2 = ratio * inputSatRSR1 + (1 - ratio) * inputSatRSR2;
+
+          inputRSR2=(*m_InputSpectralResponse)(lambda2);
+
           }
 
-        response1 = (*m_InputSpectralResponse)(lambda1) * inputSatRSR1;
-        response2 = (*m_InputSpectralResponse)(lambda2) * inputSatRSR2;
+        response1 = inputRSR1 * inputSatRSR1;
+        response2 = inputRSR2 * inputSatRSR2;
 
         ValuePrecisionType rmin = std::min(response1, response2);
         ValuePrecisionType rmax = std::max(response1, response2);
@@ -216,13 +224,17 @@ ReduceSpectralResponse<TSpectralResponse , TRSR>
    //Compute the reduce response for each band of the sensor
    for (unsigned int i=0; i<m_InputSatRSR->GetNbBands(); ++i)
    {
-      PairType pair;
-      //pair.first = center wavelength of the band
-      pair.first=( (this->m_InputSatRSR->GetRSR())[i]->GetInterval().first + (this->m_InputSatRSR->GetRSR())[i]->GetInterval().second );
-      pair.first=pair.first/2.0;
-      pair.second=(*this)(i);
-      m_ReduceResponse->GetResponse().push_back(pair);
-   }
+    m_InputSpectralResponse->SetPosGuessMin((this->m_InputSatRSR->GetRSR())[i]->GetInterval().first);
+    m_InputSpectralResponse->SetUsePosGuess(true);
+    PairType pair;
+    //pair.first = center wavelength of the band
+    pair.first = ((this->m_InputSatRSR->GetRSR())[i]->GetInterval().first
+        + (this->m_InputSatRSR->GetRSR())[i]->GetInterval().second);
+    pair.first = pair.first / 2.0;
+    pair.second = (*this)(i);
+    m_ReduceResponse->GetResponse().push_back(pair);
+    m_InputSpectralResponse->SetUsePosGuess(false);
+    }
 
 }
 

Code/Simulation/otbSpectralResponse.h

@@ -78,13 +78,16 @@ public:
 
   typedef std::pair<TPrecision, TPrecision> IntervalType;
   /** Standard macros */
-  itkNewMacro(Self)
- ; itkTypeMacro(SpectralResponse, DataObject)
- ;
+  itkNewMacro(Self);
+  itkTypeMacro(SpectralResponse, DataObject);
+
+
+  itkSetMacro(SensitivityThreshold, TPrecision);
+  itkGetConstMacro(SensitivityThreshold, TPrecision);
+
+  itkSetMacro(UsePosGuess, bool);
+  itkGetConstMacro(UsePosGuess, bool);
 
-  itkSetMacro(SensitivityThreshold, TPrecision)
- ; itkGetConstMacro(SensitivityThreshold, TPrecision)
- ;
 
   /** Clear the vector of data pairs  */
   virtual bool Clear();
@@ -115,6 +118,8 @@ public:
    */
   inline ValuePrecisionType operator()(const PrecisionType & lambda);
 
+
+
   /** Operator for comparing Pair Lambda/Response
    * Pairs are ordered by wavelength
    */
@@ -144,6 +149,10 @@ public:
     return m_Interval;
   }
 
+  /** set index in m_Response vector to accelerate () operator **/
+  void SetPosGuessMin(const PrecisionType & lambda);
+
+
 protected:
   /** Constructor */
   SpectralResponse();
@@ -153,7 +162,7 @@ protected:
   virtual ~SpectralResponse()
   {
   }
- ;
+  ;
   /** PrintSelf method */
   //void PrintSelf(std::ostream& os, itk::Indent indent) const;
 
@@ -163,8 +172,9 @@ protected:
   /** Minimum value to consider that the spectral response is not null  */
   TPrecision m_SensitivityThreshold;
   IntervalType m_Interval;
+  unsigned long m_PosGuess;
   bool m_IntervalComputed;
-
+  bool m_UsePosGuess;
   void ComputeInterval();
 
 private:

Code/Simulation/otbSpectralResponse.txx

@@ -32,6 +32,8 @@ SpectralResponse<TPrecision, TValuePrecision>::SpectralResponse()
 {
   m_SensitivityThreshold = 0.01;
   m_IntervalComputed = false;
+  m_PosGuess = 0;
+  m_UsePosGuess=false;
 }
 
 template<class TPrecision, class TValuePrecision>
@@ -85,6 +87,34 @@ unsigned int SpectralResponse<TPrecision, TValuePrecision>::Size() const
   return m_Response.size();
 }
 
+
+
+template<class TPrecision, class TValuePrecision>
+void SpectralResponse<TPrecision, TValuePrecision>::SetPosGuessMin(const PrecisionType & lambda)
+{
+  m_PosGuess = 0;
+  if (m_Response.size() <= 1)
+    {
+    itkExceptionMacro(<<"ERROR spectral response need at least 2 value to perform interpolation.");
+    }
+
+  TPrecision lambdaMax = this->GetInterval().second;
+  if (lambda > lambdaMax) return;
+  typename VectorPairType::const_iterator it = m_Response.begin();
+
+  while (((*it).first < lambda))
+    {
+    m_PosGuess++;
+    ++it;
+    if (it == (m_Response.end())) return;
+    }
+
+  if (m_PosGuess > 0) m_PosGuess--;
+  return;
+}
+
+
+
 template<class TPrecision, class TValuePrecision>
 inline typename SpectralResponse<TPrecision, TValuePrecision>::ValuePrecisionType SpectralResponse<TPrecision,
     TValuePrecision>::operator()(const PrecisionType & lambda)
@@ -108,7 +138,12 @@ inline typename SpectralResponse<TPrecision, TValuePrecision>::ValuePrecisionTyp
   if (lambda < lambdaMin) return static_cast<TValuePrecision> (0.0);
   if (lambda > lambdaMax) return static_cast<TValuePrecision> (0.0);
 
-  typename VectorPairType::const_iterator it = m_Response.begin();
+  typename VectorPairType::const_iterator it;
+
+  if(m_UsePosGuess)
+    it= beg + m_PosGuess;
+  else
+    it= beg;
 
   TPrecision lambda1 = (*beg).first;
   TValuePrecision SR1 = static_cast<TValuePrecision> (0.0);