From 8c12ae4e20676852d80363600404e727cb35571c Mon Sep 17 00:00:00 2001
From: Julien Malik <julien.malik@c-s.fr>
Date: Fri, 29 Oct 2010 23:58:35 +0200
Subject: [PATCH] ENH: normalize coordinates to avoid numerical problems, and
solve system only when necessary
---
Code/Radiometry/otbSarParametricMapFunction.h | 1 -
.../otbSarParametricMapFunction.txx | 33 ++++++++++++-------
...SarRadiometricCalibrationToImageFilter.txx | 24 +++++++++-----
...bSarParametricMapFunctionToImageFilter.cxx | 2 +-
4 files changed, 39 insertions(+), 21 deletions(-)
diff --git a/Code/Radiometry/otbSarParametricMapFunction.h b/Code/Radiometry/otbSarParametricMapFunction.h
index 445e37a8cb..4e06ce074f 100644
--- a/Code/Radiometry/otbSarParametricMapFunction.h
+++ b/Code/Radiometry/otbSarParametricMapFunction.h
@@ -96,7 +96,6 @@ public:
{
m_IsInitialize = false;
m_PointSet = val;
- EvaluateParametricCoefficient();
this->Modified();
}
diff --git a/Code/Radiometry/otbSarParametricMapFunction.txx b/Code/Radiometry/otbSarParametricMapFunction.txx
index c39ca20444..1f40fe1d66 100644
--- a/Code/Radiometry/otbSarParametricMapFunction.txx
+++ b/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]);
}
}
diff --git a/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx b/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx
index 6dc07b9363..bda03a406a 100644
--- a/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx
+++ b/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);
+
}
diff --git a/Testing/Code/Radiometry/otbSarParametricMapFunctionToImageFilter.cxx b/Testing/Code/Radiometry/otbSarParametricMapFunctionToImageFilter.cxx
index 5d841d2912..14840e4c24 100644
--- a/Testing/Code/Radiometry/otbSarParametricMapFunctionToImageFilter.cxx
+++ b/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());
--
GitLab