diff --git a/Code/Simulation/otbReduceSpectralResponse.txx b/Code/Simulation/otbReduceSpectralResponse.txx index 1dd1179a4aec4a7ba6760ca49c1cf378029ee4fb..7d4f1ff7480b7b5ff9b7af359b104ecf5a09912a 100644 --- a/Code/Simulation/otbReduceSpectralResponse.txx +++ b/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); + } } diff --git a/Code/Simulation/otbSpectralResponse.h b/Code/Simulation/otbSpectralResponse.h index 14834f7df2eca2f85109019c0f1f6113ba74569d..500a4f89e2c271311a88eef8ca7ffd6689f62c26 100644 --- a/Code/Simulation/otbSpectralResponse.h +++ b/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: diff --git a/Code/Simulation/otbSpectralResponse.txx b/Code/Simulation/otbSpectralResponse.txx index f34b9b01625bbfdaab604cbe4f4d909ed2566132..1b8c1d01dcbb78c86d76bf69b0ba1e41025ee29b 100644 --- a/Code/Simulation/otbSpectralResponse.txx +++ b/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);