1 //*****************************************************************************
2 // FILE: ossimCoarseGridModel.cc
3 //
4 // License:  See LICENSE.txt file in the top level directory.
5 //
6 // AUTHOR: Oscar Kramer
7 //
8 // DESCRIPTION:
9 //   Contains implementation of class ossimCoarseGridModel. This is an
10 //   implementation of an interpolation sensor model.
11 //
12 //   IMPORTANT: The lat/lon grid is for ground points on the ellipsoid.
13 //   The dLat/dHgt and dLon/dHgt partials therefore are used against
14 //   elevations relative to the ellipsoid.
15 //
16 //*****************************************************************************
17 //  $Id: ossimCoarseGridModel.cpp 22825 2014-07-07 23:14:52Z dburken $
18 
19 #include <ossim/projection/ossimCoarseGridModel.h>
20 
21 RTTI_DEF1(ossimCoarseGridModel, "ossimCoarseGridModel", ossimSensorModel);
22 
23 #include <ossim/base/ossimGpt.h>
24 #include <ossim/base/ossimDpt.h>
25 #include <ossim/base/ossimKeywordlist.h>
26 #include <ossim/base/ossimKeywordNames.h>
27 #include <ossim/base/ossimDatumFactory.h>
28 #include <ossim/base/ossimNotifyContext.h>
29 #include <ossim/elevation/ossimElevManager.h>
30 #include <ossim/imaging/ossimImageGeometry.h>
31 #include <ossim/support_data/ossimSupportFilesList.h>
32 #include <ossim/projection/ossimProjectionFactoryRegistry.h>
33 #include <ossim/projection/ossimBilinearProjection.h>
34 #include <cstdio>
35 #include <fstream>
36 
37 //***
38 // Define Trace flags for use within this file:
39 //***
40 #include <ossim/base/ossimTrace.h>
41 
42 using namespace std;
43 
44 static ossimTrace traceExec  ("ossimCoarseGridModel:exec");
45 static ossimTrace traceDebug ("ossimCoarseGridModel:debug");
46 
47 static const char* MODEL_TYPE = "ossimCoarseGridModel";
48 static const char* GRID_FILE_NAME_KW = "grid_file_name";
49 static const char* CROSSES_DATELINE_KW = "crosses_dateline";
50 
51 const ossimFilename DEFAULT_GEOM_FILE_EXT ("geom");
52 const ossimFilename DEFAULT_GRID_FILE_EXT ("ocg");
53 double ossimCoarseGridModel::theInterpolationError = .1;
54 ossim_int32 ossimCoarseGridModel::theMinGridSpacing     = 100;
55 
56 //*****************************************************************************
57 //  DEFAULT CONSTRUCTOR: ossimCoarseGridModel()
58 //
59 //*****************************************************************************
ossimCoarseGridModel()60 ossimCoarseGridModel::ossimCoarseGridModel()
61    :
62       ossimSensorModel(),
63       theDlatDparamGrid (0),
64       theDlonDparamGrid (0),
65       theHeightEnabledFlag(true)
66 {
67    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG)
68       << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel: entering..."
69       << std::endl;
70 
71    theLatGrid.setDomainType(ossimDblGrid::SAWTOOTH_90);
72    theLonGrid.setDomainType(ossimDblGrid::WRAP_180);
73    theLatGrid.enableExtrapolation();
74    theLonGrid.enableExtrapolation();
75 
76    setErrorStatus();
77 
78    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG)
79       << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel: returning..."
80       << std::endl;
81 }
82 
83 //*****************************************************************************
84 //  COPY CONSTRUCTOR: ossimCoarseGridModel(ossimCoarseGridModel)
85 //
86 //*****************************************************************************
ossimCoarseGridModel(const ossimCoarseGridModel & model)87 ossimCoarseGridModel::ossimCoarseGridModel(const ossimCoarseGridModel& model)
88    :
89       ossimSensorModel  (model),
90       theGridFilename   (model.theGridFilename),
91       theLatGrid        (model.theLatGrid),
92       theLonGrid        (model.theLonGrid),
93       theDlatDhGrid     (model.theDlatDhGrid),
94       theDlonDhGrid     (model.theDlonDhGrid),
95       theDlatDparamGrid (0),
96       theDlonDparamGrid (0),
97       theHeightEnabledFlag(true)
98 {
99    if (traceExec()) ossimNotify(ossimNotifyLevel_DEBUG)
100       << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel(model): entering..."
101       << std::endl;
102 
103    int numberOfParams = getNumberOfAdjustableParameters();
104    if(numberOfParams)
105    {
106       //***
107       // Allocate adjustable parameter partials grids then assign:
108       //***
109       theDlatDparamGrid = new ossimDblGrid [numberOfParams];
110       theDlonDparamGrid = new ossimDblGrid [numberOfParams];
111 
112       for (int i=0; i<numberOfParams; i++)
113       {
114          theDlatDparamGrid[i] = model.theDlatDparamGrid[i];
115          theDlonDparamGrid[i] = model.theDlonDparamGrid[i];
116       }
117    }
118 
119    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG)
120       << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel: returning..."
121       << std::endl;
122 }
123 
124 //*****************************************************************************
125 //  CONSTRUCTOR: ossimCoarseGridModel(filename)
126 //
127 //  Constructs model from geometry file
128 //
129 //*****************************************************************************
ossimCoarseGridModel(const ossimFilename & geom_file)130 ossimCoarseGridModel::ossimCoarseGridModel(const ossimFilename& geom_file)
131    :
132       ossimSensorModel(),
133       theDlatDparamGrid (0),
134       theDlonDparamGrid (0),
135       theHeightEnabledFlag(true)
136 {
137    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel(geom_file): entering..." << std::endl;
138 
139    theLatGrid.setDomainType(ossimDblGrid::SAWTOOTH_90);
140    theLonGrid.setDomainType(ossimDblGrid::WRAP_180);
141    theLatGrid.enableExtrapolation();
142    theLonGrid.enableExtrapolation();
143 
144    ossimKeywordlist kwl;
145    if(geom_file.exists()&&kwl.addFile(geom_file))
146    {
147       loadState(kwl);
148       theGridFilename = geom_file.path();
149    }
150    else
151    {
152       ++theErrorStatus;
153    }
154 
155    if (traceExec())   ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel(geom_file): returning..." << std::endl;
156    return;
157 }
158 
159 //*****************************************************************************
160 //  CONSTRUCTOR: ossimCoarseGridModel(kwl)
161 //
162 //  Constructs model from keywordlist geometry file
163 //
164 //*****************************************************************************
ossimCoarseGridModel(const ossimKeywordlist & geom_kwl)165 ossimCoarseGridModel::ossimCoarseGridModel(const ossimKeywordlist& geom_kwl)
166    :
167       ossimSensorModel(),
168       theDlatDparamGrid (0),
169       theDlonDparamGrid (0),
170       theHeightEnabledFlag(true)
171 {
172    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::ossimCoarseGridModel(geom_kwl): entering..." << std::endl;
173 
174    theLatGrid.setDomainType(ossimDblGrid::SAWTOOTH_90);
175    theLonGrid.setDomainType(ossimDblGrid::WRAP_180);
176    theLatGrid.enableExtrapolation();
177    theLonGrid.enableExtrapolation();
178 
179   // Parse keywordlist for geometry:
180    loadState(geom_kwl);
181 }
182 
183 //*************************************************************************************************
184 //! Assigns the grid data given a projection (typically a rigorous sensor model)
185 //*************************************************************************************************
buildGrid(const ossimDrect & imageBounds,ossimProjection * proj,double heightDelta,bool enableHeightFlag,bool makeAdjustableFlag)186 void ossimCoarseGridModel::buildGrid(const ossimDrect& imageBounds,
187                                      ossimProjection* proj,
188                                      double heightDelta,
189                                      bool enableHeightFlag,
190                                      bool makeAdjustableFlag)
191 {
192    ossimRefPtr<ossimImageGeometry> geom = new ossimImageGeometry();
193    geom->setProjection(proj);
194    buildGrid(imageBounds, geom.get(), heightDelta, enableHeightFlag, makeAdjustableFlag);
195 }
196 
197 
198 //*************************************************************************************************
199 //! Assigns the grid data given a geometry
200 //*************************************************************************************************
buildGrid(const ossimDrect & imageBounds,ossimImageGeometry * geom,double heightDelta,bool enableHeightFlag,bool makeAdjustableFlag)201 void ossimCoarseGridModel::buildGrid(const ossimDrect& imageBounds,
202                                      ossimImageGeometry* geom,
203                                      double heightDelta,
204                                      bool enableHeightFlag,
205                                      bool makeAdjustableFlag)
206 {
207    theHeightEnabledFlag =  enableHeightFlag;
208 
209    if (!geom->getProjection() || imageBounds.hasNans())
210       return;
211 
212    // don't let it get any smaller than 100, 100 pixels
213    // on the input projector
214    //
215    // may want this to be adjusted by outside
216    //
217    const ossimDatum* targetDatum = ossimDatumFactory::instance()->wgs84();
218    ossimIpt gridSize(2,2);
219    ossimDpt gridOrigin(0,0);
220    ossimGpt gpt;
221    ossimGpt gpt2;
222    ossimGpt bilinearGpt;
223    resizeAdjustableParameterArray(0);
224    double normSplit = 1.0;
225    ossimIpt imageSize = ossimIpt(imageBounds.width(), imageBounds.height());
226    double error = 0.0;
227 
228    ossimIpt imageOrigin = imageBounds.ul();
229    ossimDpt spacing ((double)(imageBounds.width()-1)/(gridSize.x-1),
230       (double)(imageBounds.height()-1)/(gridSize.y-1));
231 
232    if(theDlatDparamGrid)
233    {
234       delete [] theDlatDparamGrid;
235       theDlatDparamGrid = NULL;
236    }
237    if(theDlonDparamGrid)
238    {
239       delete [] theDlonDparamGrid;
240       theDlonDparamGrid = NULL;
241    }
242 
243    geom->localToWorld(imageBounds.midPoint(), gpt);
244 
245    do
246    {
247       if(traceDebug())
248       {
249          ossimNotify(ossimNotifyLevel_DEBUG) << "Checking grid size " << gridSize << std::endl;
250       }
251 
252       spacing = ossimDpt((double)(imageBounds.width()-1)/(gridSize.x-1),
253          (double)(imageBounds.height()-1)/(gridSize.y-1));
254 
255       theLatGrid.setNullValue(ossim::nan());
256       theLonGrid.setNullValue(ossim::nan());
257       theDlatDhGrid.setNullValue(0.0);
258       theDlonDhGrid.setNullValue(0.0);
259       theLatGrid.setDomainType(ossimDblGrid::SAWTOOTH_90);
260       theLonGrid.setDomainType(ossimDblGrid::WRAP_180);
261       theLatGrid.initialize(gridSize, gridOrigin, spacing);
262       theLonGrid.initialize(gridSize, gridOrigin, spacing);
263       theDlatDhGrid.initialize(gridSize, gridOrigin, spacing);
264       theDlonDhGrid.initialize(gridSize, gridOrigin, spacing);
265       ossim_int32 x, y;
266 
267       for(y = 0; y < gridSize.y; ++y)
268       {
269          for(x = 0; x < gridSize.x; ++x)
270          {
271             ossimDpt norm((double)x/(double)(gridSize.x-1),
272                (double)y/(double)(gridSize.y-1));
273 
274             ossimDpt pt(imageOrigin.x + norm.x*(imageSize.x-1),
275                imageOrigin.y + norm.y*(imageSize.y-1));
276 
277             geom->localToWorld(pt, gpt);
278             double h = gpt.height();
279             if(ossim::isnan(h))
280             {
281                h += heightDelta;
282             }
283             ossimDpt fullPt;
284             geom->rnToFull(pt, 0, fullPt);
285             geom->getProjection()->lineSampleHeightToWorld(fullPt, h, gpt2);
286             gpt.changeDatum(targetDatum);
287             gpt2.changeDatum(targetDatum);
288 
289             theLatGrid.setNode(x, y, gpt.latd());
290             theLonGrid.setNode(x, y, gpt.lond());
291 
292             theDlatDhGrid.setNode(x, y, (gpt2.latd() - gpt.latd())/heightDelta);
293             theDlonDhGrid.setNode(x, y, (gpt2.lond() - gpt.lond())/heightDelta);
294          }
295       }
296       ossim_int32 upperY = 2*gridSize.y;
297       ossim_int32 upperX = 2*gridSize.x;
298       error = 0.0;
299 
300       // Set all base-class data members needed for subsequent calls to projection code:
301       initializeModelParams(imageBounds);
302 
303       for(y = 0; ((y < upperY)&&(error < theInterpolationError)); ++y)
304       {
305          for(x = 0; ((x < upperX)&&(error<theInterpolationError)); ++x)
306          {
307             ossimDpt norm((double)x/(double)(upperX-1),
308                (double)y/(double)(upperY-1));
309 
310             ossimDpt imagePoint(imageOrigin.x + norm.x*(imageSize.x-1),
311                imageOrigin.y + norm.y*(imageSize.y-1));
312             ossimDpt testIpt;
313 
314             geom->localToWorld(imagePoint, gpt);
315             worldToLineSample(gpt, testIpt);
316             error = (testIpt-imagePoint).length();
317          }
318       }
319 
320       gridSize.x *= 2;
321       gridSize.y *= 2;
322       normSplit *= .5;
323    } while((error > theInterpolationError) &&
324            ((imageSize.x*normSplit) > theMinGridSpacing) &&
325            ((imageSize.y*normSplit) > theMinGridSpacing));
326 
327    gridSize = theLatGrid.size();
328 
329    ossimAdjustableParameterInterface* adjustableParameters =
330       PTR_CAST(ossimAdjustableParameterInterface, geom->getProjection());
331    removeAllAdjustments();
332    if(adjustableParameters&&makeAdjustableFlag)
333    {
334       if(adjustableParameters->getNumberOfAdjustableParameters() > 0)
335       {
336          newAdjustment(adjustableParameters->getNumberOfAdjustableParameters());
337 
338          int numberOfParams = getNumberOfAdjustableParameters();
339          if(numberOfParams)
340          {
341             //***
342             // Allocate adjustable parameter partials grids then assign:
343             //***
344             theDlatDparamGrid = new ossimDblGrid [numberOfParams];
345             theDlonDparamGrid = new ossimDblGrid [numberOfParams];
346             for(int paramIdx = 0; paramIdx < numberOfParams; ++ paramIdx)
347             {
348                theDlonDparamGrid[paramIdx].setNullValue(0.0);
349                theDlatDparamGrid[paramIdx].setNullValue(0.0);
350                theDlatDparamGrid[paramIdx].initialize(gridSize, gridOrigin, spacing);
351                theDlonDparamGrid[paramIdx].initialize(gridSize, gridOrigin, spacing);
352                setAdjustableParameter(paramIdx, 0.0);
353                setParameterSigma(paramIdx, adjustableParameters->getParameterSigma(paramIdx));
354                setParameterUnit(paramIdx, adjustableParameters->getParameterUnit(paramIdx));
355                setParameterCenter(paramIdx, 0.0);
356                setParameterDescription(paramIdx,
357                   adjustableParameters->getParameterDescription(paramIdx));
358 
359                double oldParameter = adjustableParameters->getAdjustableParameter(paramIdx);
360                adjustableParameters->setAdjustableParameter(paramIdx, 1.0, true);
361                double adjust = adjustableParameters->computeParameterOffset(paramIdx);
362                double deltaLat = 0;
363                double deltaLon = 0;
364                if(adjust != 0.0)
365                {
366                   for(int y = 0; y < gridSize.y; ++y)
367                   {
368                      for(int x = 0; x < gridSize.x; ++x)
369                      {
370                         ossimDpt norm((double)x/(double)(gridSize.x-1),
371                            (double)y/(double)(gridSize.y-1));
372 
373                         ossimDpt pt(imageOrigin.x + norm.x*(imageSize.x-1),
374                            imageOrigin.y + norm.y*(imageSize.y-1));
375                         geom->localToWorld(pt, gpt);
376 
377                         gpt.changeDatum(targetDatum);
378                         gpt2.latd(theLatGrid(pt));
379                         gpt2.lond(theLonGrid(pt));
380                         deltaLat = gpt.latd()-gpt2.latd();
381                         deltaLon = gpt.lond()-gpt2.lond();
382 
383                         theDlatDparamGrid[paramIdx].setNode(x, y, deltaLat/adjust);
384                         theDlonDparamGrid[paramIdx].setNode(x, y, deltaLon/adjust);
385                      }
386                   }
387 
388                   // The partials grids for this parameter are initialized, now initialize the
389                   // grid's extrapolator:
390                   theDlatDparamGrid[paramIdx].enableExtrapolation();
391                   theDlonDparamGrid[paramIdx].enableExtrapolation();
392                }
393                adjustableParameters->setAdjustableParameter(paramIdx, oldParameter, true);
394             }
395          }
396       }
397    }
398    getAdjustment(theInitialAdjustment);
399 }
400 
setInterpolationError(double error)401 void ossimCoarseGridModel::setInterpolationError(double error)
402 {
403    theInterpolationError = error;
404 }
405 
setMinGridSpacing(ossim_int32 minSpacing)406 void ossimCoarseGridModel::setMinGridSpacing(ossim_int32 minSpacing)
407 {
408    theMinGridSpacing = minSpacing;
409 }
410 
411 //*************************************************************************************************
412 //! Initializes base class data members after grids have been assigned.
413 //! It is assumed that theImageSize and the origin image point were already set.
414 //*************************************************************************************************
initializeModelParams(ossimIrect imageBounds)415 void ossimCoarseGridModel::initializeModelParams(ossimIrect imageBounds)
416 {
417    // NOTE: it is assumed that the grid size and spacing is the same for ALL grids:
418    ossimIpt gridSize (theLatGrid.size());
419    ossimDpt spacing  (theLatGrid.spacing());
420    ossimDpt v[4];
421    v[0].lat = theLatGrid.getNode(0,0);
422    v[0].lon = theLonGrid.getNode(0,0);
423    v[1].lat = theLatGrid.getNode(gridSize.x-1, 0);
424    v[1].lon = theLonGrid.getNode(gridSize.x-1, 0);
425    v[2].lat = theLatGrid.getNode(gridSize.x-1, gridSize.y-2);
426    v[2].lon = theLonGrid.getNode(gridSize.x-1, gridSize.y-2);
427    v[3].lat = theLatGrid.getNode(0, gridSize.y-2);
428    v[3].lon = theLonGrid.getNode(0, gridSize.y-2);
429 
430    // Guaranty longitude values are -180 to 180
431    for (int i=0; i<4; i++)
432    {
433       if (v[i].lon > 180.0)
434          v[i].lon -= 360.0;
435    }
436 
437    theBoundGndPolygon = ossimPolygon(4, v);
438 
439    theImageSize  = ossimDpt(imageBounds.width(), imageBounds.height());
440    theRefImgPt   = imageBounds.midPoint();
441    theRefGndPt.lat = theLatGrid(theRefImgPt);
442    theRefGndPt.lon = theLonGrid(theRefImgPt);
443 
444    ossimDpt ref_ip_dx (theRefImgPt.x+1.0, theRefImgPt.y    );
445    ossimDpt ref_ip_dy (theRefImgPt.x    , theRefImgPt.y+1.0);
446    ossimGpt ref_gp_dx (theLatGrid(ref_ip_dx), theLonGrid(ref_ip_dx));
447    ossimGpt ref_gp_dy (theLatGrid(ref_ip_dy), theLonGrid(ref_ip_dy));
448 
449    theGSD.x   = theRefGndPt.distanceTo(ref_gp_dx);
450    theGSD.y   = theRefGndPt.distanceTo(ref_gp_dy);
451 
452    theMeanGSD = (theGSD.line + theGSD.samp)/2.0;
453    theImageClipRect  = imageBounds;
454    theSubImageOffset = imageBounds.ul();
455 }
456 
457 //*****************************************************************************
458 //  DESTRUCTOR: ~ossimCoarseGridModel()
459 //
460 //*****************************************************************************
~ossimCoarseGridModel()461 ossimCoarseGridModel::~ossimCoarseGridModel()
462 {
463    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG)
464       << "DEBUG ossimCoarseGridModel::~ossimCoarseGridModel: entering..."
465       << std::endl;
466 
467    if(theDlatDparamGrid&&theDlonDparamGrid)
468    {
469       //***
470       // Deallocate memory:
471       //***
472       delete [] theDlatDparamGrid;
473       delete [] theDlonDparamGrid;
474       theDlatDparamGrid = NULL;
475       theDlonDparamGrid = NULL;
476    }
477    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG)
478       << "DEBUG ossimCoarseGridModel::~ossimCoarseGridModel: returning..."
479       << std::endl;
480 }
481 
482 //*************************************************************************************************
483 // METHOD
484 //*************************************************************************************************
lineSampleToWorld(const ossimDpt & image_point,ossimGpt & gpt) const485 void ossimCoarseGridModel::lineSampleToWorld(const ossimDpt& image_point,
486                                              ossimGpt&       gpt) const
487 {
488    if(!theHeightEnabledFlag)
489    {
490       //
491       // Extrapolate if image point is outside image:
492       //
493       if (!insideImage(image_point))
494       {
495          gpt = extrapolate(image_point);
496          return;
497       }
498 
499       lineSampleHeightToWorld(image_point, 0.0, gpt);
500    }
501    else
502    {
503       ossimSensorModel::lineSampleToWorld(image_point, gpt);
504    }
505 }
506 
507 //*****************************************************************************
508 //  METHOD: ossimCoarseGridModel::lineSampleHeightToWorld()
509 //
510 //  Establishes the ground point corresponding to the input image_point and
511 //  specified elevation above MSL
512 //
513 //*****************************************************************************
514 void
lineSampleHeightToWorld(const ossimDpt & lineSampPt,const double & arg_hgt_above_ellipsoid,ossimGpt & worldPt) const515 ossimCoarseGridModel::lineSampleHeightToWorld(const ossimDpt& lineSampPt,
516                                               const double&   arg_hgt_above_ellipsoid,
517                                               ossimGpt&       worldPt) const
518 {
519    if (traceExec()) ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::lineSampleHeightToWorld: entering..." << std::endl;
520 
521    if(theLatGrid.size().x < 1 ||
522       theLatGrid.size().y < 1)
523    {
524       worldPt.makeNan();
525       return;
526    }
527 
528    double height = (ossim::isnan(arg_hgt_above_ellipsoid)) ? 0.0 : arg_hgt_above_ellipsoid;
529 
530    // Note that there is no check for image point outside of valid image rect because this model
531    // uses the extrapolation inherent to the ossimDblGrid.
532 
533    // The image point may correspond to an offset sub-image. Need to convert
534    // to full image space before anything:
535    ossimDpt ip = lineSampPt + theSubImageOffset;
536 
537    // Establish the interpolated values from the grids:
538    worldPt.lat = theLatGrid(ip);
539    worldPt.lon = theLonGrid(ip);
540    worldPt.hgt = height;
541 
542    if(theHeightEnabledFlag)
543    {
544       // Adjust horizontally due to elevation:
545       worldPt.lat += theDlatDhGrid(ip)*height;
546       worldPt.lon += theDlonDhGrid(ip)*height;
547    }
548    int numberOfParams = getNumberOfAdjustableParameters();
549 
550    // Now add increments due to adjustable parameter deltas:
551    for (int p=0; p<numberOfParams; p++)
552    {
553        worldPt.lat += (theDlatDparamGrid[p](ip) * computeParameterOffset(p));
554        worldPt.lon += (theDlonDparamGrid[p](ip) * computeParameterOffset(p));
555    }
556 
557    worldPt.limitLonTo180();
558 
559    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::lineSampleHeightToWorld: returning..." << std::endl;
560 }
561 
562 //*************************************************************************************************
563 // METHOD
564 //*************************************************************************************************
initAdjustableParameters()565 void ossimCoarseGridModel::initAdjustableParameters()
566 {
567    if(getNumberOfAdjustableParameters() < 1)
568    {
569       addAdjustment(theInitialAdjustment, true);
570    }
571    else
572    {
573       setAdjustment(theInitialAdjustment, true);
574    }
575 }
576 
577 //*************************************************************************************************
578 // METHOD
579 //*************************************************************************************************
imagingRay(const ossimDpt & image_point,ossimEcefRay & image_ray) const580 void ossimCoarseGridModel::imagingRay(const ossimDpt& image_point,
581                                       ossimEcefRay&   image_ray) const
582 {
583    ossimSensorModel::imagingRay(image_point, image_ray);
584 }
585 
586 
587 //*****************************************************************************
588 //  METHOD: ossimCoarseGridModel::print()
589 //
590 //  Formatted dump of data members.
591 //
592 //*****************************************************************************
print(std::ostream & out) const593 std::ostream& ossimCoarseGridModel::print(std::ostream& out) const
594 {
595    out << "\nDump of ossimCoarseGridModel object at: " << this << "\n"
596        << "\n           Grid File Name: " << theGridFilename
597        << "\n                 Image ID: " << theImageID
598        << "\n                   Sensor: " << theSensorID
599        << "\n  Image Size (rows, cols): " << theImageSize
600        << "\n      Ref Pt (samp, line): " << theRefImgPt
601        << "\n   Ref Pt (lat, lon, hgt): " << theRefGndPt
602        << "\n           GSD (row, col): " << theGSD
603        << "\n  Bounding Ground Polygon: " << theBoundGndPolygon << endl;
604 //      << "\n         Number of Params: " << theNumParams << "\n"<<endl;
605 
606    char buf[256];
607    ossimIpt size (theLatGrid.size());
608    ossimDpt spacing (theLatGrid.spacing());
609    int line, samp;
610    ossimIpt node;
611 
612    out << "[ line,  samp]        lat        lon         dLat/dH      dLon/dH\n"
613        << "-------------------------------------------------------------------"
614        <<endl;
615 
616    for (node.y=0; node.y<size.y; node.y++)
617    {
618       line = (int) (node.y*spacing.y);
619 
620       for (node.x=0; node.x<size.x; node.x++)
621       {
622          samp = (int) (node.x*spacing.x);
623 
624          sprintf(buf, "[%5d, %5d]    %+9.5f  %+10.5f    %+11.4e  %+11.4e",
625                   line, samp,
626                   theLatGrid.getNode(node),
627                   theLonGrid.getNode(node),
628                   theDlatDhGrid.getNode(node),
629                   theDlonDhGrid.getNode(node));
630          out << buf << endl;
631       }
632       out <<"-----------------------------------------------------------------"
633           <<endl;
634    }
635 
636    out << "\n\nDump of lat/lon Partials w.r.t. Adjustable Parameters:"<<endl;
637    out << "\nEnd Dump of ossimCoarseGridModel.\n" <<  endl;
638    return out;
639 }
640 
641 //*****************************************************************************
642 //  METHOD: ossimCoarseGridModel::saveState()
643 //
644 //  Saves the model state to the KWL. This KWL also serves as a geometry file.
645 //
646 //  Returns true if successful.
647 //
648 //*****************************************************************************
saveState(ossimKeywordlist & kwl,const char * prefix) const649 bool ossimCoarseGridModel::saveState(ossimKeywordlist& kwl,
650                                      const char* prefix) const
651 {
652    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::saveState: entering..." << std::endl;
653 
654    kwl.add(prefix, ossimKeywordNames::TYPE_KW, "ossimCoarseGridModel");
655 
656    //---
657    // Save just the file part so that the geometry stuff is relocatable.
658    // Code was added to ossimProjectionFactoryBase::createProjectionFromGeometryFile
659    // to handle this.
660    //---
661    kwl.add( prefix, GRID_FILE_NAME_KW, theGridFilename.file() );
662 
663    kwl.add(prefix, "height_enabled_flag", theHeightEnabledFlag, true);
664    ossimSensorModel::saveState(kwl, prefix);
665    ossimString initAdjPrefix = ossimString(prefix) + "init_adjustment.";
666    theInitialAdjustment.saveState(kwl, initAdjPrefix);
667 
668    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::saveState: returning..." << std::endl;
669 
670    return true;
671 }
672 
673 //*****************************************************************************
674 //  METHOD: ossimCoarseGridModel::loadState()
675 //
676 //  Restores the model's state from the KWL. This KWL also serves as a
677 //  geometry file.
678 //
679 //  Returns true if successful.
680 //
681 //*****************************************************************************
loadState(const ossimKeywordlist & kwl,const char * prefix)682 bool ossimCoarseGridModel::loadState(const ossimKeywordlist& kwl,
683                                      const char* prefix)
684 {
685    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadState: entering..." << std::endl;
686 
687    clearErrorStatus();
688    const char* value;
689    bool success;
690 
691    //***
692    // Assure this keywordlist contains correct type info:
693    //***
694    value = kwl.find(prefix, ossimKeywordNames::TYPE_KW);
695    if (!value || (strcmp(value, "ossimCoarseGridModel")))
696    {
697       if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadState:  returning..." << std::endl;
698       theErrorStatus++;
699       return false;
700    }
701    value = kwl.find(prefix, "height_enabled_flag");
702    if(value)
703    {
704       theHeightEnabledFlag = ossimString(value).toBool();
705    }
706 
707    //***
708    // Pass on to the base-class for parsing first:
709    //***
710    success = ossimSensorModel::loadState(kwl, prefix);
711    if (!success)
712    {
713       theErrorStatus++;
714       if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadState: returning with error..." << std::endl;
715       return false;
716    }
717 
718    //***
719    // Look for geom filename or explicit grid filename to establish path to grid:
720    //***
721    theGridFilename = kwl.find(prefix, GRID_FILE_NAME_KW);
722 
723    if (!theGridFilename.isReadable())
724    {
725       //---
726       // This is set in ossimProjectionFactoryBase::createProjectionFromGeometryFile
727       // so we can derive the dot.ocg from it.
728       //---
729       ossimFilename alt_path_to_grid = kwl.find("kwl_source");
730       theGridFilename = alt_path_to_grid.setExtension(DEFAULT_GRID_FILE_EXT);
731 
732       if (!theGridFilename.isReadable())
733       {
734          ossimFilename alt_path_to_grid = kwl.find(prefix, ossimKeywordNames::GEOM_FILE_KW);
735          theGridFilename = alt_path_to_grid.setExtension(DEFAULT_GRID_FILE_EXT);
736       }
737    }
738 
739    if (!theGridFilename.isReadable())
740    {
741       ossimNotify(ossimNotifyLevel_FATAL) << "ossimCoarseGridModel::loadState() -- Error "
742          "encountered opening coarse grid file at "<< "<" <<theGridFilename << ">." << std::endl;
743       if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadState: returning with error..." << std::endl;
744       theErrorStatus++;
745       return false;
746    }
747 
748    // Load the coarse grid file specified in KWL. This method resets the
749    // theErrorStatus to OSSIM_OK if successful:
750    if (!loadCoarseGrid(theGridFilename))
751    {
752       theErrorStatus++;
753       return false;
754    }
755 
756    // crossesDateline legacy. No longer saved.
757    bool crossesDateline = false;
758    kwl.getBoolKeywordValue(crossesDateline, CROSSES_DATELINE_KW, prefix);
759    if (crossesDateline)
760       theLonGrid.setDomainType(ossimDblGrid::WRAP_360);
761 
762    // Add the coarse grid filename to list of support files being referenced for logging purposes:
763    ossimSupportFilesList::instance()->add(theGridFilename.expand());
764 
765    ossimString initAdjPrefix = ossimString(prefix) + "init_adjustment.";
766    theInitialAdjustment.loadState(kwl, initAdjPrefix.c_str());
767 
768    if((ossim::isnan(theRefGndPt.hgt)) ||
769       (theRefGndPt.hgt == 0))
770    {
771       theRefGndPt.hgt = ossimElevManager::instance()->getHeightAboveEllipsoid(theRefGndPt);
772       if(theRefGndPt.hgt < 0)
773       {
774          theRefGndPt.hgt = fabs(theRefGndPt.hgt);
775       }
776    }
777 
778    if(theInitialAdjustment.getNumberOfAdjustableParameters() < 1)
779    {
780       getAdjustment(theInitialAdjustment);
781    }
782    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadState: returning..." << std::endl;
783    if (theErrorStatus)
784       return false;
785 
786    return true;
787 }
788 
789 //*****************************************************************************
790 //  METHOD: ossimCoarseGridModel::saveCoarseGrid(cgFileName)
791 //
792 //  Saves the coarse grid to disk file.
793 //
794 //  Returns true if successful.
795 //
796 //*****************************************************************************
saveCoarseGrid(const ossimFilename & fileName) const797 bool ossimCoarseGridModel::saveCoarseGrid(const ossimFilename& fileName)const
798 {
799    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::saveCoarseGrid: entering..." << std::endl;
800 
801    // Create and open grid file as stream:
802    theGridFilename = fileName.expand();
803    theGridFilename.setExtension(DEFAULT_GRID_FILE_EXT);
804    ofstream outstream (theGridFilename.chars());
805    if (!outstream.is_open())
806    {
807       theErrorStatus++;
808       ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::saveCoarseGrid: Error "
809          "encountered creating coarse grid file <" << theGridFilename<< ">. Check that directory "
810          "exists and is writable." << std::endl;
811       if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::saveCoarseGrid: returning with error..." << std::endl;
812       return false;
813    }
814 
815    // Let each grid object write itself to the output file:
816    theLatGrid.save(outstream, "Latitude Grid");
817    theLonGrid.save(outstream, "Longitude Grid");
818    theDlatDhGrid.save(outstream, "dLat/dH Grid");
819    theDlonDhGrid.save(outstream, "dLon_dH Grid");
820 
821    ossimString descr;
822    int numberOfParams = getNumberOfAdjustableParameters();
823    for (int p=0; p<numberOfParams; p++)
824    {
825       descr = getParameterDescription(p) + " dLat_dParam Grid";
826       theDlatDparamGrid[p].save(outstream, descr.chars());
827       descr = getParameterDescription(p) + " dLon_dParam Grid";
828       theDlonDparamGrid[p].save(outstream, descr.chars());
829    }
830 
831    // Since the geom file is needed in the same path as the grid file, take this opportunity to
832    // write the geom file out as well:
833    ossimFilename geom_file (theGridFilename);
834    geom_file.setExtension(DEFAULT_GEOM_FILE_EXT);
835    ossimKeywordlist kwl;
836    saveState(kwl);
837    kwl.write(geom_file);
838 
839    // Add to the list of support files referenced (though technically it has not yet been
840    // referenced, but will be next time this image is opened):
841    ossimSupportFilesList::instance()->add(geom_file);
842    ossimSupportFilesList::instance()->add(theGridFilename);
843 
844    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::saveCoarseGrid: returning..." << std::endl;
845    return true;
846 }
847 
848 
849 //*****************************************************************************
850 //  METHOD: ossimCoarseGridModel::loadCoarseGrid(cgFileName)
851 //
852 //  Loads the coarse grid from disk file.
853 //
854 //  Returns true if successful.
855 //
856 //*****************************************************************************
loadCoarseGrid(const ossimFilename & cgFileName)857 bool ossimCoarseGridModel::loadCoarseGrid(const ossimFilename& cgFileName)
858 {
859    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadCoarseGrid: entering..." << std::endl;
860 
861    ossimDpt v[4];
862    ossimIpt grid_size;
863 
864    //***
865    // Open existing grid file:
866    //***
867    ifstream instream (cgFileName.chars());
868    if (!instream.is_open())
869    {
870       theErrorStatus++;
871       if(traceDebug())
872       {
873          ossimNotify(ossimNotifyLevel_FATAL) << "FATAL ossimCoarseGridModel::loadCoarseGrid: Error encountered opening coarse grid file <" << cgFileName
874          << ">. Check that the file exists and is readable." << std::endl;
875       }
876       if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "CEBUG ossimCoarseGridModel::loadCoarseGrid: returning with error..." << std::endl;
877       return false;
878    }
879    theGridFilename = cgFileName;
880    if(theDlatDparamGrid)
881      {
882        delete [] theDlatDparamGrid;
883        theDlatDparamGrid = NULL;
884      }
885    if(theDlonDparamGrid)
886      {
887        delete [] theDlonDparamGrid;
888        theDlonDparamGrid = NULL;
889      }
890    //
891    // Reallocate memory:
892    //
893    int numberOfParams = getNumberOfAdjustableParameters();
894    if(numberOfParams)
895    {
896       theDlatDparamGrid = new ossimDblGrid [numberOfParams];
897       theDlonDparamGrid = new ossimDblGrid [numberOfParams];
898    }
899    //***
900    // Let each grid object read itself from the input file:
901    //***
902    if (!theLatGrid.load(instream))
903    {
904       ++theErrorStatus;
905       return false;
906    }
907    if (!theLonGrid.load(instream))
908    {
909       ++theErrorStatus;
910       return false;
911    }
912    if (!theDlatDhGrid.load(instream))
913    {
914       ++theErrorStatus;
915       return false;
916    }
917    if (!theDlonDhGrid.load(instream))
918    {
919       ++theErrorStatus;
920       return false;
921    }
922 
923    for (int p=0; p<numberOfParams; p++)
924    {
925       if (!theDlatDparamGrid[p].load(instream))
926       {
927          ++theErrorStatus;
928          return false;
929       }
930       if (!theDlonDparamGrid[p].load(instream))
931       {
932          ++theErrorStatus;
933          return false;
934       }
935    }
936 
937    //***
938    // Initialize the bounding ground rectangle:
939    //***
940    grid_size = theLatGrid.size();
941 
942    v[0].lat = theLatGrid(0,0);
943    v[0].lon = theLonGrid(0,0);
944    v[1].lat = theLatGrid(theImageSize.x-1, 0);
945    v[1].lon = theLonGrid(theImageSize.x-1, 0);
946    v[2].lat = theLatGrid(theImageSize.x-1, theImageSize.y-1);
947    v[2].lon = theLonGrid(theImageSize.x-1, theImageSize.y-1);
948    v[3].lat = theLatGrid(0, theImageSize.y-1);
949    v[3].lon = theLonGrid(0, theImageSize.y-1);
950 
951    theBoundGndPolygon = ossimPolygon(4, v);
952 
953    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::loadCoarseGrid: returning..." << std::endl;
954    return true;
955 }
956 
957 
958 //*****************************************************************************
959 //  METHOD: ossimCoarseGridModel::reallocateGrid()
960 //
961 //  Deletes existing grid arrays and allocates new ones.
962 //
963 //*****************************************************************************
reallocateGrid(const ossimIpt & grid_size)964 void ossimCoarseGridModel::reallocateGrid(const ossimIpt& grid_size)
965 {
966    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::reallocateGrid:entering..." << endl;
967 
968 
969    //***
970    // Deallocate memory:
971    //***
972    if(theDlatDparamGrid)
973      {
974        delete [] theDlatDparamGrid;
975        theDlatDparamGrid = NULL;
976      }
977    if(theDlonDparamGrid)
978      {
979        delete [] theDlonDparamGrid;
980        theDlonDparamGrid = NULL;
981      }
982    //***
983    // determine grid spacing given new info:
984    //***
985    ossimDpt spacing ((double)(theImageSize.x-1)/(double)(grid_size.x-1),
986                      (double)(theImageSize.y-1)/(double)(grid_size.y-1));
987 
988    //***
989    // Allocate all:
990    //***
991    ossimDpt grid_origin(0.0, 0.0);
992    theLatGrid.setNullValue(ossim::nan());
993    theLonGrid.setNullValue(ossim::nan());
994    theDlatDhGrid.setNullValue(0.0);
995    theDlonDhGrid.setNullValue(0.0);
996    theLatGrid.initialize(grid_size, grid_origin, spacing);
997    theLonGrid.initialize(grid_size, grid_origin, spacing);
998    theDlatDhGrid.initialize(grid_size, grid_origin, spacing);
999    theDlonDhGrid.initialize(grid_size, grid_origin, spacing);
1000 
1001    int numberOfParams = getNumberOfAdjustableParameters();
1002    if(numberOfParams)
1003    {
1004 
1005       theDlatDparamGrid = new ossimDblGrid [numberOfParams];
1006       theDlonDparamGrid = new ossimDblGrid [numberOfParams];
1007    }
1008    for (int p=0; p<numberOfParams; p++)
1009    {
1010       theDlonDparamGrid[p].setNullValue(0.0);
1011       theDlatDparamGrid[p].setNullValue(0.0);
1012       theDlatDparamGrid[p].initialize(grid_size, grid_origin, spacing);
1013       theDlonDparamGrid[p].initialize(grid_size, grid_origin, spacing);
1014    }
1015 
1016    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::reallocateGrid: returning..." << std::endl;
1017    return;
1018 }
1019 
1020 //*****************************************************************************
1021 // STATIC METHOD: ossimCoarseGridModel::writeGeomTemplate
1022 //
1023 //  Writes a sample geometry KWL to the output stream.
1024 //
1025 //*****************************************************************************
writeGeomTemplate(ostream & os)1026 void ossimCoarseGridModel::writeGeomTemplate(ostream& os)
1027 {
1028    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::writeGeomTemplate: entering..." << std::endl;
1029 
1030    os <<
1031       "//**************************************************************\n"
1032       "// Template for OCG model kewordlist\n"
1033       "//**************************************************************\n"
1034       << ossimKeywordNames::TYPE_KW << ": " << MODEL_TYPE << endl;
1035 
1036    ossimSensorModel::writeGeomTemplate(os);
1037 
1038    os << "//\n"
1039       << "// Derived-class ossimCoarseGridModel Keywords:\n"
1040       << "//\n"
1041       << GRID_FILE_NAME_KW << ": <string>\n" << endl;
1042 
1043    if (traceExec())  ossimNotify(ossimNotifyLevel_DEBUG) << "DEBUG ossimCoarseGridModel::writeGeomTemplate: returning..." << std::endl;
1044    return;
1045 }
1046 
1047 //*************************************************************************************************
1048 // Overrides base-class extrapolation code. Uses extrapolation inherent to ossimDbleGrid
1049 //*************************************************************************************************
extrapolate(const ossimDpt & local_ip,const double & height) const1050 ossimGpt ossimCoarseGridModel::extrapolate(const ossimDpt& local_ip, const double& height) const
1051 {
1052    ossimGpt gpt;
1053    lineSampleHeightToWorld(local_ip, height, gpt);
1054    return gpt;
1055 }
1056 
isAffectedByElevation() const1057 bool ossimCoarseGridModel::isAffectedByElevation() const
1058 {
1059    return theHeightEnabledFlag;
1060 }
1061