File failed to load: https://isis.astrogeology.usgs.gov/9.0.0/Object/assets/jax/output/NativeMML/config.js
Isis 3 Programmer Reference
ProjectionFactory.cpp
1
5
6/* SPDX-License-Identifier: CC0-1.0 */
7#include "ProjectionFactory.h"
8
9
10#include <cfloat>
11#include <cmath>
12#include <iomanip>
13
14#include "Camera.h"
15#include "Cube.h"
16#include "Displacement.h"
17#include "Distance.h"
18#include "FileName.h"
19#include "IException.h"
20#include "Plugin.h"
21#include "Projection.h"
22#include "RingPlaneProjection.h"
23#include "TProjection.h"
24
25using namespace std;
26
27namespace Isis {
28 Plugin ProjectionFactory::m_projPlugin;
29
52 bool allowDefaults) {
53 // Try to load a plugin file in the current working directory and then
54 // load the system file
55 Plugin p;
56
57 if (m_projPlugin.fileName() == "") {
58 FileName localFile("Projection.plugin");
59 if (localFile.fileExists())
60 m_projPlugin.read(localFile.expanded());
61
62 FileName systemFile("$ISISROOT/lib/Projection.plugin");
63 if (systemFile.fileExists())
64 m_projPlugin.read(systemFile.expanded());
65 }
66
67 try {
68 // Look for info in the mapping group
69 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
70 QString proj = mapGroup["ProjectionName"];
71
72 // Now get the plugin for the projection
73 QFunctionPointer ptr;
74 try {
75 ptr = m_projPlugin.GetPlugin(proj);
76 }
77 catch(IException &e) {
78 QString msg = "Unsupported projection, unable to find plugin for [" +
79 proj + "]";
80 throw IException(e, IException::Unknown, msg, _FILEINFO_);
81 }
82
83 // Now cast that pointer in the proper way
84 Isis::TProjection * (*plugin)(Isis::Pvl & label, bool flag);
85 // plugin = (Isis::TProjection * ( *)(Isis::Pvl & label, bool flag)) ptr;
86 plugin = (Isis::TProjection * ( *)(Isis::Pvl & label, bool flag)) ptr;
87 // Create the projection as requested
88 return (Isis::Projection *) (*plugin)(label, allowDefaults);
89 }
90 catch(IException &e) {
91 QString message = "Unable to initialize Projection information ";
92 message += "from group [Mapping]";
93 throw IException(e, IException::Io, message, _FILEINFO_);
94 }
95 }
96
97
120 bool allowDefaults) {
121 // Try to load a plugin file in the current working directory and then
122 // load the system file
123 Plugin p;
124
125 if (m_projPlugin.fileName() == "") {
126 FileName localFile("Projection.plugin");
127 if (localFile.fileExists())
128 m_projPlugin.read(localFile.expanded());
129
130 FileName systemFile("$ISISROOT/lib/Projection.plugin");
131 if (systemFile.fileExists())
132 m_projPlugin.read(systemFile.expanded());
133 }
134
135 try {
136 // Look for info in the mapping group
137 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
138 QString proj = mapGroup["ProjectionName"];
139
140 // Now get the plugin for the projection
141 QFunctionPointer ptr;
142 try {
143 ptr = m_projPlugin.GetPlugin(proj);
144 }
145 catch(IException &e) {
146 QString msg = "Unsupported projection, unable to find plugin for [" +
147 proj + "]";
148 throw IException(e, IException::Unknown, msg, _FILEINFO_);
149 }
150
151 // Now cast that pointer in the proper way
152 Isis::RingPlaneProjection * (*plugin)(Isis::Pvl & label, bool flag);
153 plugin = (Isis::RingPlaneProjection * ( *)(Isis::Pvl & label, bool flag)) ptr;
154 // Create the projection as requested
155 return (Projection *) (*plugin)(label, allowDefaults);
156 }
157 catch(IException &e) {
158 QString message = "Unable to initialize Projection information ";
159 message += "from group [Mapping]";
160 throw IException(e, IException::Io, message, _FILEINFO_);
161 }
162 }
163
164
189 int &samples, int &lines,
190 bool sizeMatch) {
191 // Create a temporary projection and get the radius at the special latitude
192 Isis::TProjection *proj = (Isis::TProjection *) Create(label, true);
193 double trueScaleLat = proj->TrueScaleLatitude();
194 double localRadius = proj->LocalRadius(trueScaleLat);
195 delete proj;
196
197 IException errors;
198 try {
199 // Try to get the pixel resolution and then compute the scale
200 double scale, pixelResolution;
201 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
202 try {
203 pixelResolution = mapGroup["PixelResolution"];
204 scale = (2.0 * Isis::PI * localRadius) / (360.0 * pixelResolution);
205 }
206
207 // If not get the scale and then compute the pixel resolution
208 catch(IException &e) {
209 errors.append(e);
210
211 scale = mapGroup["Scale"];
212 pixelResolution = (2.0 * Isis::PI * localRadius) / (360.0 * scale);
213 }
214 // Write out the scale and resolution with units and truescale latitude
215 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution", toString(pixelResolution), "meters/pixel"),
216 Isis::Pvl::Replace);
217 mapGroup.addKeyword(Isis::PvlKeyword("Scale", toString(scale), "pixels/degree"), Isis::Pvl::Replace);
218 //mapGroup.addKeyword(Isis::PvlKeyword ("TrueScaleLatitude", trueScaleLat),
219 // Isis::Pvl::Replace);
220
221 // Get the upper left corner from the labels if possible
222 // This forces an exact match of projection parameters for
223 // output cubes
224 bool sizeFound = false;
225 double upperLeftX = Null, upperLeftY = Null;
226 if (label.hasObject("IsisCube")) {
227 Isis::PvlGroup &dims = label.findGroup("Dimensions", Isis::Pvl::Traverse);
228 samples = dims["Samples"];
229 lines = dims["Lines"];
230
231 upperLeftX = mapGroup["UpperLeftCornerX"];
232 upperLeftY = mapGroup["UpperLeftCornerY"];
233 sizeFound = true;
234 }
235 if (!sizeMatch) sizeFound = false;
236
237 // Initialize the rest of the projection
238 proj = (Isis::TProjection *) Create(label, true);
239
240 // Couldn't find the cube size from the labels so compute it
241 if (!sizeFound) {
242 if (!proj->HasGroundRange()) {
243 QString msg = "Invalid ground range [MinimumLatitude,MaximumLatitude,";
244 msg += "MinimumLongitude,MaximumLongitude] missing or invalid";
245 throw IException(IException::Unknown, msg, _FILEINFO_);
246 }
247
248 double minX, maxX, minY, maxY;
249 if (!proj->XYRange(minX, maxX, minY, maxY)) {
250 QString msg = "Invalid ground range [MinimumLatitude,MaximumLatitude,";
251 msg += "MinimumLongitude,MaximumLongitude] cause invalid computation ";
252 msg += "of image size";
253 throw IException(IException::Unknown, msg, _FILEINFO_);
254 }
255
256
257 // Convert upperleft coordinate to units of pixel
258 // Truncate it to the nearest whole pixel (floor/ceil)
259 // Convert it back to meters. But don't do this if
260 // the X/Y position is already close to a whole pixel because
261 // the floor/ceil function could cause an extra pixel to be added
262 // just due to machine precision issues
263
264 bool flipX = false;
265 bool flipY = false;
266
267 double minXFlipped = minX;
268 double maxXFlipped = maxX;
269 double minYFlipped = minY;
270 double maxYFlipped = maxY;
271
272 //New range is (-1)*[0,maxX] = [-maxX, 0]
273
274 if (minX == 0) {
275 minXFlipped = -maxX;
276 maxXFlipped = 0;
277 flipX = true;
278
279 }
280
281 if (fabs(fmod(minXFlipped, pixelResolution)) > 1.0e-6) {
282 if (pixelResolution - fabs(fmod(minXFlipped, pixelResolution)) > 1.0e-6) {
283 double sampleOffset = floor(minXFlipped / pixelResolution);
284
285 minXFlipped = sampleOffset * pixelResolution;
286
287 }
288 }
289
290
291 // make sure that the distance from minX to maxX is at least one pixel wide
292 // so we have at least one sample in the created cube
293
294
295 if (maxXFlipped < minXFlipped + pixelResolution) {
296 maxXFlipped = minXFlipped + pixelResolution;
297 }
298
299
300 //New range is (-1)*[minY,0] = [0,-minY]
301 if (maxY == 0) {
302 maxYFlipped = -minY;
303 minYFlipped = 0;
304 flipY = true;
305
306 }
307
308 if (fabs(fmod(maxYFlipped, pixelResolution)) > 1.0e-6) {
309 if (abs(pixelResolution - fabs(fmod(maxYFlipped, pixelResolution))) > 1.0e-6) {
310 double lineOffset = ceil(maxYFlipped / pixelResolution);
311 maxYFlipped = lineOffset * pixelResolution;
312
313 }
314 }
315
316
317 // make sure that the distance from minY to maxY is at least one pixel wide
318 // so we have at least one line in the created cube
319 if (minYFlipped > maxYFlipped - pixelResolution) {
320 minYFlipped = maxYFlipped - pixelResolution;
321 }
322
323 // Determine the number of samples and lines
324
325 samples = (int)((maxXFlipped - minXFlipped) / pixelResolution + 0.5);
326 lines = (int)((maxYFlipped - minYFlipped) / pixelResolution + 0.5);
327
328
329 // Set the upper left corner and add to the labels
330 if (flipX) {
331 upperLeftX = 0;
332
333 }
334 else {
335 upperLeftX = minXFlipped;
336 }
337 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX", toString(upperLeftX)),
338 Isis::Pvl::Replace);
339
340 if (flipY) {
341 upperLeftY = 0;
342
343 }
344 else {
345 upperLeftY = maxYFlipped;
346
347 }
348
349 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY", toString(upperLeftY)),
350 Isis::Pvl::Replace);
351
352 // Write it in pixel units as well
353#if 0
354 lineOffset += 0.5; // This matches the PDS definition
355 sampleOffset += 0.5; // of the offsets (center of pixel). This statement is questionable!
356 mapGroup.addKeyword(Isis::PvlKeyword("LineProjectionOffset", lineOffset),
357 Isis::Pvl::Replace);
358 mapGroup.addKeyword(Isis::PvlKeyword("SampleProjectionOffset", sampleOffset),
359 Isis::Pvl::Replace);
360#endif
361 }
362
363
364 // Make sure labels have good units
365 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution",
366 (QString) mapGroup["PixelResolution"],
367 "meters/pixel"), Isis::Pvl::Replace);
368
369 mapGroup.addKeyword(Isis::PvlKeyword("Scale",
370 (QString) mapGroup["Scale"],
371 "pixels/degree"), Isis::Pvl::Replace);
372
373 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX",
374 (QString) mapGroup["UpperLeftCornerX"],
375 "meters"), Isis::Pvl::Replace);
376
377 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY",
378 (QString) mapGroup["UpperLeftCornerY"],
379 "meters"), Isis::Pvl::Replace);
380
381 mapGroup.addKeyword(Isis::PvlKeyword("EquatorialRadius",
382 (QString) mapGroup["EquatorialRadius"],
383 "meters"), Isis::Pvl::Replace);
384
385 mapGroup.addKeyword(Isis::PvlKeyword("PolarRadius",
386 (QString) mapGroup["PolarRadius"],
387 "meters"), Isis::Pvl::Replace);
388
389 // Add the mapper from pixel coordinates to projection coordinates
390 PFPixelMapper *pixelMapper = new PFPixelMapper(pixelResolution, upperLeftX, upperLeftY);
391 proj->SetWorldMapper(pixelMapper);
392
395 }
396 catch(IException &e) {
397 QString msg = "Unable to create projection";
398 if (label.fileName() != "") msg += " from file [" + label.fileName() + "]";
399 IException finalError(IException::Unknown, msg, _FILEINFO_);
400 finalError.append(errors);
401 finalError.append(e);
402 throw finalError;
403 }
404 return (Isis::Projection *) proj;
405 }
406
407
433 int &samples, int &lines, bool sizeMatch) {
434
435 // Create a temporary projection and get the radius at the special radius
436 // NOTE: by "special radius" we mean that radius where the projection is
437 // not distorted
439 double localRadius = proj->TrueScaleRingRadius();
440 delete proj;
441
442 IException errors;
443 try {
444 // Try to get the pixel resolution and then compute the scale
445 double scale, pixelResolution;
446 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
447 try {
448 pixelResolution = mapGroup["PixelResolution"];
449 scale = (2.0 * Isis::PI * localRadius) / (360.0 * pixelResolution);
450 }
451 catch(IException &e) {
452 // If this fails, use the scale to compute the pixel resolution
453 errors.append(e);
454
455 scale = mapGroup["Scale"];
456 pixelResolution = (2.0 * Isis::PI * localRadius) / (360.0 * scale);
457 }
458
459 // Write out the scale and resolution with units and truescale radius
460 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution", toString(pixelResolution),
461 "meters/pixel"),
462 Isis::Pvl::Replace);
463 mapGroup.addKeyword(Isis::PvlKeyword("Scale", toString(scale), "pixels/degree"),
464 Isis::Pvl::Replace);
465
466 //mapGroup.AddKeyword(Isis::PvlKeyword ("TrueScaleRadius", trueScaleRadius),
467 // Isis::Pvl::Replace);
468
469 // Get the upper left corner from the labels if possible
470 // This forces an exact match of projection parameters for
471 // output cubes
472 bool sizeFound = false;
473 double upperLeftX = Null, upperLeftY = Null;
474 if (label.hasObject("IsisCube")) {
475 Isis::PvlGroup &dims = label.findGroup("Dimensions", Isis::Pvl::Traverse);
476 samples = dims["Samples"];
477 lines = dims["Lines"];
478
479 upperLeftX = mapGroup["UpperLeftCornerX"];
480 upperLeftY = mapGroup["UpperLeftCornerY"];
481 sizeFound = true;
482 }
483
484 if (!sizeMatch) {
485 sizeFound = false;
486 }
487
488 // Initialize the rest of the projection
489 proj = (Isis::RingPlaneProjection *) RingsCreate(label, true);
490
491 // Couldn't find the cube size from the labels so compute it
492 if (!sizeFound) {
493 if (!proj->HasGroundRange()) {
494 QString msg = "Invalid ring range [MinimumRingRadius,MaximumRingRadius,";
495 msg += "MinimumRingLongitude,MaximumRingLongitude] missing or invalid";
496 throw IException(IException::Unknown, msg, _FILEINFO_);
497 }
498
499 double minX, maxX, minY, maxY;
500 if (!proj->XYRange(minX, maxX, minY, maxY)) {
501 QString msg = "Invalid ring range [MinimumRingRadius,MaximumRingRadius,";
502 msg += "MinimumRingLongitude,MaximumRingLongitude] cause invalid computation ";
503 msg += "of image size";
504 throw IException(IException::Unknown, msg, _FILEINFO_);
505 }
506
507 // Convert upperleft coordinate to pixel units
508 // Truncate it to the nearest whole pixel (floor/ceil)
509 // Convert it back to meters. But don't do this if
510 // the X/Y position is already close to a whole pixel because
511 // the floor/ceil function could cause an extra pixel to be added
512 // just due to machine precision issues
513 if (fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
514 if (pixelResolution - fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
515 double sampleOffset = floor(minX / pixelResolution);
516 minX = sampleOffset * pixelResolution;
517 }
518 }
519 // make sure that the distance from minX to maxX is at least one pixel wide
520 // so we have at least one sample in the created cube
521 if (maxX < minX + pixelResolution) {
522 maxX = minX + pixelResolution;
523 }
524 if (fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
525 if (pixelResolution - fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
526 double lineOffset = -1.0 * ceil(maxY / pixelResolution);
527 maxY = -1.0 * lineOffset * pixelResolution;
528 }
529 }
530 // make sure that the distance from minY to maxY is at least one pixel wide
531 // so we have at least one line in the created cube
532 if (minY > maxY - pixelResolution) {
533 minY = maxY - pixelResolution;
534 }
535
536 // Determine the number of samples and lines
537 samples = (int)((maxX - minX) / pixelResolution + 0.5);
538 lines = (int)((maxY - minY) / pixelResolution + 0.5);
539
540 // Set the upper left corner and add to the labels
541 upperLeftX = minX;
542 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX", toString(upperLeftX)),
543 Isis::Pvl::Replace);
544
545 upperLeftY = maxY;
546 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY", toString(upperLeftY)),
547 Isis::Pvl::Replace);
548
549 // Write it in pixel units as well
550#if 0
551 lineOffset += 0.5; // This matches the PDS definition
552 sampleOffset += 0.5; // of the offsets (center of pixel). This statement is questionable!
553 mapGroup.AddKeyword(Isis::PvlKeyword("LineProjectionOffset", lineOffset),
554 Isis::Pvl::Replace);
555 mapGroup.AddKeyword(Isis::PvlKeyword("SampleProjectionOffset", sampleOffset),
556 Isis::Pvl::Replace);
557#endif
558 }
559
560
561 // Make sure labels have good units
562 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution",
563 (QString) mapGroup["PixelResolution"],
564 "meters/pixel"), Isis::Pvl::Replace);
565
566 mapGroup.addKeyword(Isis::PvlKeyword("Scale",
567 (QString) mapGroup["Scale"],
568 "pixels/degree"), Isis::Pvl::Replace);
569
570 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX",
571 (QString) mapGroup["UpperLeftCornerX"],
572 "meters"), Isis::Pvl::Replace);
573
574 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY",
575 (QString) mapGroup["UpperLeftCornerY"],
576 "meters"), Isis::Pvl::Replace);
577
578 // Add the mapper from pixel coordinates to projection coordinates
579 PFPixelMapper *pixelMapper = new PFPixelMapper(pixelResolution, upperLeftX, upperLeftY);
580 proj->SetWorldMapper(pixelMapper);
581
584 }
585 catch(IException &e) {
586 QString msg = "Unable to create projection";
587 if (label.fileName() != "") msg += " from file [" + label.fileName() + "]";
588 IException finalError(IException::Unknown, msg, _FILEINFO_);
589 finalError.append(errors);
590 finalError.append(e);
591 throw finalError;
592 }
593 return (Isis::Projection *) proj;
594 }
595
596
622 int &samples, int &lines,
623 Camera &cam) {
624 // Create a temporary projection and get the radius at the special latitude
625 Isis::TProjection *proj = (Isis::TProjection *) Create(label, true);
626 double trueScaleLat = proj->TrueScaleLatitude();
627 double localRadius = proj->LocalRadius(trueScaleLat);
628 delete proj;
629
630 try {
631 // Try to get the pixel resolution and then compute the scale
632 double scale, pixelResolution;
633 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
634 try {
635 pixelResolution = mapGroup["PixelResolution"];
636 scale = (2.0 * Isis::PI * localRadius) / (360.0 * pixelResolution);
637 }
638
639 // If not get the scale and then compute the pixel resolution
640 catch(IException &) {
641 scale = mapGroup["Scale"];
642 pixelResolution = (2.0 * Isis::PI * localRadius) / (360.0 * scale);
643 }
644 // Write out the scale and resolution with units and truescale latitude
645 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution", toString(pixelResolution), "meters/pixel"),
646 Isis::Pvl::Replace);
647 mapGroup.addKeyword(Isis::PvlKeyword("Scale", toString(scale), "pixels/degree"), Isis::Pvl::Replace);
648 //mapGroup.AddKeyword(Isis::PvlKeyword ("TrueScaleLatitude", trueScaleLatitude),
649 // Isis::Pvl::Replace);
650
651 // Initialize the rest of the projection
652 proj = (Isis::TProjection *) Create(label, true);
653 double minX = DBL_MAX;
654 double maxX = -DBL_MAX;
655 double minY = DBL_MAX;
656 double maxY = -DBL_MAX;
657
658 // Walk the boundaries of the camera to determine the x/y range
659 int eband = cam.Bands();
660 if (cam.IsBandIndependent()) eband = 1;
661 for(int band = 1; band <= eband; band++) {
662 cam.SetBand(band);
663
664 // Loop for each line testing the left and right sides of the image
665 for(int line = 0; line <= cam.Lines(); line++) {
666 // Look for the first good lat/lon on the left edge of the image
667 // If it is the first or last line then test the whole line
668 int samp;
669 for(samp = 0; samp <= cam.Samples(); samp++) {
670 if (cam.SetImage((double)samp + 0.5, (double)line + 0.5)) {
671 double lat = cam.UniversalLatitude();
672 double lon = cam.UniversalLongitude();
673 proj->SetUniversalGround(lat, lon);
674 if (proj->IsGood()) {
675 if (proj->XCoord() < minX) minX = proj->XCoord();
676 if (proj->XCoord() > maxX) maxX = proj->XCoord();
677 if (proj->YCoord() < minY) minY = proj->YCoord();
678 if (proj->YCoord() > maxY) maxY = proj->YCoord();
679 if ((line != 0) && (line != cam.Lines())) break;
680 }
681 }
682 }
683
684 // Look for the first good lat/lon on the right edge of the image
685 if (samp < cam.Samples()) {
686 for(samp = cam.Samples(); samp >= 0; samp--) {
687 if (cam.SetImage((double)samp + 0.5, (double)line + 0.5)) {
688 double lat = cam.UniversalLatitude();
689 double lon = cam.UniversalLongitude();
690 proj->SetUniversalGround(lat, lon);
691 if (proj->IsGood()) {
692 if (proj->XCoord() < minX) minX = proj->XCoord();
693 if (proj->XCoord() > maxX) maxX = proj->XCoord();
694 if (proj->YCoord() < minY) minY = proj->YCoord();
695 if (proj->YCoord() > maxY) maxY = proj->YCoord();
696 break;
697 }
698 }
699 }
700 }
701 }
702
703 // Special test for ground range to see if either pole is in the image
704 if (cam.SetUniversalGround(90.0, 0.0)) {
705 if (cam.Sample() >= 0.5 && cam.Line() >= 0.5 &&
706 cam.Sample() <= cam.Samples() + 0.5 && cam.Line() <= cam.Lines() + 0.5) {
707 double lat = cam.UniversalLatitude();
708 double lon = cam.UniversalLongitude();
709 proj->SetUniversalGround(lat, lon);
710 if (proj->IsGood()) {
711 if (proj->XCoord() < minX) minX = proj->XCoord();
712 if (proj->XCoord() > maxX) maxX = proj->XCoord();
713 if (proj->YCoord() < minY) minY = proj->YCoord();
714 if (proj->YCoord() > maxY) maxY = proj->YCoord();
715 }
716 }
717 }
718
719 if (cam.SetUniversalGround(-90.0, 0.0)) {
720 if (cam.Sample() >= 0.5 && cam.Line() >= 0.5 &&
721 cam.Sample() <= cam.Samples() + 0.5 && cam.Line() <= cam.Lines() + 0.5) {
722 double lat = cam.UniversalLatitude();
723 double lon = cam.UniversalLongitude();
724 proj->SetUniversalGround(lat, lon);
725 if (proj->IsGood()) {
726 if (proj->XCoord() < minX) minX = proj->XCoord();
727 if (proj->XCoord() > maxX) maxX = proj->XCoord();
728 if (proj->YCoord() < minY) minY = proj->YCoord();
729 if (proj->YCoord() > maxY) maxY = proj->YCoord();
730 }
731 }
732 }
733
734#if 0
735 // Another special test for ground range as we could have the
736 // 0-360 seam running right through the image so
737 // test it as well (the increment may not be fine enough !!!)
738 for(double lat = p_minlat; lat <= p_maxlat; lat += (p_maxlat - p_minlat) / 10.0) {
739 if (SetUniversalGround(lat, 0.0)) {
740 if (Sample() >= 0.5 && Line() >= 0.5 &&
741 Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
742 p_minlon = 0.0;
743 p_maxlon = 360.0;
744 break;
745 }
746 }
747 }
748
749 // Another special test for ground range as we could have the
750 // -180-180 seam running right through the image so
751 // test it as well (the increment may not be fine enough !!!)
752 for(double lat = p_minlat; lat <= p_maxlat; lat += (p_maxlat - p_minlat) / 10.0) {
753 if (SetUniversalGround(lat, 180.0)) {
754 if (Sample() >= 0.5 && Line() >= 0.5 &&
755 Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
756 p_minlon180 = -180.0;
757 p_maxlon180 = 180.0;
758 break;
759 }
760 }
761 }
762#endif
763 }
764
765 // Convert upperleft coordinate to units of pixel
766 // Truncate it to the nearest whole pixel (floor/ceil)
767 // Convert it back to meters. But don't do this if
768 // the X/Y position is already close to a whole pixel because
769 // the floor/ceil function could cause an extra pixel to be added
770 // just due to machine precision issues
771 if (fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
772 if (pixelResolution - fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
773 double sampleOffset = floor(minX / pixelResolution);
774 minX = sampleOffset * pixelResolution;
775 }
776 }
777 // make sure that the distance from minX to maxX is at least one pixel wide
778 // so we have at least one sample in the created cube
779 if (maxX < minX + pixelResolution) {
780 maxX = minX + pixelResolution;
781 }
782 if (fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
783 if (pixelResolution - fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
784 double lineOffset = -1.0 * ceil(maxY / pixelResolution);
785 maxY = -1.0 * lineOffset * pixelResolution;
786 }
787 }
788 // make sure that the distance from minY to maxY is at least one pixel wide
789 // so we have at least one line in the created cube
790 if (minY > maxY - pixelResolution) {
791 minY = maxY - pixelResolution;
792 }
793
794 // Determine the number of samples and lines
795 samples = (int)((maxX - minX) / pixelResolution + 0.5);
796 lines = (int)((maxY - minY) / pixelResolution + 0.5);
797
798 // Set the upper left corner and add to the labels
799 double upperLeftX = minX;
800 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX", toString(upperLeftX)),
801 Isis::Pvl::Replace);
802
803 double upperLeftY = maxY;
804 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY", toString(upperLeftY)),
805 Isis::Pvl::Replace);
806
807 // Make sure labels have good units
808 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution",
809 (QString) mapGroup["PixelResolution"],
810 "meters/pixel"), Isis::Pvl::Replace);
811
812 mapGroup.addKeyword(Isis::PvlKeyword("Scale",
813 (QString) mapGroup["Scale"],
814 "pixels/degree"), Isis::Pvl::Replace);
815
816 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX",
817 (QString) mapGroup["UpperLeftCornerX"],
818 "meters"), Isis::Pvl::Replace);
819
820 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY",
821 (QString) mapGroup["UpperLeftCornerY"],
822 "meters"), Isis::Pvl::Replace);
823
824 mapGroup.addKeyword(Isis::PvlKeyword("EquatorialRadius",
825 (QString) mapGroup["EquatorialRadius"],
826 "meters"), Isis::Pvl::Replace);
827
828 mapGroup.addKeyword(Isis::PvlKeyword("PolarRadius",
829 (QString) mapGroup["PolarRadius"],
830 "meters"), Isis::Pvl::Replace);
831
832 // Add the mapper from pixel coordinates to projection coordinates
833 PFPixelMapper *pixelMapper = new PFPixelMapper(pixelResolution, upperLeftX, upperLeftY);
834 proj->SetWorldMapper(pixelMapper);
835
838 }
839 catch(IException &e) {
840 QString msg = "Unable to create projection";
841 if (label.fileName() != "") msg += " from file [" + label.fileName() + "]";
842 throw IException(e, IException::Unknown, msg, _FILEINFO_);
843 }
844 return (Isis::Projection *) proj;
845 }
846
847
873 int &samples, int &lines, Camera &cam) {
874
875 // Create a temporary projection
877 double localRadius = proj->TrueScaleRingRadius();
878 delete proj;
879
880 IException errors;
881 try {
882 // Try to get the pixel resolution and then compute the scale
883 double scale, pixelResolution;
884 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
885 try {
886 pixelResolution = mapGroup["PixelResolution"];
887 scale = (2.0 * Isis::PI * localRadius) / (360.0 * pixelResolution);
888 // scale = proj->Scale();
889 }
890
891 // If not get the scale and then compute the pixel resolution
892 catch(IException &e) {
893 errors.append(e);
894 scale = mapGroup["Scale"];
895 pixelResolution = (2.0 * Isis::PI * localRadius) / (360.0 * scale);
896 }
897 // Write out the scale and resolution with units and truescale radius
898 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution", toString(pixelResolution),
899 "meters/pixel"), Isis::Pvl::Replace);
900 mapGroup.addKeyword(Isis::PvlKeyword("Scale", toString(scale), "pixels/degree"),
901 Isis::Pvl::Replace);
902
903 // Initialize the rest of the projection
904 proj = (Isis::RingPlaneProjection *) RingsCreate(label, true);
905 double minX = DBL_MAX;
906 double maxX = -DBL_MAX;
907 double minY = DBL_MAX;
908 double maxY = -DBL_MAX;
909
910 // Walk the boundaries of the camera to determine the x/y range
911 int eband = cam.Bands();
912 if (cam.IsBandIndependent()) eband = 1;
913 for(int band = 1; band <= eband; band++) {
914 cam.SetBand(band);
915
916 // Loop for each line testing the left and right sides of the image
917 for(int line = 0; line <= cam.Lines(); line++) {
918 // Look for the first good rad/az on the left edge of the image
919 // If it is the first or last line then test the whole line
920 int samp;
921 for(samp = 0; samp <= cam.Samples(); samp++) {
922 if (cam.SetImage((double)samp + 0.5, (double)line + 0.5)) {
923 double radius = cam.LocalRadius().meters();
924 double az = cam.UniversalLongitude();
925 proj->SetGround(radius, az);
926 if (proj->IsGood()) {
927 if (proj->XCoord() < minX) minX = proj->XCoord();
928 if (proj->XCoord() > maxX) maxX = proj->XCoord();
929 if (proj->YCoord() < minY) minY = proj->YCoord();
930 if (proj->YCoord() > maxY) maxY = proj->YCoord();
931 if ((line != 0) && (line != cam.Lines())) break;
932 }
933 }
934 }
935
936 // Look for the first good rad/az on the right edge of the image
937 if (samp < cam.Samples()) {
938 for(samp = cam.Samples(); samp >= 0; samp--) {
939 if (cam.SetImage((double)samp + 0.5, (double)line + 0.5)) {
940 double radius = cam.LocalRadius().meters();
941 double az = cam.UniversalLongitude();
942 proj->SetGround(radius, az);
943 if (proj->IsGood()) {
944 if (proj->XCoord() < minX) minX = proj->XCoord();
945 if (proj->XCoord() > maxX) maxX = proj->XCoord();
946 if (proj->YCoord() < minY) minY = proj->YCoord();
947 if (proj->YCoord() > maxY) maxY = proj->YCoord();
948 break;
949 }
950 }
951 }
952 }
953 }
954
955#if 0
956 // Another special test for ground range as we could have the
957 // 0-360 seam running right through the image so
958 // test it as well (the increment may not be fine enough !!!)
959 for(double rad = p_minRadius; rad <= p_maxRadius; rad += (p_maxRadius - p_minRadius) / 10.0) {
960 if (SetUniversalGround(rad, 0.0)) {
961 if (Sample() >= 0.5 && Line() >= 0.5 &&
962 Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
963 p_minaz = 0.0;
964 p_maxaz = 360.0;
965 break;
966 }
967 }
968 }
969
970 // Another special test for ground range as we could have the
971 // -180-180 seam running right through the image so
972 // test it as well (the increment may not be fine enough !!!)
973 for(double rad = p_minrad; rad <= p_maxrad; rad += (p_maxrad - p_minrad) / 10.0) {
974 if (SetUniversalGround(rad, 180.0)) {
975 if (Sample() >= 0.5 && Line() >= 0.5 &&
976 Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
977 p_minaz180 = -180.0;
978 p_maxaz180 = 180.0;
979 break;
980 }
981 }
982 }
983#endif
984 }
985
986 // Convert upperleft coordinate to units of pixel
987 // Truncate it to the nearest whole pixel (floor/ceil)
988 // Convert it back to meters. But don't do this if
989 // the X/Y position is already close to a whole pixel because
990 // the floor/ceil function could cause an extra pixel to be added
991 // just due to machine precision issues
992 if (fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
993 if (pixelResolution - fabs(fmod(minX, pixelResolution)) > 1.0e-6) {
994 double sampleOffset = floor(minX / pixelResolution);
995 minX = sampleOffset * pixelResolution;
996 }
997 }
998 // make sure that the distance from minX to maxX is at least one pixel wide
999 // so we have at least one sample in the created cube
1000 if (maxX < minX + pixelResolution) {
1001 maxX = minX + pixelResolution;
1002 }
1003
1004 if (fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
1005 if (pixelResolution - fabs(fmod(maxY, pixelResolution)) > 1.0e-6) {
1006 double lineOffset = -1.0 * ceil(maxY / pixelResolution);
1007 maxY = -1.0 * lineOffset * pixelResolution;
1008 }
1009 }
1010 // make sure that the distance from minY to maxY is at least one pixel wide
1011 // so we have at least one line in the created cube
1012 if (minY > maxY - pixelResolution) {
1013 minY = maxY - pixelResolution;
1014 }
1015
1016 // Determine the number of samples and lines
1017 samples = (int)((maxX - minX) / pixelResolution + 0.5);
1018 lines = (int)((maxY - minY) / pixelResolution + 0.5);
1019
1020 // Set the upper left corner and add to the labels
1021 double upperLeftX = minX;
1022 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX", toString(upperLeftX)),
1023 Isis::Pvl::Replace);
1024
1025 double upperLeftY = maxY;
1026 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY", toString(upperLeftY)),
1027 Isis::Pvl::Replace);
1028
1029 // Make sure labels have good units
1030 mapGroup.addKeyword(Isis::PvlKeyword("PixelResolution",
1031 (QString) mapGroup["PixelResolution"],
1032 "meters/pixel"), Isis::Pvl::Replace);
1033
1034 mapGroup.addKeyword(Isis::PvlKeyword("Scale",
1035 (QString) mapGroup["Scale"],
1036 "pixels/degree"), Isis::Pvl::Replace);
1037
1038 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerX",
1039 (QString) mapGroup["UpperLeftCornerX"],
1040 "meters"), Isis::Pvl::Replace);
1041
1042 mapGroup.addKeyword(Isis::PvlKeyword("UpperLeftCornerY",
1043 (QString) mapGroup["UpperLeftCornerY"],
1044 "meters"), Isis::Pvl::Replace);
1045
1046 // Add the mapper from pixel coordinates to projection coordinates
1047 PFPixelMapper *pixelMapper = new PFPixelMapper(pixelResolution, upperLeftX, upperLeftY);
1048 proj->SetWorldMapper(pixelMapper);
1049
1051 Displacement(upperLeftY, Displacement::Meters));
1052 }
1053 catch(IException &e) {
1054 QString msg = "Unable to create projection";
1055 if (label.fileName() != "") msg += " from file [" + label.fileName() + "]";
1056 throw IException(e, IException::Unknown, msg, _FILEINFO_);
1057 }
1058 return (Isis::Projection *) proj;
1059 }
1060
1061
1072
1073
1084
1085
1096 Isis::TProjection *proj;
1097 try {
1098 // Get the pixel resolution
1099 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
1100 double pixelResolution = mapGroup["PixelResolution"];
1101
1102 // Get the upper left corner
1103 double upperLeftX = mapGroup["UpperLeftCornerX"];
1104 double upperLeftY = mapGroup["UpperLeftCornerY"];
1105
1106 // Initialize the rest of the projection
1107 proj = (Isis::TProjection *) Create(label, true);
1108
1109 // Create a mapper to transform pixels into projection x/y and vice versa
1110 PFPixelMapper *pixelMapper = new PFPixelMapper(pixelResolution, upperLeftX, upperLeftY);
1111 proj->SetWorldMapper(pixelMapper);
1112
1114 Displacement(upperLeftY, Displacement::Meters));
1115 }
1116 catch (IException &e) {
1117 QString msg = "Unable to initialize cube projection";
1118 if (label.fileName() != "") msg += " from file [" + label.fileName() + "]";
1119 throw IException(e, IException::Unknown, msg, _FILEINFO_);
1120 }
1121 return (Isis::Projection *) proj;
1122 }
1123
1124
1136 try {
1137 // Get the pixel resolution
1138 Isis::PvlGroup &mapGroup = label.findGroup("Mapping", Isis::Pvl::Traverse);
1139 double pixelResolution = mapGroup["PixelResolution"];
1140
1141 // Get the upper left corner
1142 double upperLeftX = mapGroup["UpperLeftCornerX"];
1143 double upperLeftY = mapGroup["UpperLeftCornerY"];
1144
1145 // Initialize the rest of the projection
1146 proj = (Isis::RingPlaneProjection *) RingsCreate(label, true);
1147
1148 // Create a mapper to transform pixels into projection x/y and vice versa
1149 PFPixelMapper *pixelMapper = new PFPixelMapper(pixelResolution, upperLeftX, upperLeftY);
1150 proj->SetWorldMapper(pixelMapper);
1151
1153 Displacement(upperLeftY, Displacement::Meters));
1154 }
1155 catch (IException &e) {
1156 QString msg = "Unable to initialize cube projection";
1157 if (label.fileName() != "") msg += " from file [" + label.fileName() + "]";
1158 throw IException(e, IException::Unknown, msg, _FILEINFO_);
1159 }
1160 return (Isis::Projection *) proj;
1161 }
1162} //end namespace isis
virtual double Line() const
Returns the current line number.
Definition Camera.cpp:2740
int Lines() const
Returns the number of lines in the image.
Definition Camera.cpp:2816
virtual double Sample() const
Returns the current sample number.
Definition Camera.cpp:2720
virtual bool SetImage(const double sample, const double line)
Sets the sample/line values of the image to get the lat/lon values.
Definition Camera.cpp:156
virtual bool SetUniversalGround(const double latitude, const double longitude)
Sets the lat/lon values to get the sample/line values.
Definition Camera.cpp:382
virtual bool IsBandIndependent()
Virtual method that checks if the band is independent.
Definition Camera.cpp:2679
int Samples() const
Returns the number of samples in the image.
Definition Camera.cpp:2806
virtual void SetBand(const int band)
Virtual method that sets the band number.
Definition Camera.cpp:2710
int Bands() const
Returns the number of bands in the image.
Definition Camera.cpp:2826
IO Handler for Isis Cubes.
Definition Cube.h:168
Pvl * label() const
Returns a pointer to the IsisLabel object associated with the cube.
Definition Cube.cpp:1734
Displacement is a signed length, usually in meters.
@ Meters
The distance is being specified in meters.
double meters() const
Get the distance in meters.
Definition Distance.cpp:85
File name manipulation and expansion.
Definition FileName.h:100
bool fileExists() const
Returns true if the file exists; false otherwise.
Definition FileName.cpp:449
QString expanded() const
Returns a QString of the full file name including the file path, excluding the attributes.
Definition FileName.cpp:196
Isis exception class.
Definition IException.h:91
@ Unknown
A type of error that cannot be classified as any of the other error types.
Definition IException.h:118
@ Io
A type of error that occurred when performing an actual I/O operation.
Definition IException.h:155
void append(const IException &exceptionSource)
Appends the given exception (and its list of previous exceptions) to this exception's causational exc...
Loads plugins from a shared library.
Definition Plugin.h:55
static Isis::Projection * CreateFromCube(Isis::Cube &cube)
This method is a helper method.
static Isis::Projection * RingsCreate(Isis::Pvl &label, bool allowDefaults=false)
This method returns a pointer to a RingPlaneProjection object.
static Isis::Projection * Create(Isis::Pvl &label, bool allowDefaults=false)
This method returns a pointer to a Projection object.
static Isis::Projection * RingsCreateForCube(Isis::Pvl &label, int &samples, int &lines, bool sizeMatch)
This method creates a projection for a cube to a ring plane given a label.
static Isis::Projection * RingsCreateFromCube(Isis::Cube &cube)
This method is a helper method.
static Isis::Projection * CreateForCube(Isis::Pvl &label, int &ns, int &nl, bool sizeMatch=true)
This method creates a map projection for a cube given a label.
Base class for Map Projections.
Definition Projection.h:155
void SetWorldMapper(WorldMapper *mapper)
If desired the programmer can use this method to set a world mapper to be used in the SetWorld,...
double XCoord() const
This returns the projection X provided SetGround, SetCoordinate, SetUniversalGround,...
virtual bool HasGroundRange() const
This indicates if the longitude direction type is positive west (as opposed to postive east).
void SetUpperLeftCorner(const Displacement &x, const Displacement &y)
This method is used to find the XY range for oblique aspect projections (non-polar projections) by "w...
double YCoord() const
This returns the projection Y provided SetGround, SetCoordinate, SetUniversalGround,...
bool IsGood() const
This indicates if the last invocation of SetGround, SetCoordinate, SetUniversalGround,...
QString fileName() const
Returns the filename used to initialise the Pvl object.
void addKeyword(const PvlKeyword &keyword, const InsertMode mode=Append)
Add a keyword to the container.
Contains multiple PvlContainers.
Definition PvlGroup.h:41
Container for cube-like labels.
Definition Pvl.h:119
A single keyword-value pair.
Definition PvlKeyword.h:87
bool hasObject(const QString &name) const
Returns a boolean value based on whether the object exists in the current PvlObject or not.
Definition PvlObject.h:323
@ Traverse
Search child objects.
Definition PvlObject.h:158
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition PvlObject.h:129
Base class for Map Projections of plane shapes.
virtual double TrueScaleRingRadius() const
This method returns the radius of true scale.
virtual bool XYRange(double &minX, double &maxX, double &minY, double &maxY)
This method is used to determine the x/y range which completely covers the area of interest specified...
virtual bool SetGround(const double ringRadius, const double ringLongitude)
This method is used to set the ring radius/longitude (assumed to be of the correct LatitudeType,...
virtual double UniversalLatitude() const
Returns the planetocentric latitude, in degrees, at the surface intersection point in the body fixed ...
Definition Sensor.cpp:212
virtual double UniversalLongitude() const
Returns the positive east, 0-360 domain longitude, in degrees, at the surface intersection point in t...
Definition Sensor.cpp:235
Distance LocalRadius() const
Returns the local radius at the intersection point.
Definition Sensor.cpp:269
Base class for Map TProjections.
virtual bool XYRange(double &minX, double &maxX, double &minY, double &maxY)
This method is used to determine the x/y range which completely covers the area of interest specified...
double LocalRadius(double lat) const
This method returns the local radius in meters at the specified latitude position.
virtual bool SetUniversalGround(const double lat, const double lon)
This method is used to set the latitude/longitude which must be Planetocentric (latitude) and Positiv...
virtual double TrueScaleLatitude() const
This method returns the latitude of true scale.
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition IString.cpp:211
const double Null
Value for an Isis Null pixel.
const double PI
The mathematical constant PI.
Definition Constants.h:40
Namespace for the standard library.