ossimRpcSolver Class Reference

This currently only support Rational poilynomial B format. More...

#include <ossimRpcSolver.h>

Inheritance diagram for ossimRpcSolver:

Public Member Functions
	ossimRpcSolver (bool useElevation=false, bool useHeightAboveMSLFlag=false)
	The use elvation flag will deterimne if we force the height t be 0. More...
virtual	~ossimRpcSolver ()
void	solveCoefficients (const ossimDrect &imageBounds, ossimProjection *imageProj, ossim_uint32 xSamples=8, ossim_uint32 ySamples=8)
	This will convert any projector to an RPC model. More...
void	solveCoefficients (const ossimDrect &imageBounds, ossimImageGeometry *geom, ossim_uint32 xSamples=8, ossim_uint32 ySamples=8)
bool	solve (const ossimDrect &aoiBounds, ossimImageGeometry *geom, const double &pixel_tolerance=0.5)
	Similar to the other solve methods except that the final grid size is established iteratively so that the error at the midpoints between grid nodes falls below tolerance. More...
bool	solve (const ossimFilename &imageFilename, const double &pixel_tolerance=0.5)
	Performs iterative solve using the other solve method, but uses an image filename to initialize, and computes RPC over entire image rect. More...
void	solveCoefficients (const std::vector< ossimDpt > &imagePoints, const std::vector< ossimGpt > &groundControlPoints)
	takes associated image points and ground points and solves the coefficents for the rational polynomial for line and sample calculations from world points. More...
const ossimRefPtr< ossimRpcModel >	getRpcModel () const
	Fetches the solved-for RPC model. More...
double	getRmsError () const
double	getMaxError () const
ossimRefPtr< ossimNitfRegisteredTag >	getNitfRpcBTag () const
void	setValidImageRect (const ossimIrect &imageRect)
	Sets the image rect over which to compute the RPCs. More...
Public Member Functions inherited from ossimReferenced
	ossimReferenced ()
	ossimReferenced (const ossimReferenced &)
ossimReferenced &	operator= (const ossimReferenced &)
void	ref () const
	increment the reference count by one, indicating that this object has another pointer which is referencing it. More...
void	unref () const
	decrement the reference count by one, indicating that a pointer to this object is referencing it. More...
void	unref_nodelete () const
	decrement the reference count by one, indicating that a pointer to this object is referencing it. More...
int	referenceCount () const

Protected Member Functions
virtual void	solveInitialCoefficients (NEWMAT::ColumnVector &coeff, const std::vector< double > &f, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z) const
virtual void	solveCoefficients (NEWMAT::ColumnVector &coeff, const std::vector< double > &f, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z) const
double	eval (const std::vector< double > &coeff, const double &x, const double &y, const double &z) const
void	evalPoint (const ossimGpt &gpt, ossimDpt &ipt) const
NEWMAT::Matrix	invert (const NEWMAT::Matrix &m) const
	Inverts using the SVD method. More...
void	setupSystemOfEquations (NEWMAT::Matrix &equations, const NEWMAT::ColumnVector &f, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z) const
void	setupWeightMatrix (NEWMAT::DiagonalMatrix &result, const NEWMAT::ColumnVector &coefficients, const NEWMAT::ColumnVector &f, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z) const
Protected Member Functions inherited from ossimReferenced
virtual	~ossimReferenced ()

Protected Attributes
bool	theUseElevationFlag
bool	theHeightAboveMSLFlag
ossim_float64	theMeanResidual
ossim_float64	theMaxResidual
ossimRefPtr< ossimImageGeometry >	theRefGeom
ossimRefPtr< ossimRpcModel >	theRpcModel

Detailed Description

This currently only support Rational poilynomial B format.

This can be found in the NITF registered commercial tag document.

Note

x=longitude, y=latitude, z=height

Format is:
 coeff[ 0]       + coeff[ 1]*x     + coeff[ 2]*y     + coeff[ 3]*z     +
 coeff[ 4]*x*y   + coeff[ 5]*x*z   + coeff[ 6]*y*z   + coeff[ 7]*x*x   +
 coeff[ 8]*y*y   + coeff[ 9]*z*z   + coeff[10]*x*y*z + coeff[11]*x*x*x +
 coeff[12]*x*y*y + coeff[13]*x*z*z + coeff[14]*x*x*y + coeff[15]*y*y*y +
 coeff[16]*y*z*z + coeff[17]*x*x*z + coeff[18]*y*y*z + coeff[19]*z*z*z;
 where coeff is one of XNum, XDen, YNum, and YDen.  So there are 80
 coefficients all together.

Currently we use a linear least squares fit to solve the coefficients.
This is the simplest to implement.  We probably relly need a nonlinear
minimizer to fit the coefficients but I don't have time to experiment.
Levenberg Marquardt might be a solution to look into.

HOW TO USE:
   ossimRpcSolver solver;
   solver.solveCoefficients(rect,
                            *proj.get());

We can also call solve coefficients with a list of ground control points.
First is the list of image points followed by the ground points.
NOTE: Thes must be equal in size.
   solver.solveCoefficients(imagePoints,
                            groundPoints);

Once you call solveCoefficients you can create the projector:
   ossimRefPtr<ossimRpcProjection> rpc = solver.createRpcProjection();

Note that a sub-image bounding rect can be passed into the solve methods. This
constrains the solution fit to cover only that rectangle in the original image space,
but the image coordinates used are still based on the full image. If the intent is to
generate an RPC that will work for an image chip in that chip's local image coordinate
system (i.e., the UL corner of the chip is 0, 0), then you'll need to call
rpcModel->setImageOffset(chip_offset) on the output RPC model.

Definition at line 72 of file ossimRpcSolver.h.

Constructor & Destructor Documentation

ossimRpcSolver()

ossimRpcSolver::ossimRpcSolver	(	bool	useElevation = false,
		bool	useHeightAboveMSLFlag = false
	)

The use elvation flag will deterimne if we force the height t be 0.

If the elevation is enabled then we use the height field of the control points to determine the coefficients of the RPC00 polynomial. If its false then we will ignore the height by setting the height field to 0.0.

Note: even if the elevation is enabled all NAN heights are set to 0.0.

Definition at line 21 of file ossimRpcSolver.cpp.

:  theUseElevationFlag(useElevation),
   theHeightAboveMSLFlag(useHeightAboveMSLFlag),
   theMeanResidual(0),
   theMaxResidual(0)
{
}

~ossimRpcSolver()

virtual ossimRpcSolver::~ossimRpcSolver ( )

inlinevirtual

Definition at line 86 of file ossimRpcSolver.h.

{}

Member Function Documentation

eval()

double ossimRpcSolver::eval ( const std::vector< double > &  coeff, const double &  x, const double &  y, const double &  z  ) const

protected

evalPoint()

void ossimRpcSolver::evalPoint ( const ossimGpt &  gpt, ossimDpt &  ipt  ) const

protected

Definition at line 677 of file ossimRpcSolver.cpp.

References ossimDpt::makeNan(), theRpcModel, and ossimRpcModel::worldToLineSample().

Referenced by solveCoefficients().

{
   if (!theRpcModel)
   {
      ipt.makeNan();
      return;
   }

   theRpcModel->worldToLineSample(gpt, ipt);
}

getMaxError()

double ossimRpcSolver::getMaxError

(

)

const

Definition at line 426 of file ossimRpcSolver.cpp.

References theMaxResidual.

{
   return theMaxResidual;
}

getNitfRpcBTag()

ossimRefPtr< ossimNitfRegisteredTag > ossimRpcSolver::getNitfRpcBTag

(

)

const

Returns

ossimRefPtr<ossimNitfRegisteredTag>

Note

one of the solve methods should have been called prior to this.

Definition at line 689 of file ossimRpcSolver.cpp.

References ossimNitfRpcBase::setRpcModelParams(), and theRpcModel.

Referenced by ossimNitfWriterBase::addRpcbTag().

{
   ossimNitfRpcBTag* rpcbTag = new ossimNitfRpcBTag();
   rpcbTag->setRpcModelParams(theRpcModel);
   
   // Return it as an ossimRefPtr<ossimNitfRegisteredTag>...
   ossimRefPtr<ossimNitfRegisteredTag> tag = rpcbTag;
   
   return tag;
}

getRmsError()

double ossimRpcSolver::getRmsError

(

)

const

Definition at line 421 of file ossimRpcSolver.cpp.

References theMeanResidual.

Referenced by ossimRpcProjection::optimizeFit().

{
   return theMeanResidual;
}

getRpcModel()

const ossimRefPtr<ossimRpcModel> ossimRpcSolver::getRpcModel ( ) const

inline

Fetches the solved-for RPC model.

See note in header above on setting the image offset if this model will be applied to a sub-image chip.

Definition at line 138 of file ossimRpcSolver.h.

Referenced by ossimDemTool::doASP(), ossimSubImageTool::execute(), and ossimRpcProjection::optimizeFit().

{ return theRpcModel; }

invert()

NEWMAT::Matrix ossimRpcSolver::invert ( const NEWMAT::Matrix &  m ) const

protected

Inverts using the SVD method.

Definition at line 544 of file ossimRpcSolver.cpp.

References FLT_EPSILON, and SVD().

Referenced by solveCoefficients(), and solveInitialCoefficients().

{
   ossim_uint32 idx = 0;
   NEWMAT::DiagonalMatrix d;
   NEWMAT::Matrix u;
   NEWMAT::Matrix v;

   // decompose m.t*m which is stored in Temp into the singular values and vectors.
   //
   NEWMAT::SVD(m, d, u, v, true, true);

   // invert the diagonal
   // this is just doing the reciprical fo all diagonal components and store back int
   // d.  ths compute d inverse.
   //
   for(idx=0; idx < (ossim_uint32)d.Ncols(); ++idx)
   {
      if(d[idx] > FLT_EPSILON)
      {
         d[idx] = 1.0/d[idx];
      }
      else
      {
         d[idx] = 0.0;
      }
   }

   //compute inverse of decomposed m;
   return v*d*u.t();
}

setupSystemOfEquations()

void ossimRpcSolver::setupSystemOfEquations ( NEWMAT::Matrix &  equations, const NEWMAT::ColumnVector &  f, const std::vector< double > &  x, const std::vector< double > &  y, const std::vector< double > &  z  ) const

protected

Definition at line 576 of file ossimRpcSolver.cpp.

References x, and y.

Referenced by solveCoefficients(), and solveInitialCoefficients().

{
   ossim_uint32 idx;
   equations.ReSize(f.Nrows(), 39);
   
   for(idx = 0; idx < (ossim_uint32)f.Nrows();++idx)
   {
      equations[idx][0]  = 1;
      equations[idx][1]  = x[idx];
      equations[idx][2]  = y[idx];
      equations[idx][3]  = z[idx];
      equations[idx][4]  = x[idx]*y[idx];
      equations[idx][5]  = x[idx]*z[idx];
      equations[idx][6]  = y[idx]*z[idx];
      equations[idx][7]  = x[idx]*x[idx];
      equations[idx][8]  = y[idx]*y[idx];
      equations[idx][9]  = z[idx]*z[idx];
      equations[idx][10] = x[idx]*y[idx]*z[idx];
      equations[idx][11] = x[idx]*x[idx]*x[idx];
      equations[idx][12] = x[idx]*y[idx]*y[idx];
      equations[idx][13] = x[idx]*z[idx]*z[idx];
      equations[idx][14] = x[idx]*x[idx]*y[idx];
      equations[idx][15] = y[idx]*y[idx]*y[idx];
      equations[idx][16] = y[idx]*z[idx]*z[idx];
      equations[idx][17] = x[idx]*x[idx]*z[idx];
      equations[idx][18] = y[idx]*y[idx]*z[idx];
      equations[idx][19] = z[idx]*z[idx]*z[idx];
      equations[idx][20] = -f[idx]*x[idx];
      equations[idx][21] = -f[idx]*y[idx];
      equations[idx][22] = -f[idx]*z[idx];
      equations[idx][23] = -f[idx]*x[idx]*y[idx];
      equations[idx][24] = -f[idx]*x[idx]*z[idx];
      equations[idx][25] = -f[idx]*y[idx]*z[idx];
      equations[idx][26] = -f[idx]*x[idx]*x[idx];
      equations[idx][27] = -f[idx]*y[idx]*y[idx];
      equations[idx][28] = -f[idx]*z[idx]*z[idx];
      equations[idx][29] = -f[idx]*x[idx]*y[idx]*z[idx];
      equations[idx][30] = -f[idx]*x[idx]*x[idx]*x[idx];
      equations[idx][31] = -f[idx]*x[idx]*y[idx]*y[idx];
      equations[idx][32] = -f[idx]*x[idx]*z[idx]*z[idx];
      equations[idx][33] = -f[idx]*x[idx]*x[idx]*y[idx];
      equations[idx][34] = -f[idx]*y[idx]*y[idx]*y[idx];
      equations[idx][35] = -f[idx]*y[idx]*z[idx]*z[idx];
      equations[idx][36] = -f[idx]*x[idx]*x[idx]*z[idx];
      equations[idx][37] = -f[idx]*y[idx]*y[idx]*z[idx];
      equations[idx][38] = -f[idx]*z[idx]*z[idx]*z[idx];
   }
}

setupWeightMatrix()

void ossimRpcSolver::setupWeightMatrix ( NEWMAT::DiagonalMatrix &  result, const NEWMAT::ColumnVector &  coefficients, const NEWMAT::ColumnVector &  f, const std::vector< double > &  x, const std::vector< double > &  y, const std::vector< double > &  z  ) const

protected

Definition at line 629 of file ossimRpcSolver.cpp.

References FLT_EPSILON, x, and y.

Referenced by solveCoefficients().

{
   result.ReSize(f.Nrows());
   ossim_uint32 idx = 0;
   ossim_uint32 idx2 = 0;
   NEWMAT::RowVector row(coefficients.Nrows());
   
    for(idx = 0; idx < (ossim_uint32)f.Nrows(); ++idx)
    {
       row[0]  = 1;
       row[1]  = x[idx];
       row[2]  = y[idx];
       row[3]  = z[idx];
       row[4]  = x[idx]*y[idx];
       row[5]  = x[idx]*z[idx];
       row[6]  = y[idx]*z[idx];
       row[7]  = x[idx]*x[idx];
       row[8]  = y[idx]*y[idx];
       row[9]  = z[idx]*z[idx];
       row[10] = x[idx]*y[idx]*z[idx];
       row[11] = x[idx]*x[idx]*x[idx];
       row[12] = x[idx]*y[idx]*y[idx];
       row[13] = x[idx]*z[idx]*z[idx];
       row[14] = x[idx]*x[idx]*y[idx];
       row[15] = y[idx]*y[idx]*y[idx];
       row[16] = y[idx]*z[idx]*z[idx];
       row[17] = x[idx]*x[idx]*z[idx];
       row[18] = y[idx]*y[idx]*z[idx];
       row[19] = z[idx]*z[idx]*z[idx];

      result[idx] = 0.0;
      for(idx2 = 0; idx2 < (ossim_uint32)row.Ncols(); ++idx2)
      {
         result[idx] += row[idx2]*coefficients[idx2];
      }

      if(result[idx] > FLT_EPSILON)
      {
         result[idx] = 1.0/result[idx];
      }
    }
}

setValidImageRect()

void ossimRpcSolver::setValidImageRect

(

const ossimIrect &

imageRect

)

Sets the image rect over which to compute the RPCs.

The Resulting RPC will only be valid over that range of image space.

solve() [1/2]

bool ossimRpcSolver::solve	(	const ossimDrect &	aoiBounds,
		ossimImageGeometry *	geom,
		const double &	pixel_tolerance = 0.5
	)

Similar to the other solve methods except that the final grid size is established iteratively so that the error at the midpoints between grid nodes falls below tolerance.

The RPC is computed for the specified image bounds range only, not the full image. The expectation here (when the imageBounds is less than the full valid image rect) is to generate an RPC to accompany a subimage that will be written to disk.

Parameters

imageBounds	The AOI in image space for the RPC computation.
geom	Represents the geometry of the input image
pixel_tolerance	Maximum error in pixels (typically fraction of a pixel) to achieve.

Returns

true if solution converged below pixel tolerance.

Definition at line 290 of file ossimRpcSolver.cpp.

References ossimDrect::height(), theRefGeom, ossimDrect::ul(), and ossimDrect::width().

Referenced by ossimDemTool::doASP(), ossimSubImageTool::execute(), and solve().

{
   static const char* MODULE = "ossimRpcSolver::solve()  ";

   if (!geom)
      return false;

   theRefGeom = geom;
   ossimDpt ul = imageBounds.ul();
   ossim_float64 w = imageBounds.width();
   ossim_float64 h = imageBounds.height();
   ossimDpt ipt, irpc;
   ossimGpt gpt;

   // Start at the minimum grid size:
   ossim_uint32 xSamples = STARTING_GRID_SIZE;
   ossim_uint32 ySamples = STARTING_GRID_SIZE;

   // Loop until error is below threshold:
   bool converged = false;
   while (!converged)
   {
      double residual = 0;
      double sumResiduals = 0;
      int numResiduals = 0;
      theMaxResidual = 0;

      converged = true; // hope for the best and get proved otherwise below
      solveCoefficients(imageBounds, geom, xSamples, ySamples);

      // Sample along x and y directions to accumulate errors:
      double deltaX = w/(xSamples-1);
      double deltaY = h/(ySamples-1);

      // Sample the midpoints between image grid used to compute RPC:
      for (ossim_uint32 y=0; y<ySamples-1; ++y)
      {
         ipt.y = deltaY*((double)y + 0.5) + ul.y;
         for (ossim_uint32 x=0; x<xSamples-1; ++x)
         {
            // Forward projection using input model:
            ipt.x = deltaX*((double)x + 0.5) + ul.x;
            if (theUseElevationFlag)
               geom->localToWorld(ipt, gpt);
            else
               geom->localToWorld(ipt, 0, gpt);
            if(theHeightAboveMSLFlag)
            {
               double h = ossimElevManager::instance()->getHeightAboveMSL(gpt);
               if(ossim::isnan(h) == false)
                  gpt.height(h);
            }

            // Reverse projection using RPC:
            evalPoint(gpt, irpc);

            // Compute residual and accumulate:
            residual = (ipt-irpc).length();
            if (residual > theMaxResidual)
               theMaxResidual = residual;
            sumResiduals += residual;
            ++numResiduals;
         }
      }

      theMeanResidual = sumResiduals/numResiduals;
      if (theMaxResidual > tolerance)
         converged = false;

#if 1
      { //### DEBUG BLOCK ###
         ossimNotify(ossimNotifyLevel_INFO)<<MODULE
               <<"\n   sampling grid size: ("<<xSamples<<", "<<ySamples<<")"
               <<"\n        mean residual: "<<theMeanResidual
               <<"\n         max residual: "<<theMaxResidual
               <<"\n            converged: "<<ossimString::toString(converged)<<endl;
      }
#endif

      if ((xSamples >= ENDING_GRID_SIZE) && (ySamples >= ENDING_GRID_SIZE))
         break;

      if (!converged)
      {
         xSamples *= 2;
         if (xSamples > ENDING_GRID_SIZE)
            xSamples = ENDING_GRID_SIZE;
         ySamples *= 2;
         if (ySamples > ENDING_GRID_SIZE)
            ySamples = ENDING_GRID_SIZE;
      }
   }

   // Was testing the max residual. But use the mean for final test:
   if (theMeanResidual <= tolerance)
         converged = true;

   if (!converged)
   {
      ossimNotify(ossimNotifyLevel_WARN)
            << "WARNING: Unable to converge on desired error tolerance ("<<tolerance<<" p).\n"
            <<"    RMS residual error: " << theMeanResidual << "\n"
            <<"    Max residual error: " << theMaxResidual<<std::endl;
   }
   else if (theMaxResidual > tolerance)
   {
      ossimNotify(ossimNotifyLevel_WARN)
            << "WARNING: While the RPC solution did converge, at least one residual ("
            <<theMaxResidual<<") is larger than the desired error tolerance ("<<tolerance<<" p)."
            << std::endl;
   }
   return converged;
}

solve() [2/2]

bool ossimRpcSolver::solve	(	const ossimFilename &	imageFilename,
		const double &	pixel_tolerance = 0.5
	)

Performs iterative solve using the other solve method, but uses an image filename to initialize, and computes RPC over entire image rect.

Definition at line 406 of file ossimRpcSolver.cpp.

References ossimRefPtr< T >::get(), ossimImageHandler::getBoundingRect(), ossimImageHandler::getImageGeometry(), ossimImageHandlerRegistry::instance(), ossimImageHandlerRegistry::open(), solve(), and ossimRefPtr< T >::valid().

{
   // Establish input geometry:
   ossimRefPtr<ossimImageHandler> h = ossimImageHandlerRegistry::instance()->open(imageFilename);
   if(!h.valid())
      return false;

   ossimRefPtr<ossimImageGeometry> geom = h->getImageGeometry();
   ossimDrect imageRect (h->getBoundingRect());

   return solve(imageRect, geom.get(), pixel_tolerance);
}

solveCoefficients() [1/4]

void ossimRpcSolver::solveCoefficients	(	const ossimDrect &	imageBounds,
		ossimProjection *	imageProj,
		ossim_uint32	xSamples = 8,
		ossim_uint32	ySamples = 8
	)

This will convert any projector to an RPC model.

Definition at line 29 of file ossimRpcSolver.cpp.

References ossimRefPtr< T >::get(), and ossimImageGeometry::setProjection().

Referenced by ossimNitfWriterBase::addRpcbTag(), ossimRpcProjection::optimizeFit(), and solveCoefficients().

{
   ossimRefPtr<ossimImageGeometry> geom = new ossimImageGeometry();
   geom->setProjection(proj);
   solveCoefficients(imageBounds, geom.get(), xSamples, ySamples);
}

solveCoefficients() [2/4]

void ossimRpcSolver::solveCoefficients	(	const ossimDrect &	imageBounds,
		ossimImageGeometry *	geom,
		ossim_uint32	xSamples = 8,
		ossim_uint32	ySamples = 8
	)

Definition at line 39 of file ossimRpcSolver.cpp.

References ossimImageGeometry::getProjection(), theRefGeom, x, and y.

{
   if (!geom || !(geom->getProjection()))
      return;
   theRefGeom = geom;

   std::vector<ossimGpt> groundPoints;
   std::vector<ossimDpt> imagePoints;
   ossim_uint32 x,y;
   ossimGpt gpt;
   ossimGpt defaultGround;
   if (ySamples <= 1)
      ySamples = STARTING_GRID_SIZE;
   if (xSamples <= 1)
      xSamples = STARTING_GRID_SIZE;
   srand(time(0));
   double Dx = imageBounds.width()/(xSamples-1);
   double Dy = imageBounds.height()/(ySamples-1);
   ossimDpt dpt;
   for(y = 0; y < ySamples; ++y)
   {
      dpt.y = y*Dy + imageBounds.ul().y;
      for(x = 0; x < xSamples; ++x)
      {
         dpt.x = x*Dx + imageBounds.ul().x;
         if (theUseElevationFlag)
            geom->localToWorld(dpt, gpt);
         else
            geom->localToWorld(dpt, 0, gpt);

         if (gpt.isLatLonNan())
            continue;

         if(gpt.isHgtNan())
            gpt.height(0.0);

         gpt.changeDatum(defaultGround.datum());
         if(theHeightAboveMSLFlag)
         {
            double h = ossimElevManager::instance()->getHeightAboveMSL(gpt);
            if(ossim::isnan(h) == false)
               gpt.height(h);
         }

         imagePoints.push_back(dpt);
         groundPoints.push_back(gpt);
      }
   }
   solveCoefficients(imagePoints, groundPoints);
}

solveCoefficients() [3/4]

void ossimRpcSolver::solveCoefficients	(	const std::vector< ossimDpt > &	imagePoints,
		const std::vector< ossimGpt > &	groundControlPoints
	)

takes associated image points and ground points and solves the coefficents for the rational polynomial for line and sample calculations from world points.

Note: All data will be normalized between -1 and 1 for numerical robustness.

Definition at line 93 of file ossimRpcSolver.cpp.

References ossimRpcModel::B, DBL_EPSILON, evalPoint(), FLT_EPSILON, ossimGpt::height(), ossimDrect::height(), ossim::isnan(), ossimDrect::midPoint(), solveCoefficients(), ossimRpcModel::theHgtOffset, ossimRpcModel::theHgtScale, ossimRpcModel::theLatOffset, ossimRpcModel::theLatScale, ossimRpcModel::theLineDenCoef, ossimRpcModel::theLineNumCoef, ossimRpcModel::theLineOffset, ossimRpcModel::theLineScale, ossimRpcModel::theLonOffset, ossimRpcModel::theLonScale, theMaxResidual, theMeanResidual, ossimRpcModel::thePolyType, theRpcModel, ossimRpcModel::theSampDenCoef, ossimRpcModel::theSampNumCoef, ossimRpcModel::theSampOffset, ossimRpcModel::theSampScale, theUseElevationFlag, ossimDrect::width(), ossimDpt::x, x, ossimDpt::y, and y.

{
   if((imagePoints.size() != groundControlPoints.size()))
      return;

   // we will first create f which holds the result of f(x,y,z).
   // This basically holds the cooresponding image point for each
   // ground control point.  One for x and a second array for y
   int numPoints = imagePoints.size();
   std::vector<double> fx, fy;

   //  Holds the x, y, z vectors
   //
   std::vector<double> x;
   std::vector<double> y;
   std::vector<double> z;
   ossim_uint32 c = 0;
   fx.resize(imagePoints.size());
   fy.resize(imagePoints.size());
   x.resize(imagePoints.size());
   y.resize(imagePoints.size());
   z.resize(imagePoints.size());

   // compute the image bounds for the given image points
   ossimDrect rect(imagePoints);
   ossimDpt centerImagePoint  = rect.midPoint();

   // get the width and height that will be used in data normalization
   ossim_float64 w = rect.width();
   ossim_float64 h = rect.height();

   double latSum=0.0;
   double lonSum=0.0;
   double heightSum=0.0;

   // find the center ground  Use elevation only if its enabled
   for(c = 0; c < groundControlPoints.size();++c)
   {
      if(ossim::isnan(groundControlPoints[c].latd()) == false)
      {
         latSum += groundControlPoints[c].latd();
      }
      if(ossim::isnan(groundControlPoints[c].lond()) == false)
      {
         lonSum += groundControlPoints[c].lond();
      }
      if(!groundControlPoints[c].isHgtNan())
      {
         if(theUseElevationFlag)
         {
            heightSum += groundControlPoints[c].height();
         }
      }
   }

   // set the center ground for the offset
   //
   ossimGpt centerGround(latSum/groundControlPoints.size(),
                         lonSum/groundControlPoints.size(),
                         heightSum/groundControlPoints.size());

   // set up ground scales and deltas for normalization
   ossim_float64 deltaLat       = 0.0;
   ossim_float64 deltaLon       = 0.0;
   ossim_float64 deltaHeight    = 0.0;
   ossim_float64 maxDeltaLat    = 0.0;
   ossim_float64 maxDeltaLon    = 0.0;
   ossim_float64 maxDeltaHeight = 0.0;
   ossim_float64 heightTest       = 0.0;
   for(c = 0; c < groundControlPoints.size(); ++c)
   {
      deltaLat = (groundControlPoints[c].latd()-centerGround.latd());
      deltaLon = (groundControlPoints[c].lond()-centerGround.lond());
      if(!groundControlPoints[c].isHgtNan())
      {
         if(theUseElevationFlag)
         {
            deltaHeight = groundControlPoints[c].height() - centerGround.height();
            heightTest  = groundControlPoints[c].height();
         }
         else
         {
            deltaHeight = 0.0;
            heightTest  = 0.0;
         }
      }
      else
      {
         deltaHeight = 0.0;
      }
      fx[c] = (imagePoints[c].x - centerImagePoint.x)/(w/2.0);
      fy[c] = (imagePoints[c].y - centerImagePoint.y)/(h/2.0);
      
      x[c] = deltaLon;
      y[c] = deltaLat;
      z[c] = deltaHeight;

      if(fabs(deltaLat) > maxDeltaLat)
         maxDeltaLat = fabs(deltaLat);
      if(fabs(deltaLon) > maxDeltaLon)
         maxDeltaLon = fabs(deltaLon);
      if(fabs(heightTest) > maxDeltaHeight)
         maxDeltaHeight = fabs(heightTest);
   }

   bool elevationEnabled = theUseElevationFlag;

   // if max delta is less than a degree set it to 1 degree.
//   if(maxDeltaLat < 1.0)
//      maxDeltaLat = 1.0;
//   if(maxDeltaLon < 1.0)
//      maxDeltaLon = 1.0;

   if(maxDeltaHeight < FLT_EPSILON)
      elevationEnabled = false;
   if(maxDeltaHeight < 1.0)
      maxDeltaHeight = 1.0;

   // set the height scale to something pretty large
   if(!elevationEnabled)
   {
      maxDeltaHeight = 1.0/DBL_EPSILON;
      maxDeltaHeight = 10000;
      centerGround.height(0.0);
   }
   // normalize the ground points
   for(c = 0; c < groundControlPoints.size(); ++c)
   {
      x[c] /= maxDeltaLon;
      y[c] /= maxDeltaLat;
      z[c] /= maxDeltaHeight;
   }

   theRpcModel = new ossimRpcModel;
   theRpcModel->thePolyType = ossimRpcModel::B;
   theRpcModel->theLineOffset = centerImagePoint.y;
   theRpcModel->theSampOffset = centerImagePoint.x;
   theRpcModel->theLineScale  = h/2.0;
   theRpcModel->theSampScale  = w/2.0;

   theRpcModel->theLatScale = maxDeltaLat;
   theRpcModel->theLonScale = maxDeltaLon;
   theRpcModel->theHgtScale = maxDeltaHeight;

   theRpcModel->theLatOffset = centerGround.lat;
   theRpcModel->theLonOffset = centerGround.lon;
   theRpcModel->theHgtOffset = centerGround.hgt;

   if(ossim::isnan(theRpcModel->theHgtOffset))
      theRpcModel->theHgtOffset = 0.0;

   NEWMAT::ColumnVector coeffxVec;
   NEWMAT::ColumnVector coeffyVec;

   // perform a least squares fit for sample values found in f
   // given the world values with variables x, y, z
   solveCoefficients(coeffxVec, fx, x, y, z);

   // perform a least squares fit for line values found in f
   // given the world values with variables x, y, z
   solveCoefficients(coeffyVec, fy, x, y, z);

   // there are 20 numerator coefficients and 19 denominator coefficients.
   // I believe that the very first one for the denominator coefficients is fixed at 1.
   theRpcModel->theLineNumCoef[0] = coeffyVec[0];
   theRpcModel->theLineDenCoef[0] = 1.0;
   theRpcModel->theSampNumCoef[0] = coeffxVec[0];
   theRpcModel->theSampDenCoef[0] = 1.0;
   for (int i=1; i<20; i++)
   {
      theRpcModel->theLineNumCoef[i] = coeffyVec[i];
      theRpcModel->theLineDenCoef[i] = coeffyVec[i+19];
      theRpcModel->theSampNumCoef[i] = coeffxVec[i];
      theRpcModel->theSampDenCoef[i] = coeffxVec[i+19];
   }

   // now lets compute the RMSE for the given control points by feeding it
   // back through the modeled RPC
   ossim_float64  sumSquareError = 0.0;
   ossim_uint32 idx = 0;

   theMaxResidual = 0;
   for (idx = 0; idx<imagePoints.size(); idx++)
   {
      ossimDpt evalPt;
      evalPoint(groundControlPoints[idx], evalPt);
      ossim_float64 len = (evalPt - imagePoints[idx]).length();
      if (len > theMaxResidual)
         theMaxResidual = len;
      sumSquareError += (len*len);
   }

   // set the error
   theMeanResidual = sqrt(sumSquareError/imagePoints.size());
}

solveCoefficients() [4/4]

void ossimRpcSolver::solveCoefficients ( NEWMAT::ColumnVector &  coeff, const std::vector< double > &  f, const std::vector< double > &  x, const std::vector< double > &  y, const std::vector< double > &  z  ) const

protectedvirtual

Definition at line 449 of file ossimRpcSolver.cpp.

References FLT_EPSILON, invert(), setupSystemOfEquations(), setupWeightMatrix(), x, and y.

{
   // this is an iterative  linear least square fit.  We really pobably need
   // a nonlinear fit instead
   //
   ossim_uint32 idx = 0;
   NEWMAT::Matrix m;

   NEWMAT::ColumnVector r((int)f.size());

   for(idx = 0; idx < f.size(); ++idx)
   {
      r[idx] = f[idx];
   }

   NEWMAT::ColumnVector tempCoeff;
   NEWMAT::DiagonalMatrix weights((int)f.size());
   NEWMAT::ColumnVector denominator(20);

   // initialize the weight matrix to the identity
   //
   for(idx = 0; idx < f.size(); ++idx)
   {
      weights[idx] = 1.0;
   }

   double residualValue = 1.0/FLT_EPSILON;
   ossim_uint32 iterations = 0;
   NEWMAT::Matrix w2;
   do
   {
      w2 = weights*weights;

#if 0
      { //### DEBUG ###
         cout<<"\nw2 = \n"<<w2<<endl;
         cout<<"\nr = "<<r<<endl;
      }
#endif

      // sets up the matrix to hold the system of equations
      setupSystemOfEquations(m, r, x, y, z);

      // solve the least squares solution.  Note: the invert is used
      // to do a Singular Value Decomposition for the inverse since the
      // matrix is more than likely singular.  Slower but more robust
#if 0
      { //### DEBUG ###
         NEWMAT::Matrix mt = m.t();
         cout<<"\nm = "<<m<<endl;
         cout<<"\nmt = "<<mt<<endl;

         NEWMAT::Matrix mtw2 = m.t()*w2;
         cout<<"\nmtw2 = "<<mtw2<<endl;
         NEWMAT::Matrix mtw2m = mtw2*m;
         cout<<"\nmtw2m = "<<mtw2m<<endl;
         NEWMAT::Matrix mtw2r = mtw2*r;
         cout<<"\nmtw2r = "<<mtw2r<<endl;

         NEWMAT::Matrix mtw2m_inv = invert(mtw2m);
         cout<<"\nmtw2m_inv = "<<mtw2m_inv<<endl;
         tempCoeff = mtw2m_inv * mtw2r;
         cout<<"\ntempCoeff = "<<tempCoeff<<endl;
      }
#else
      tempCoeff = invert(m.t()*w2*m)*m.t()*w2*r;
#endif

      // set up the weight matrix by using the denominator
      for(idx = 0; idx < 19; ++idx)
      {
         denominator[idx+1] = tempCoeff[20+idx];
      }
      denominator[0] = 1.0;
      
      setupWeightMatrix(weights, denominator, r, x, y, z);

      // compute the residual
      NEWMAT::ColumnVector residual = m.t()*w2*(m*tempCoeff-r);

      // now get the innerproduct
      NEWMAT::Matrix tempRes = (residual.t()*residual);
      residualValue = sqrt(tempRes[0][0]);

      ++iterations;

   } while ((residualValue > FLT_EPSILON) && (iterations < 10));
   coeff = tempCoeff;

}

solveInitialCoefficients()

void ossimRpcSolver::solveInitialCoefficients ( NEWMAT::ColumnVector &  coeff, const std::vector< double > &  f, const std::vector< double > &  x, const std::vector< double > &  y, const std::vector< double > &  z  ) const

protectedvirtual

Definition at line 431 of file ossimRpcSolver.cpp.

References invert(), setupSystemOfEquations(), x, and y.

{
   ossim_uint32 idx = 0;
   NEWMAT::Matrix m;
   NEWMAT::ColumnVector r((int)f.size());
   for(idx = 0; idx < f.size(); ++idx)
   {
      r[idx] = f[idx];
   }
   setupSystemOfEquations(m, r, x, y, z);
   
   coeff = invert(m.t()*m)*m.t()*r;
}

Member Data Documentation

theHeightAboveMSLFlag

bool ossimRpcSolver::theHeightAboveMSLFlag

protected

Definition at line 192 of file ossimRpcSolver.h.

theMaxResidual

ossim_float64 ossimRpcSolver::theMaxResidual

protected

Definition at line 194 of file ossimRpcSolver.h.

Referenced by getMaxError(), and solveCoefficients().

theMeanResidual

ossim_float64 ossimRpcSolver::theMeanResidual

protected

Definition at line 193 of file ossimRpcSolver.h.

Referenced by getRmsError(), and solveCoefficients().

theRefGeom

ossimRefPtr<ossimImageGeometry> ossimRpcSolver::theRefGeom

protected

Definition at line 195 of file ossimRpcSolver.h.

Referenced by solve(), and solveCoefficients().

theRpcModel

ossimRefPtr<ossimRpcModel> ossimRpcSolver::theRpcModel

protected

Definition at line 196 of file ossimRpcSolver.h.

Referenced by evalPoint(), getNitfRpcBTag(), and solveCoefficients().

theUseElevationFlag

bool ossimRpcSolver::theUseElevationFlag

protected

Definition at line 191 of file ossimRpcSolver.h.

Referenced by solveCoefficients().

---

The documentation for this class was generated from the following files:

• ossimRpcSolver.h
• ossimRpcSolver.cpp