1 //*******************************************************************
2 //
3 // License: See top level LICENSE.txt file.
4 //
5 // Author: Oscar Kramer
6 //
7 // Description:
8 //
9 // Sensor Model for Radarsat1 SAR sensor.
10 //
11 //*******************************************************************
12 // $Id:$
13
14 #include <ossim/projection/ossimRS1SarModel.h>
15 #include <cstdio>
16 #include <fstream>
17 #include <iomanip>
18 #include <iostream>
19 #include <sstream>
20 #include <string>
21 #include <ossim/base/ossimTrace.h>
22 #include <ossim/matrix/newmat.h>
23 #include <ossim/base/ossimLsrPoint.h>
24 #include <ossim/elevation/ossimElevManager.h>
25
26
27 using namespace std;
28 using namespace ossim;
29
30 static ossimTrace traceDebug("ossimRS1SarModel:debug");
31
32 static const double SEC_PER_DAY = 86400.0;
33 static const double NOMINAL_ORBIT_ELEVATION = 798000.0; // m
34 static const double DEG_PER_SEC = 0.004178074;
35 static const double INTRACK_OFFSET_SIGMA = 1000.0; // meters
36 static const double CRTRACK_OFFSET_SIGMA = 300.0;
37 static const double RADIAL_OFFSET_SIGMA = 300.0;
38 static const double LINE_SCALE_SIGMA = 0.01;
39 static const double SKEW_SIGMA = 0.1; // deg
40 static const double ORIENTATION_SIGMA = 0.2; // deg
41 static const double EARTH_ANGULAR_VELOCITY = 7.2921151467e-05; // radians/sec
42 static const char* IMAGING_MODE_ID[] =
43 {"UNKNOWN", "SCN", "SCW", "SGC", "SGF", "SGX", "SLC", "SPG", "SSG","RAW","ERS"};
44
45 //******************************************************************************
46 // Default.
47 //******************************************************************************
ossimRS1SarModel()48 ossimRS1SarModel::ossimRS1SarModel()
49 : theCeosData (0),
50 theImagingMode (UNKNOWN_MODE)
51 {
52 setErrorStatus();
53 }
54
55 //******************************************************************************
56 // Takes a ossimFilename& directory containing support data, initializes the model,
57 // and writes out the initial geometry file.
58 //******************************************************************************
ossimRS1SarModel(const ossimFilename & imageFile)59 ossimRS1SarModel::ossimRS1SarModel(const ossimFilename& imageFile)
60 :
61 theCeosData (0),
62 theImagingMode (UNKNOWN_MODE)
63 {
64 static const char MODULE[] = "Constructor ossimRS1SarModel(ossimFilename)";
65 if (traceDebug()) CLOG << "entering..." << endl;
66
67 clearErrorStatus();
68
69 // Parse the CEOS data files:
70 ossimFilename dataDirName (imageFile.path());
71 initFromCeos(dataDirName);
72 if(getErrorStatus() != ossimErrorCodes::OSSIM_OK)
73 {
74 return;
75 }
76 if (theImagingMode == UNKNOWN_MODE)
77 return;
78
79 // Parse the image data file for the local ORPs interpolator:
80 if ((theImagingMode == SCN) || (theImagingMode == SCW))
81 establishOrpGrid();
82 else
83 establishOrpInterp();
84
85 // Establish the ephemeris interpolator:
86 establishEphemeris();
87
88 initAdjParms();
89 establishVehicleSpace();
90
91 if (traceDebug()) CLOG << "returning..." << endl;
92 return;
93 }
94
95
96 //******************************************************************************
97 // DESTRUCTOR: ~ossimRS1SarModel
98 //******************************************************************************
~ossimRS1SarModel()99 ossimRS1SarModel::~ossimRS1SarModel()
100 {
101 deallocateMemory();
102 }
103
104 //******************************************************************************
105 // loadSate()
106 //******************************************************************************
loadState(const ossimKeywordlist &,const char *)107 bool ossimRS1SarModel::loadState(const ossimKeywordlist& /* kwl */, const char* /* prefix */)
108 {
109 // NOT YET IMPLEMENTED
110 setErrorStatus();
111 return false;
112 }
113
114 //******************************************************************************
115 // saveState()
116 //******************************************************************************
saveState(ossimKeywordlist &,const char *) const117 bool ossimRS1SarModel::saveState(ossimKeywordlist& /* kwl */, const char* /* prefix */) const
118 {
119 // NOT YET IMPLEMENTED
120 setErrorStatus();
121 return false;
122 }
123
124 //*************************************************************************************************
125 //! Given an image point, returns a ray originating at some arbitrarily high
126 //! point (in this model at the sensor position) and pointing towards the target.
127 //*************************************************************************************************
imagingRay(const ossimDpt & image_point,ossimEcefRay &) const128 void ossimRS1SarModel::imagingRay(const ossimDpt& image_point, ossimEcefRay& /* image_ray */) const
129 {
130 const char* MODULE = "ossimRS1SarModel::imagingRay()";
131
132 // Apply offset, scale, skew, and rotation to image point:
133 ossimDpt offset (theImageSize.line/2.0, theImageSize.samp/2.0);
134 ossimDpt ip1 = image_point - offset;
135 ossimDpt ip2;
136 ip2.line = ip1.line*theCosOrientation*theLineScale +
137 ip1.samp*(theCosOrientation*theSinSkew+theSinOrientation);
138 ip2.samp = -ip1.line*theSinOrientation*theLineScale -
139 ip1.samp*(theSinOrientation*theSinSkew-theCosOrientation);
140
141 // Offset point back to origin of image:
142 ip1 = ip2 + offset;
143
144 // Given new line number, obtain the interpolated ORP:
145 ossimEcefPoint localORP;
146 if ((theImagingMode == SCN) || (theImagingMode == SCW))
147 {
148 ossimGpt gpt (theLatGrid(ip1.x, ip1.y), theLonGrid(ip1.x, ip1.y), theRefHeight);
149 localORP = ossimEcefPoint(gpt);
150 }
151 else
152 {
153 NEWMAT::ColumnVector orpPos(3);
154 theLocalOrpInterp->interpolate(ip1.line, orpPos);
155 localORP = ossimEcefPoint(orpPos[0], orpPos[1], orpPos[2]);
156 }
157
158 // Establish imaging line time (zero-Doppler time):
159 double arpTime = theFirstLineTime + ip1.line*theTimePerLine;
160
161 NEWMAT::ColumnVector arpEph(3);
162 theArpPosInterp->interpolate(arpTime, arpEph);
163 ossimEcefPoint ecfArpPos (arpEph[0], arpEph[1], arpEph[2]);
164
165 theArpVelInterp->interpolate(arpTime, arpEph);
166 ossimEcefVector ecfArpVel (arpEph[0], arpEph[1], arpEph[2]);
167
168 // Need to correct the velocity vector by the earth rotational velocity:
169 ossimEcefVector earthVel
170 (-ecfArpPos.y()*EARTH_ANGULAR_VELOCITY, ecfArpPos.x()*EARTH_ANGULAR_VELOCITY, 0.0);
171 ecfArpVel = ecfArpVel - earthVel;
172
173 // Compute a vehicle space with Z near intrack, X normal to Z and range
174 // vector to ORP (slant range normal), and Y in range direction.
175 // Note the application of the position adjustable parameter to the space's
176 // origin (after the range vector was established):
177 ossimEcefVector rangeVector = localORP - ecfArpPos;
178 ossimEcefVector rangeNormal = rangeVector.cross(ecfArpVel);
179 // ossimEcefVector zDirection = rangeNormal.cross(rangeVector);
180 ossimLsrSpace localVehicleSpace (ossimEcefPoint(ecfArpPos + thePosCorrection), // origin
181 rangeNormal, // X-direction
182 rangeVector, 0); // Y-direction (Z not specified
183 double local_orp_range = rangeVector.length();
184
185 // Compute slant range distance to pixel in question:
186 double slant_range = local_orp_range;
187 switch (theImagingMode)
188 {
189 case SLC:
190 if (theDirectionFlag == ASCENDING)
191 slant_range += thePixelSpacing.samp*ip1.samp;
192 else
193 slant_range += (theImageSize.samp-1.0-ip1.samp) * thePixelSpacing.samp;
194 break;
195
196 case SGC:
197 case SGF:
198 case SGX:
199 case ERS:
200 {
201 // Compute the slant range as a polynomial expansion given ground range
202 double ground_range;
203 if (theDirectionFlag == ASCENDING)
204 ground_range = ip1.samp * thePixelSpacing.samp;
205 else
206 ground_range = (theImageSize.samp - 1.0 - ip1.samp) *
207 thePixelSpacing.samp;
208 // double delta_g_i = 1.0;
209 // for (int 0=1; i<6; i++)
210 double delta_g_i = ground_range;
211 for (int i=1; i<6; i++) // NOTE: not using offset (i=0) term
212 {
213 slant_range += theSrGrCoeff[i] * delta_g_i;
214 delta_g_i *= ground_range;
215 }
216 }
217 break;
218
219 case SCN:
220 case SCW:
221 slant_range = local_orp_range;
222 break;
223
224 default:
225 CLOG << "ERROR: Invalid imaging mode encountered." << endl;
226 setErrorStatus();
227 }
228 }
229
230 //*************************************************************************************************
231 // PUBLIC METHOD: ossimRS1SarModel::lineSampleHeightToWorld(image_point, height, is_inside_image)
232 //
233 // Performs image to ground projection.
234 //*************************************************************************************************
lineSampleHeightToWorld(const ossimDpt & image_point,const double & height_ellip,ossimGpt & worldPt) const235 void ossimRS1SarModel::lineSampleHeightToWorld(const ossimDpt& image_point,
236 const double& height_ellip,
237 ossimGpt& worldPt) const
238 {
239 static const char MODULE[] = "ossimRS1SarModel::lineSampleHeightToWorld()";
240 if (traceDebug()) CLOG << "entering..." << endl;
241
242 // static const int MAX_NUM_ITERS = 50;
243 // static const double MAX_ELEV_DIFF = 0.001;
244
245 ossimEcefRay imaging_ray;
246 imagingRay(image_point, imaging_ray);
247 worldPt = imaging_ray.intersectAboveEarthEllipsoid(height_ellip);
248
249 if (traceDebug()) CLOG << "returning..." << endl;
250 }
251
252 //******************************************************************************
253 // This method sets theImagingMode data member given a character string
254 // abbreviation of the mode.
255 //******************************************************************************
setImagingMode(char * modeStr)256 void ossimRS1SarModel::setImagingMode(char* modeStr)
257 {
258 static const char MODULE[] = "ossimRS1SarModel::setImagingMode(modeStr)";
259 if (traceDebug()) CLOG << "entering..." << endl;
260
261 bool illegal_mode = false;
262
263 if (strcmp(modeStr, "SCN") == 0)
264 {
265 theImagingMode = SCN;
266 }
267 else if (strcmp(modeStr, "SCW") == 0)
268 {
269 theImagingMode = SCW;
270 }
271 else if (strcmp(modeStr, "SGC") == 0)
272 {
273 theImagingMode = SGC;
274 }
275 else if (strcmp(modeStr, "SGF") == 0)
276 {
277 theImagingMode = SGF;
278 }
279 else if (strcmp(modeStr, "SGX") == 0)
280 {
281 theImagingMode = SGX;
282 }
283 else if (strcmp(modeStr, "SLC") == 0)
284 {
285 theImagingMode = SLC;
286 }
287 else if (strcmp(modeStr, "SPG") == 0)
288 {
289 theImagingMode = SPG;
290 illegal_mode = true;
291 }
292 else if (strcmp(modeStr, "SSG") == 0)
293 {
294 theImagingMode = SSG;
295 illegal_mode = true;
296 }
297 else if (strcmp(modeStr, "RAW") == 0)
298 {
299 theImagingMode = RAW;
300 illegal_mode = true;
301 }
302 else if (strcmp(modeStr, "ERS") == 0)
303 {
304 theImagingMode = ERS;
305 }
306 else
307 {
308 theImagingMode = UNKNOWN_MODE;
309 illegal_mode = true;
310 }
311
312 // Test for ERS product:
313 if ((theImagingMode == UNKNOWN_MODE) && theCeosData.valid())
314 {
315 char buf[] = "xxx";
316 strncpy(buf, &(theCeosData->textRec()->product_type[8]), 3);
317 if (strcmp(buf, "ERS") == 0)
318 {
319 theImagingMode = ERS;
320 illegal_mode = false;
321 }
322 }
323
324 // Filter the imaging modes that are not handled:
325 if (illegal_mode)
326 {
327 CLOG << "\n\t ERROR: The imaging mode <"
328 << IMAGING_MODE_ID[(int) theImagingMode]
329 << "> is currently not supported. \n" << endl;
330 theImagingMode = UNKNOWN_MODE;
331 }
332
333 if (traceDebug()) CLOG << "returning..." << endl;
334 }
335
336
337 //******************************************************************************
338 // Initialize various data members from the CEOS records.
339 //******************************************************************************
initFromCeos(const ossimFilename & fname)340 void ossimRS1SarModel::initFromCeos(const ossimFilename& fname)
341 {
342 static const char MODULE[] = "ossimRS1SarModel::initFromCeos()";
343 if (traceDebug()) CLOG << "entering..." << endl;
344
345 // Instantiate a CeosData object:
346 theCeosData = new ossimCeosData(fname);
347
348 if(theCeosData->errorStatus() != ossimErrorCodes::OSSIM_OK) return;
349 const dataset_sum_rec* dsr = theCeosData->dataSetSumRec();
350 const proc_parm_rec* ppr = theCeosData->procParmRec();
351 char buf[1024];
352
353 // Set imaging mode:
354 strncpy(buf, &(theCeosData->volDescRec()->logvol_id[11]), 3);
355 buf[3] = '\0';
356 setImagingMode(buf);
357
358 // Determine whether ascending or descending path acquisition. This flag
359 // dictates whether the first or last pixel of a line is of shorter range:
360 if (theImagingMode == ERS)
361 {
362 theDirectionFlag = DESCENDING;
363 }
364 else
365 {
366 if (dsr->asc_des[0] == 'A')
367 theDirectionFlag = ASCENDING;
368 else if (dsr->asc_des[0] == 'D')
369 theDirectionFlag = DESCENDING;
370 else
371 {
372 CLOG << "ERROR: Direction Flag: " << dsr->asc_des << " not supported"
373 << endl;
374 theDirectionFlag = UNKNOWN_DIRECTION;
375 setErrorStatus();
376 if (traceDebug()) CLOG << "returning with error..." << endl;
377 return;
378 }
379 }
380
381 // Image (Product) ID:
382 buf[8] = '\0';
383 strncpy(buf, theCeosData->volDescRec()->product_id, 8);
384 theImageID = buf;
385
386 // pixel spacing in range:
387 buf[16] = '\0';
388 strncpy(buf, dsr->pix_spacing, 16);
389 thePixelSpacing.samp = atof(buf);
390
391 // pixel spacing in azimuth:
392 strncpy(buf, dsr->line_spacing, 16);
393 buf[16] = '\0';
394 thePixelSpacing.line = atof(buf);
395
396 // The ground reference point (ORP) latitude:
397 ossimGpt grp;
398 strncpy(buf, dsr->pro_lat, 16);
399 buf[16] = '\0';
400 grp.lat = atof(buf);
401
402 // The ground reference point (ORP) longitude:
403 strncpy(buf, dsr->pro_long, 16);
404 buf[16] = '\0';
405 grp.lon = atof(buf);
406
407 // The ground reference point (ORP) elevation:
408 strncpy(buf, dsr->terrain_h, 16);
409 buf[16] = '\0';
410 grp.hgt = atof(buf);
411 theORP = ossimEcefPoint(grp);
412 theRefHeight = grp.hgt;
413
414 // Illumination elevation:
415 strncpy(buf, dsr->incident_ang, 8);
416 buf[8] = '\0';
417 double incidence = atof(buf);
418 theIllumElevation = 90.0 - incidence;
419
420 // Illumination azimuth -- need to determine whether sensor is in normal
421 // right-looking mode or "antarctic" (left-looking) mode:
422 strncpy(buf, dsr->plat_head, 8);
423 theIllumAzimuth = atof(buf);
424 char sensor_orientation[] = "123456789";
425 strncpy(sensor_orientation, ppr->sens_orient, 9);
426 if (strcmp(sensor_orientation, "ANTARCTIC") == 0)
427 {
428 theIllumAzimuth -= 90.0;
429 if (theIllumAzimuth < 0.0) theIllumAzimuth += 360.0;
430 }
431 else
432 {
433 theIllumAzimuth += 90.0;
434 if (theIllumAzimuth >= 360.0) theIllumAzimuth -= 360.0;
435 }
436
437 // Determine the number of lines per image and pixels per line. NOTE: if
438 // image is scan mode, the number of lines must be computed indirectly.
439 if ((theImagingMode == SCN) || (theImagingMode == SCW))
440 {
441 FILE* fptr = fopen(theCeosData->imageFile().chars(), "r");
442 if (!fptr)
443 {
444 CLOG << " ERROR:\n\tCannot open CEOS image file: "
445 << theCeosData->imageFile() << endl;
446 setErrorStatus();
447 return;
448 }
449 fseek(fptr, 0, SEEK_END);
450 long byte_count = ftell(fptr);
451 buf[6] = '\0';
452 strncpy(buf, theCeosData->imopDescRec()->l_dataset, 6);
453 int lengthOfRec = atoi(buf);
454 theImageSize.line = (byte_count - sizeof(*theCeosData->imopDescRec()))
455 /lengthOfRec;
456 fclose(fptr);
457 }
458 else
459 {
460 strncpy(buf, theCeosData->imopDescRec()->n_dataset, 6);
461 buf[6] = '\0';
462 theImageSize.line = atoi(buf);
463 }
464 strncpy(buf, theCeosData->imopDescRec()->nleft, 4);
465 buf[4] = '\0';
466 int left_pixels = atoi(buf);
467 strncpy(buf, theCeosData->imopDescRec()->ngrp, 8);
468 buf[8] = '\0';
469 int data_pixels = atoi(buf);
470 strncpy(buf, theCeosData->imopDescRec()->nright, 4);
471 buf[4] = '\0';
472 int right_pixels = atoi(buf);
473 theImageSize.samp = left_pixels + data_pixels + right_pixels;
474
475 // If ground range product, need to read the ground-to-slant range
476 // conversion coefficients from the proc_parm_rec. NOTE: only single-look
477 // products handled properly. ScanSAR requires reading multiple coefficient
478 // sets. See MDA Detailed Processing Parameter Record Description sec.3.84.
479 if ((theImagingMode==SGC) ||
480 (theImagingMode==SGF) ||
481 (theImagingMode==SGX) ||
482 (theImagingMode==ERS))
483 {
484 buf[16] = '\0';
485 for (int i=0; i<6; i++)
486 {
487 strncpy(buf, ppr->srgr_coefset[0].srgr_coef[i], 16);
488 theSrGrCoeff[i] = atof(buf);
489 }
490 }
491 else
492 {
493 theSrGrCoeff[0] = 0.0;
494 theSrGrCoeff[1] = 1.0; // note 1.0 here (should never be accessed)
495 theSrGrCoeff[2] = 0.0;
496 theSrGrCoeff[3] = 0.0;
497 theSrGrCoeff[4] = 0.0;
498 theSrGrCoeff[5] = 0.0;
499 }
500
501 if (traceDebug())
502 {
503 CLOG << "DEBUG -- "
504 << "\n\t theDirectionFlag = " << theDirectionFlag
505 << "\n\t thePixelSpacing = " << thePixelSpacing
506 << "\n\t theORP = " << theORP
507 << "\n\t theRefHeight = " << theRefHeight
508 << "\n\t theImageSize = " << theImageSize
509 << "\n\t sensor_orientation = " << sensor_orientation
510 << "\n\t theIllumElevation = " << theIllumElevation
511 << "\n\t theIllumAzimuth = " << theIllumAzimuth
512 << endl;
513 }
514
515 if (traceDebug()) CLOG << "returning..." << endl;
516 }
517
518
519 //******************************************************************************
520 // This method sets the rigorous model adjustable parameters to
521 // their initial values.
522 //******************************************************************************
initAdjParms()523 void ossimRS1SarModel::initAdjParms()
524 {
525 // Adjustable model not yet implemented
526 #if 0
527
528 static const char* MODULE = "ossimRS1SarModel::initAdjParms()";
529 if (traceDebug()) CLOG << "entering..." << endl;
530
531 theInTrackOffset = (init_adj_parm[INTRACK_OFFSET]+adj_parm[INTRACK_OFFSET])
532 * adj_sigma[INTRACK_OFFSET];
533
534 theCrTrackOffset = (init_adj_parm[CRTRACK_OFFSET]+adj_parm[CRTRACK_OFFSET])
535 * adj_sigma[CRTRACK_OFFSET];
536
537 theRadialOffset = (init_adj_parm[RADIAL_OFFSET]+adj_parm[RADIAL_OFFSET])
538 * adj_sigma[RADIAL_OFFSET];
539
540 theLineScale = (init_adj_parm[LINE_SCALE]+adj_parm[LINE_SCALE])
541 * adj_sigma[LINE_SCALE];
542
543 theSkew = (init_adj_parm[SKEW]+adj_parm[SKEW])
544 * adj_sigma[SKEW];
545
546 theOrientation = (init_adj_parm[ORIENTATION]+adj_parm[ORIENTATION]);
547
548 // Initialize base-class initial adjustable parameter array:
549 init_adj_parm[INTRACK_OFFSET] = theInTrackOffset/adj_sigma[INTRACK_OFFSET];
550 init_adj_parm[CRTRACK_OFFSET] = theCrTrackOffset/adj_sigma[CRTRACK_OFFSET];
551 init_adj_parm[RADIAL_OFFSET] = theRadialOffset/adj_sigma[RADIAL_OFFSET];
552 init_adj_parm[LINE_SCALE] = theLineScale/adj_sigma[LINE_SCALE];
553 init_adj_parm[SKEW] = theSkew/adj_sigma[SKEW];
554 init_adj_parm[ORIENTATION] = theOrientation/adj_sigma[ORIENTATION];
555
556 // Initialize sensor adjustable parameter description strings:
557 strncpy (adj_desc[INTRACK_OFFSET], "intrack_offset", MAX_DESC_CHARS);
558 strncpy (adj_desc[CRTRACK_OFFSET], "crtrack_offset", MAX_DESC_CHARS);
559 strncpy (adj_desc[RADIAL_OFFSET], "radial_offset", MAX_DESC_CHARS);
560 // strncpy (adj_desc[SAMP_SCALE], "samp_scale", MAX_DESC_CHARS);
561 strncpy (adj_desc[LINE_SCALE], "line_scale", MAX_DESC_CHARS);
562 strncpy (adj_desc[SKEW], "image_skew", MAX_DESC_CHARS);
563 strncpy (adj_desc[ORIENTATION], "image_orientation",MAX_DESC_CHARS);
564
565 // Initialize the adj_parms (parameter corrections) to 0:
566 for (int i=0; i<num_adj_parm; i++)
567 adj_parm[i] = 0.0;
568
569 if (traceDebug()) CLOG << "returning..." << endl;
570 #endif
571 }
572
573
574 //******************************************************************************
575 // This method parses the SAR header information associated with the
576 // ephemeris state vectors, and instantiates a Lagrange interpolator object
577 // for interpolating intermediate ehemeris states.
578 //******************************************************************************
establishEphemeris()579 void ossimRS1SarModel::establishEphemeris()
580 {
581 static const char MODULE[] = "ossimRS1SarModel::establishEphemeris()";
582 if (traceDebug()) CLOG << "entering..." << endl;
583
584 char buf[23];
585
586 // Establish pointer to platform position data record:
587 const pos_data_rec* pdr = theCeosData->posDataRec();
588
589 // Establish the Greenwich mean hour angle at time of first sample:
590 char gha_str [23];
591 strncpy(gha_str, pdr->hr_angle, 22);
592 gha_str[22] = '\0';
593 theGHA = atof(gha_str);
594
595 // Fetch number of points in record and allocate memory:
596 strncpy(buf, pdr->ndata, 4);
597 buf[4] = '\0';
598 int numPoints = atoi(buf);
599
600 // Fetch sampling period:
601 strncpy(buf, pdr->data_int, 22);
602 buf[22] = '\0';
603 double sampling_period = atof(buf);
604
605 // Establish the first sample time's day:
606 strncpy(buf, pdr->gmt_day, 4);
607 buf[4] = '\0';
608 int first_day = atoi(buf);
609 double time_offset = (double) (first_day - theFirstLineDay)*SEC_PER_DAY;
610
611 // Establish the first sample time in seconds from beginning of day:
612 strncpy(buf, pdr->gmt_sec, 22);
613 buf[22] = '\0';
614 theEphFirstSampTime = atof(buf) + time_offset;
615
616 int i, j;
617 double sample_time = theEphFirstSampTime;
618 buf[22] = '\0';
619 NEWMAT::ColumnVector eciArpPos(3);
620 NEWMAT::ColumnVector eciArpVel(3);
621 NEWMAT::ColumnVector ecfArpPos(3);
622 NEWMAT::ColumnVector ecfArpVel(3);
623 NEWMAT::ColumnVector earthVel(3);
624 NEWMAT::Matrix xform;
625 theArpPosInterp = new ossimLagrangeInterpolator;
626 theArpVelInterp = new ossimLagrangeInterpolator;
627
628 // Loop over each point, converting them to vectors:
629 for (i=0; i<numPoints; i++)
630 {
631 for (j=0; j<3; j++)
632 {
633 strncpy(buf, pdr->pos_vect[i].pos[j], 22);
634 eciArpPos[j] = atof(buf);
635 strncpy(buf, pdr->pos_vect[i].vel[j], 22);
636 eciArpVel[j] = atof(buf)/1000.0;
637 }
638
639 // Convert pos and vel vectors from ECI to ECF:
640 eciToEcfXform(sample_time, xform);
641 ecfArpPos = xform*eciArpPos;
642 ecfArpVel = xform*eciArpVel;
643 theArpPosInterp->addData(sample_time, ecfArpPos);
644
645 // Correct velocity for earth rotation:
646 earthVel[0] = -ecfArpPos[1]*EARTH_ANGULAR_VELOCITY;
647 earthVel[1] = ecfArpPos[0]*EARTH_ANGULAR_VELOCITY;
648 earthVel[2] = 0.0;
649 ecfArpVel = ecfArpVel - earthVel;
650 theArpVelInterp->addData(sample_time, ecfArpVel);
651
652 sample_time += sampling_period;
653 }
654
655 if (traceDebug())
656 {
657 CLOG << "DEBUG -- "
658 << "\n\t theGHA: " << theGHA
659 << "\n\t numPoints: " << numPoints
660 << "\n\t sampling_period: " << sampling_period
661 << "\n\t theEphFirstSampTime: " << theEphFirstSampTime << endl;
662 }
663
664 if (traceDebug()) CLOG << "returning..." << endl;
665 }
666
667
668 //******************************************************************************
669 // Returns 3x3 transform to rotate ECI into ECF. The argument is
670 // in seconds of the day.
671 //******************************************************************************
eciToEcfXform(const double & acq_time,NEWMAT::Matrix & xform) const672 void ossimRS1SarModel::eciToEcfXform(const double& acq_time, NEWMAT::Matrix& xform) const
673 {
674 static const char MODULE[] = "ossimRS1SarModel::eciToEcfXform(acq_time)";
675 if (traceDebug()) CLOG << "entering..." << endl;
676
677 // Determine time elapsed since ephemeris first sample time:
678 double delta_t = acq_time - theEphFirstSampTime;
679
680 // Compute GHA at time of this image line (in degrees):
681 double gha = theGHA + delta_t*DEG_PER_SEC;
682 if (gha >= 360.0)
683 gha -=360.0;
684 else if (gha < 0.0)
685 gha += 360.0;
686
687 // Establish rotation:
688 double cos_gha = ossim::cosd(gha);
689 double sin_gha = ossim::sind(gha);
690
691 // Establish rotation matrix:
692 xform = ossimMatrix3x3::create(cos_gha, sin_gha, 0.0,
693 -sin_gha, cos_gha, 0.0,
694 0.0, 0.0, 1.0);
695 if (traceDebug())
696 {
697 CLOG << "DEBUG -- "
698 << "\n\t acq_time = " << acq_time
699 << "\n\t delta_t = " << delta_t
700 << "\n\t gha = " << gha << endl;
701 }
702
703 if (traceDebug()) CLOG << "returning..." << endl;
704 }
705
706
707 //******************************************************************************
708 // PRIVATE METHOD: ossimRS1SarModel::establishOrpInterp()
709 //
710 // This method reads the image data file for obtaining the local ORPs (taken
711 // to be the lat-lon of the nearest pixel for each line), and the zero-doppler
712 // time for each line. These ORPs are used for computing SPNs.
713 //
714 // The processed data records correspond to lines of imagery. Accompanying
715 // each data line is the lat-lon (at the ellipsoid) of the first, middle, and
716 // last sample of the line, the time of zero-doppler (broadside) imaging for
717 // the line.
718 //
719 // A Lagrange interpolator is established with N data points representing the
720 // local ORPs and their corresponding image line numbers
721 // uniformly distributed over the image.
722 //
723 // If the image is a ground range product, need to also read the ground-to-
724 // slant range conversion coefficients.
725 //
726 //******************************************************************************
establishOrpInterp()727 void ossimRS1SarModel::establishOrpInterp()
728 {
729 static const char MODULE[] = "ossimRS1SarModel::establishOrpInterp()";
730 if (traceDebug()) CLOG << "entering..." << endl;
731
732 if (traceDebug())
733 CLOG << "DEBUG -- " << endl;
734
735 static const int NUM_DATA_POINTS = 11;
736
737 desc_rec descRec;
738 pdr_prefix_rec prefix;
739 int sizeOfDescRec = sizeof(descRec);
740 int sizeOfPrefRec = sizeof(prefix);
741 int headerSize = sizeOfDescRec + sizeOfPrefRec;
742 char buf[] = "123456";
743 ossimGpt localOrp (0.0, 0.0, theRefHeight);
744 std::vector<double> line_numbers_list;
745 std::vector<NEWMAT::ColumnVector> orpVectorList; // X, Y, Z
746
747 if (traceDebug())
748 {
749 CLOG << "DEBUG:"
750 << "\nsizeOfDescRec: " << sizeOfDescRec
751 << "\nsizeOfPrefRec: " << sizeOfPrefRec << endl;
752 }
753
754 // Open the image data file given the image directory name and seek to first
755 // SAR data record:
756 FILE* fptr = fopen(theCeosData->imageFile().chars(), "r");
757 if (!fptr)
758 {
759 CLOG << "ERROR: Could not open data file <" << theCeosData->imageFile()
760 << ">" << endl;
761 setErrorStatus();
762 return;
763 }
764 fseek(fptr, theCeosData->imopDescRec()->desc.length, SEEK_SET);
765
766 // Declare some constants and variables used in loop:
767 int last_line = (int) theImageSize.line - 1;
768 int delta_line = (int)ceil(last_line/((double)NUM_DATA_POINTS-1.0));
769 int line_number = 0;
770
771 // Compute the number of SAR data bytes to skip to reach prefix of next
772 // record of interest:
773 strncpy(buf, theCeosData->imopDescRec()->l_dataset, 6);
774 int record_size = atoi(buf);
775
776 // Loop over each imaging line:
777 int index;
778 theLocalOrpInterp = new ossimLagrangeInterpolator;
779 for (index=0; index<NUM_DATA_POINTS; index++)
780 {
781 if (traceDebug())
782 clog << "\t Processing ORP for line number " << line_number << endl;
783
784 // Read the image line prefix info:
785 fread(&descRec, sizeOfDescRec, 1, fptr);
786 fread(&prefix, sizeOfPrefRec, 1, fptr);
787
788 // Verify correct line number being read:
789 if ((prefix.line_num-1) != line_number)
790 {
791 CLOG << "\n\tERROR: Synchronization error reading image file. "
792 << "Expected line number " << line_number+1
793 << " but read line number " << prefix.line_num << "." << endl;
794 setErrorStatus();
795 fclose(fptr);
796 fptr = 0;
797 return;
798 }
799
800 // Latch the imaging time for the first line collected:
801 if (index == 0)
802 {
803 theFirstLineDay = prefix.acq_date.acq_day;
804 theFirstLineTime = (double) prefix.acq_date.acq_msec/1000.0;
805 }
806
807 // Read the local ORP for this line and convert to ECF:
808 if (theDirectionFlag == ASCENDING)
809 {
810 localOrp.lat = (double) prefix.lat_first / 1.0e6;
811 localOrp.lon = (double) prefix.long_first / 1.0e6;
812 }
813 else
814 {
815 localOrp.lat = (double) prefix.lat_last / 1.0e6;
816 localOrp.lon = (double) prefix.long_last / 1.0e6;
817 }
818
819 theLocalOrpInterp->addData(line_number, ossimEcefPoint(localOrp).toVector());
820
821 // Update the line number for next record, insuring we don't go passed
822 // the last line:
823 line_number += delta_line;
824 if (line_number > last_line)
825 {
826 delta_line -= line_number - last_line; // adjust for fseek below
827 line_number = last_line;
828 }
829
830 // Advance to read the next prefix record of interest:
831 fseek(fptr, (delta_line)*record_size-headerSize, SEEK_CUR);
832 }
833
834 // Compute time interval between each line, handle possible day rollover:
835 double lastLineTime = (double) prefix.acq_date.acq_msec/1000.0;
836 if ((lastLineTime-theFirstLineTime) > (SEC_PER_DAY/2.0))
837 lastLineTime -= SEC_PER_DAY;
838 else if ((theFirstLineTime-lastLineTime) > (SEC_PER_DAY/2.0))
839 lastLineTime += SEC_PER_DAY;
840 theTimePerLine = (lastLineTime - theFirstLineTime)/line_number;
841
842 if (traceDebug())
843 {
844 clog << "\n\t line_number (last) = " << line_number
845 << "\n\t last_line (in image) = " << last_line
846 << "\n\t theFirstLineTime = " << theFirstLineTime
847 << "\n\t lastLineTime = " << lastLineTime
848 << "\n\t theTimePerLine = " << theTimePerLine
849 << endl;
850 }
851
852 if (traceDebug()) CLOG << "returning..." << endl;
853 }
854
855 //******************************************************************************
856 // This method establishes a common mean intrack-crosstrack-radial LSR
857 // space at the ARP.
858 //******************************************************************************
establishVehicleSpace()859 void ossimRS1SarModel::establishVehicleSpace()
860 {
861 static const char MODULE[] = "ossimRS1SarModel::establishVehicleSpace()";
862 if (traceDebug()) CLOG << "entering..." << endl;
863
864 // Obtain the imaging time when the vehicle is at the ARP:
865 double line = theImageSize.line/2.0;
866 double arpTime = theFirstLineTime + line*theTimePerLine;
867
868 // Obtain the ephemeris at this time:
869 NEWMAT::ColumnVector arpPos(3);
870 theArpPosInterp->interpolate(arpTime, arpPos);
871 NEWMAT::ColumnVector arpVel(3);
872 theArpVelInterp->interpolate(arpTime, arpVel);
873
874 // Convert ephemeris to ECF:
875 NEWMAT::Matrix xform;
876 eciToEcfXform(arpTime, xform);
877 ossimEcefPoint ecfArpPos (xform * arpPos);
878 ossimEcefVector ecfArpVel (xform * arpVel);
879
880 // Need to correct the velocity vector by the earth rotational velocity:
881 ossimEcefVector earthVel(-ecfArpPos.y()*EARTH_ANGULAR_VELOCITY,
882 ecfArpPos.x()*EARTH_ANGULAR_VELOCITY,
883 0.0);
884 ecfArpVel = ecfArpVel - earthVel;
885
886 // Establish LSR space for intrack-crosstrack-radial at vehicle:
887 ossimEcefVector intrackDir (ecfArpVel);
888 ossimEcefVector crtrackDir (ecfArpPos.data().cross(intrackDir.data()));
889 theVehicleSpace = ossimLsrSpace (ecfArpPos, intrackDir, crtrackDir, 0);
890
891 if (traceDebug())
892 {
893 CLOG << "DEBUG -- "
894 << "\n\t theVehicleSpace: " << theVehicleSpace << endl;
895 }
896
897 if (traceDebug()) CLOG << "returning..." << endl;
898 }
899
900
901 //******************************************************************************
902 // For scan-mode imagery, this method interpolates the lat/lon given the image pixel.
903 //******************************************************************************
interpolatedScanORP(const ossimDpt & orp,ossimEcefPoint & orp_ecf) const904 void ossimRS1SarModel::interpolatedScanORP(const ossimDpt& orp, ossimEcefPoint& orp_ecf) const
905 {
906 static const char MODULE[] = "ossimRS1SarModel::interpolatedScanORP(gDblPoint)";
907 if (traceDebug()) CLOG << "entering..." << endl;
908
909 ossimGpt gpt;
910 gpt.lat = theLatGrid(orp);
911 gpt.lon = theLonGrid(orp);
912
913 // Convert to ECF:
914 orp_ecf = ossimEcefPoint(gpt);
915
916 if (traceDebug()) CLOG << "returning..." << endl;
917 }
918
919
920 //******************************************************************************
921 // PRIVATE METHOD: ossimRS1SarModel::establishOrpGrid()
922 //
923 // This method establishes a grid of lat-lon across the image for interpolating
924 // the geographic point on the ellipsoid given a pixel value. This is used for
925 // scan-mode imagery for obtaining a local ORP.
926 //
927 //******************************************************************************
establishOrpGrid()928 void ossimRS1SarModel::establishOrpGrid()
929 {
930 static const char MODULE[] = "ossimRS1SarModel::establishOrpGrid()";
931 if (traceDebug()) CLOG << "entering..." << endl;
932
933 if (traceDebug())
934 {
935 CLOG << "DEBUG -- " << endl;
936 }
937
938 static const int NUM_GRID_POINTS_U = 11;
939
940 desc_rec descRec;
941 pdr_prefix_rec prefix;
942 int sizeOfDescRec = sizeof(descRec);
943 int sizeOfPrefRec = sizeof(prefix);
944 int headerSize = sizeOfDescRec + sizeOfPrefRec;
945 char buf[] = "123456";
946 ossimDpt gridSize ((double) NUM_GRID_POINTS_U, 3.0);
947 ossimDpt cellSize (theImageSize.line/gridSize.u, theImageSize.samp/gridSize.v);
948
949 // Because we are doing integer arithmetic, it will be necessary to
950 // consider the last grid cells as being of a slightly different size:
951 //theLastGridCellSize.line = theImageSize.line - 1.0 -
952 // (theGridSize.u-2.0)*theGridCellSize.u;
953 //theLastGridCellSize.samp = theImageSize.samp - 1.0 -
954 // (theGridSize.v-2.0)*theGridCellSize.v;
955
956 // Open the image data file given the image directory name and seek to first
957 // SAR data record:
958 FILE* fptr = fopen(theCeosData->imageFile().chars(), "r");
959 if (!fptr)
960 {
961 CLOG << "ERROR: Could not open data file <" << theCeosData->imageFile()
962 << ">" << endl;
963 setErrorStatus();
964 return;
965 }
966 fseek(fptr, theCeosData->imopDescRec()->desc.length, SEEK_SET);
967
968 // Declare some constants and variables used in loop:
969 int last_line = (int) theImageSize.line - 1;
970 int delta_line = (int) cellSize.u;
971 int line_number = 0;
972
973 // Allocate memory for coarse grid:
974 theLatGrid.initialize(gridSize, ossimDpt(0,0), ossimDpt(1,1));
975 theLonGrid.initialize(gridSize, ossimDpt(0,0), ossimDpt(1,1));
976
977 // Compute the number of SAR data bytes to skip to reach prefix of next
978 // record of interest:
979 strncpy(buf, theCeosData->imopDescRec()->l_dataset, 6);
980 int record_size = atoi(buf);
981 int test_line = (int) gridSize.u - 2;
982
983 // Loop over each imaging line:
984 for (int u=0; u<(int)gridSize.u; u++)
985 {
986 if (traceDebug())
987 {
988 clog << "\t Processing grid line number " << u << endl;
989 }
990
991 // Read the image line prefix info:
992 fread(&descRec, sizeOfDescRec, 1, fptr);
993 fread(&prefix, sizeOfPrefRec, 1, fptr);
994
995 // Verify correct line number being read:
996 if ((prefix.line_num-1) != line_number)
997 {
998 CLOG << "\n\tERROR: Synchronization error reading image file. "
999 << "Expected line number " << line_number+1
1000 << " but read line number " << prefix.line_num << "." << endl;
1001 setErrorStatus();
1002 fclose(fptr);
1003 fptr = 0;
1004 return;
1005 }
1006
1007 // Latch the imaging time for the first line collected:
1008 if (u == 0)
1009 {
1010 theFirstLineDay = prefix.acq_date.acq_day;
1011 theFirstLineTime = (double) prefix.acq_date.acq_msec/1000.0;
1012 }
1013
1014 // Read the ground points for this line:
1015 theLatGrid.setNode(u, 0, (double) prefix.lat_first /1.0e6);
1016 theLatGrid.setNode(u, 1, (double) prefix.lat_mid /1.0e6);
1017 theLatGrid.setNode(u, 2, (double) prefix.lat_last /1.0e6);
1018 theLonGrid.setNode(u, 0, (double) prefix.long_first/1.0e6);
1019 theLonGrid.setNode(u, 1, (double) prefix.long_mid /1.0e6);
1020 theLonGrid.setNode(u, 2, (double) prefix.long_last /1.0e6);
1021
1022 // Update the line number for next record, insuring we don't go passed
1023 // the last line:
1024 if (u == test_line)
1025 {
1026 delta_line = last_line - line_number;
1027 line_number = last_line;
1028 }
1029 else
1030 {
1031 line_number += delta_line;
1032 }
1033
1034 //***
1035 // Advance to read the next prefix record of interest:
1036 //***
1037 fseek(fptr, (delta_line)*record_size-headerSize, SEEK_CUR);
1038 }
1039
1040 // Compute time interval between each line, handle possible day rollover:
1041 double lastLineTime = (double) prefix.acq_date.acq_msec/1000.0;
1042 if ((lastLineTime-theFirstLineTime) > (SEC_PER_DAY/2.0))
1043 {
1044 lastLineTime -= SEC_PER_DAY;
1045 }
1046 else if ((theFirstLineTime-lastLineTime) > (SEC_PER_DAY/2.0))
1047 {
1048 lastLineTime += SEC_PER_DAY;
1049 }
1050 theTimePerLine = (lastLineTime - theFirstLineTime)/line_number;
1051
1052 if (traceDebug())
1053 {
1054 clog << "\n\t line_number (last) = " << line_number
1055 << "\n\t last_line (in image) = " << last_line
1056 << "\n\t theFirstLineTime = " << theFirstLineTime
1057 << "\n\t lastLineTime = " << lastLineTime
1058 << "\n\t theTimePerLine = " << theTimePerLine
1059 << endl;
1060 }
1061
1062 if (traceDebug()) CLOG << "returning..." << endl;
1063 }
1064
1065 //******************************************************************************
1066 // PRIVATE METHOD: deallocateMemory()
1067 //
1068 // This method permits the partial destruction of the object so that it may
1069 // be reconstructed without invoking the destructor.
1070 //
1071 //******************************************************************************
deallocateMemory()1072 void ossimRS1SarModel::deallocateMemory()
1073 {
1074 static const char MODULE[] = "ossimRS1SarModel::deallocateMemory()";
1075 if (traceDebug()) CLOG << "entering..." << endl;
1076
1077 theArpPosInterp = 0;
1078 theArpVelInterp = 0;
1079 theLocalOrpInterp = 0;
1080 theCeosData = 0;
1081 theLatGrid.clear();
1082 theLonGrid.clear();
1083
1084 if (traceDebug()) CLOG << "returning..." << endl;
1085 }
1086
1087
1088