ENH: normalize coordinates to avoid numerical problems, and solve system only when necessary

Code/Radiometry/otbSarParametricMapFunction.h

@@ -96,7 +96,6 @@ public:
   {
     m_IsInitialize = false;
     m_PointSet = val;
-    EvaluateParametricCoefficient();
     this->Modified();
   }
 

Code/Radiometry/otbSarParametricMapFunction.txx

@@ -82,10 +82,7 @@ SarParametricMapFunction<TInputImage, TCoordRep>
       ++pointId;
       }
     }
-  if (m_PointSet->GetNumberOfPoints() > 0)
-    {
-    EvaluateParametricCoefficient();
-    }
+
   this->Modified();
 }
 
@@ -94,8 +91,7 @@ void
 SarParametricMapFunction<TInputImage, TCoordRep>
 ::EvaluateParametricCoefficient()
 {
-  PointSetPointer pointSet;
-  pointSet = this->GetPointSet();
+  PointSetPointer pointSet = this->GetPointSet();
 
   PointType coef;
   PointType point;
@@ -118,6 +114,12 @@ SarParametricMapFunction<TInputImage, TCoordRep>
     }
   else
     {
+    InputImageType* inputImage = const_cast<InputImageType*>(this->GetInputImage());
+    //std::cout << inputImage << std::endl;
+    typename InputImageType::RegionType region = inputImage->GetLargestPossibleRegion();
+    typename InputImageType::IndexType  origin = region.GetIndex();
+    typename InputImageType::SizeType   size = region.GetSize();
+
     // Perform the plane least square estimation
     unsigned int nbRecords = pointSet->GetNumberOfPoints();
     unsigned int nbCoef = m_Coeff->GetNumberOfPoints();
@@ -141,8 +143,8 @@ SarParametricMapFunction<TInputImage, TCoordRep>
         PointType powerCoef;
         powerCoef.Fill(0);
         this->GetCoeff()->GetPoint(pointId, &powerCoef);
-        a(i, pointId) = vcl_pow(point[0], powerCoef[0]);
-        a(i, pointId) *= vcl_pow(point[1], powerCoef[1]);
+        a(i, pointId) = vcl_pow( (point[0] - origin[0]) / size[0], powerCoef[0]);
+        a(i, pointId) *= vcl_pow( (point[1] - origin[1]) / size[1], powerCoef[1]);
         //std::cout << "a(" << i << "," << pointId << ") = " << a(i, pointId) << std::endl;
         }
       }
@@ -188,13 +190,18 @@ SarParametricMapFunction<TInputImage, TCoordRep>
     return (itk::NumericTraits<RealType>::max());
     }
 
-  if (m_IsInitialize == false )
+  if (!m_IsInitialize)
     {
     itkExceptionMacro(<< "must estimate parameters before evaluating ");
     }
 
-  if(m_UsingClosestPointMethod == false )
+  if(!m_UsingClosestPointMethod)
     {
+    const InputImageType* inputImage = this->GetInputImage();
+    typename InputImageType::RegionType region = inputImage->GetLargestPossibleRegion();
+    typename InputImageType::IndexType  origin = region.GetIndex();
+    typename InputImageType::SizeType   size = region.GetSize();
+
     for(unsigned int pointId = 0; pointId < m_Coeff->GetNumberOfPoints(); ++pointId)
       {
       PointType  powerCoef;
@@ -203,7 +210,11 @@ SarParametricMapFunction<TInputImage, TCoordRep>
       this->GetCoeff()->GetPoint(pointId, &powerCoef);
       this->GetCoeff()->GetPointData(pointId, &pointValue);
 
-      result += pointValue * vcl_pow(index[0],powerCoef[0]) * vcl_pow(index[1],powerCoef[1]);
+      PointType normalized_point;
+      normalized_point[0] = static_cast<typename PointType::ValueType>(index[0] - origin[0]) / size[0];
+      normalized_point[1] = static_cast<typename PointType::ValueType>(index[1] - origin[1]) / size[1];
+
+      result += pointValue * vcl_pow(normalized_point[0],powerCoef[0]) * vcl_pow(normalized_point[1],powerCoef[1]);
       }
     }
 

Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx

@@ -43,10 +43,13 @@ void
 SarRadiometricCalibrationToImageFilter<TInputImage, TOutputImage>
 ::BeforeThreadedGenerateData()
 {
+  // will SetInputImage on the function
+  Superclass::BeforeThreadedGenerateData();
+
   SarImageMetadataInterface::Pointer imageMetadataInterface = SarImageMetadataInterfaceFactory::CreateIMI(
       this->GetInput()->GetMetaDataDictionary());
 
-  FunctionPointer function = FunctionType::New();
+  FunctionPointer function = this->GetFunction();
 
   function->SetScale(imageMetadataInterface->GetRadiometricCalibrationScale());
   
@@ -59,30 +62,35 @@ SarRadiometricCalibrationToImageFilter<TInputImage, TOutputImage>
   noise = function->GetNoise();
   noise->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationNoise());
   noise->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationNoisePolynomialDegree());
-  function->SetNoise(noise);
+  noise->EvaluateParametricCoefficient();
+  //function->SetNoise(noise);
   
   antennaPatternNewGain = function->GetAntennaPatternNewGain();
   antennaPatternNewGain->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternNewGain());
   antennaPatternNewGain->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternNewGainPolynomialDegree());
-  function->SetAntennaPatternNewGain(antennaPatternNewGain);
+  antennaPatternNewGain->EvaluateParametricCoefficient();
+  //function->SetAntennaPatternNewGain(antennaPatternNewGain);
   
   antennaPatternOldGain = function->GetAntennaPatternOldGain();
   antennaPatternOldGain->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternOldGain());
   antennaPatternOldGain->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternOldGainPolynomialDegree());
-  function->SetAntennaPatternNewGain(antennaPatternOldGain);
+  antennaPatternOldGain->EvaluateParametricCoefficient();
+  //function->SetAntennaPatternNewGain(antennaPatternOldGain);
   
   incidenceAngle = function->GetIncidenceAngle();
   incidenceAngle->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationIncidenceAngle());
   incidenceAngle->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationIncidenceAnglePolynomialDegree());
-  function->SetAntennaPatternNewGain(incidenceAngle);
+  incidenceAngle->EvaluateParametricCoefficient();
+  //function->SetAntennaPatternNewGain(incidenceAngle);
   
   rangeSpreadLoss = function->GetRangeSpreadLoss();
   rangeSpreadLoss->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationRangeSpreadLoss());
   rangeSpreadLoss->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationRangeSpreadLossPolynomialDegree());
-  function->SetAntennaPatternNewGain(rangeSpreadLoss);
+  rangeSpreadLoss->EvaluateParametricCoefficient();
+  //function->SetAntennaPatternNewGain(rangeSpreadLoss);
   
-  this->SetFunction(function);
-  Superclass::BeforeThreadedGenerateData();
+  //this->SetFunction(function);
+
 }
 
 

Testing/Code/Radiometry/otbSarParametricMapFunctionToImageFilter.cxx

@@ -80,9 +80,9 @@ int otbSarParametricMapFunctionToImageFilter(int argc, char * argv[])
   polynomalSize[1] = 0;
   function->SetPolynomalSize(polynomalSize);
 
+  function->EvaluateParametricCoefficient();
   std::cout << function << std::endl;     
 
-
   filter->SetFunction(function);
 
   writer->SetInput(filter->GetOutput());