Official websites use .gov
A .gov website belongs to an official government organization in the United States.

Secure .gov websites use HTTPS
A lock ( ) or https:// means you’ve safely connected to the .gov website. Share sensitive information only on official, secure websites.

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
51 Isis::Projection *ProjectionFactory::Create(Isis::Pvl &label,
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
119 Isis::Projection *ProjectionFactory::RingsCreate(Isis::Pvl &label,
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
188 Isis::Projection *ProjectionFactory::CreateForCube(Isis::Pvl &label,
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
393 proj->SetUpperLeftCorner(Displacement(upperLeftX, Displacement::Meters),
394 Displacement(upperLeftY, Displacement::Meters));
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
432 Isis::Projection *ProjectionFactory::RingsCreateForCube(Isis::Pvl &label,
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
438 Isis::RingPlaneProjection *proj = (Isis::RingPlaneProjection *) RingsCreate(label, true);
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
582 proj->SetUpperLeftCorner(Displacement(upperLeftX, Displacement::Meters),
583 Displacement(upperLeftY, Displacement::Meters));
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
621 Isis::Projection *ProjectionFactory::CreateForCube(Isis::Pvl &label,
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
836 proj->SetUpperLeftCorner(Displacement(upperLeftX, Displacement::Meters),
837 Displacement(upperLeftY, Displacement::Meters));
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
872 Isis::Projection *ProjectionFactory::RingsCreateForCube(Isis::Pvl &label,
873 int &samples, int &lines, Camera &cam) {
874
875 // Create a temporary projection
876 Isis::RingPlaneProjection *proj = (Isis::RingPlaneProjection *) RingsCreate(label, true);
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
1050 proj->SetUpperLeftCorner(Displacement(upperLeftX, Displacement::Meters),
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
1095 Isis::Projection *ProjectionFactory::CreateFromCube(Isis::Pvl &label) {
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
1113 proj->SetUpperLeftCorner(Displacement(upperLeftX, Displacement::Meters),
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
1134 Isis::Projection *ProjectionFactory::RingsCreateFromCube(Isis::Pvl &label) {
1135 Isis::RingPlaneProjection *proj;
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
1152 proj->SetUpperLeftCorner(Displacement(upperLeftX, Displacement::Meters),
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
Isis::Camera::Line
virtual double Line() const
Returns the current line number.
Definition Camera.cpp:2740
Isis::Camera::Lines
int Lines() const
Returns the number of lines in the image.
Definition Camera.cpp:2816
Isis::Camera::Sample
virtual double Sample() const
Returns the current sample number.
Definition Camera.cpp:2720
Isis::Camera::SetImage
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
Isis::Camera::SetUniversalGround
virtual bool SetUniversalGround(const double latitude, const double longitude)
Sets the lat/lon values to get the sample/line values.
Definition Camera.cpp:382
Isis::Camera::IsBandIndependent
virtual bool IsBandIndependent()
Virtual method that checks if the band is independent.
Definition Camera.cpp:2679
Isis::Camera::Samples
int Samples() const
Returns the number of samples in the image.
Definition Camera.cpp:2806
Isis::Camera::SetBand
virtual void SetBand(const int band)
Virtual method that sets the band number.
Definition Camera.cpp:2710
Isis::Camera::Bands
int Bands() const
Returns the number of bands in the image.
Definition Camera.cpp:2826
Isis::Cube
IO Handler for Isis Cubes.
Definition Cube.h:168
Isis::Cube::label
Pvl * label() const
Returns a pointer to the IsisLabel object associated with the cube.
Definition Cube.cpp:1975
Isis::Displacement
Displacement is a signed length, usually in meters.
Definition Displacement.h:31
Isis::Displacement::Meters
@ Meters
The distance is being specified in meters.
Definition Displacement.h:40
Isis::Distance::meters
double meters() const
Get the distance in meters.
Definition Distance.cpp:85
Isis::Plugin
Loads plugins from a shared library.
Definition Plugin.h:55
Isis::ProjectionFactory::CreateFromCube
static Isis::Projection * CreateFromCube(Isis::Cube &cube)
This method is a helper method.
Definition ProjectionFactory.cpp:1069
Isis::ProjectionFactory::RingsCreate
static Isis::Projection * RingsCreate(Isis::Pvl &label, bool allowDefaults=false)
This method returns a pointer to a RingPlaneProjection object.
Definition ProjectionFactory.cpp:119
Isis::ProjectionFactory::Create
static Isis::Projection * Create(Isis::Pvl &label, bool allowDefaults=false)
This method returns a pointer to a Projection object.
Definition ProjectionFactory.cpp:51
Isis::ProjectionFactory::RingsCreateForCube
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.
Definition ProjectionFactory.cpp:432
Isis::ProjectionFactory::RingsCreateFromCube
static Isis::Projection * RingsCreateFromCube(Isis::Cube &cube)
This method is a helper method.
Definition ProjectionFactory.cpp:1081
Isis::ProjectionFactory::CreateForCube
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.
Definition ProjectionFactory.cpp:188
Isis::Projection
Base class for Map Projections.
Definition Projection.h:154
Isis::Projection::SetWorldMapper
void SetWorldMapper(WorldMapper *mapper)
If desired the programmer can use this method to set a world mapper to be used in the SetWorld,...
Definition Projection.cpp:468
Isis::Projection::XCoord
double XCoord() const
This returns the projection X provided SetGround, SetCoordinate, SetUniversalGround,...
Definition Projection.cpp:381
Isis::Projection::HasGroundRange
virtual bool HasGroundRange() const
This indicates if the longitude direction type is positive west (as opposed to postive east).
Definition Projection.cpp:343
Isis::Projection::SetUpperLeftCorner
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...
Definition Projection.cpp:1367
Isis::Projection::YCoord
double YCoord() const
This returns the projection Y provided SetGround, SetCoordinate, SetUniversalGround,...
Definition Projection.cpp:394
Isis::Projection::IsGood
bool IsGood() const
This indicates if the last invocation of SetGround, SetCoordinate, SetUniversalGround,...
Definition Projection.cpp:368
Isis::RingPlaneProjection
Base class for Map Projections of plane shapes.
Definition RingPlaneProjection.h:147
Isis::RingPlaneProjection::TrueScaleRingRadius
virtual double TrueScaleRingRadius() const
This method returns the radius of true scale.
Definition RingPlaneProjection.cpp:199
Isis::RingPlaneProjection::XYRange
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...
Definition RingPlaneProjection.cpp:656
Isis::RingPlaneProjection::SetGround
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,...
Definition RingPlaneProjection.cpp:456
Isis::Sensor::UniversalLatitude
virtual double UniversalLatitude() const
Returns the planetocentric latitude, in degrees, at the surface intersection point in the body fixed ...
Definition Sensor.cpp:210
Isis::Sensor::UniversalLongitude
virtual double UniversalLongitude() const
Returns the positive east, 0-360 domain longitude, in degrees, at the surface intersection point in t...
Definition Sensor.cpp:233
Isis::Sensor::LocalRadius
Distance LocalRadius() const
Returns the local radius at the intersection point.
Definition Sensor.cpp:267
Isis::TProjection
Base class for Map TProjections.
Definition TProjection.h:165
Isis::TProjection::XYRange
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...
Definition TProjection.cpp:985
Isis::TProjection::LocalRadius
double LocalRadius(double lat) const
This method returns the local radius in meters at the specified latitude position.
Definition TProjection.cpp:318
Isis::TProjection::SetUniversalGround
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...
Definition TProjection.cpp:831
Isis::TProjection::TrueScaleLatitude
virtual double TrueScaleLatitude() const
This method returns the latitude of true scale.
Definition TProjection.cpp:360
Isis
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
Isis::toString
QString toString(const LinearAlgebra::Vector &vector, int precision)
A global function to format LinearAlgebra::Vector as a QString with the given precision.
Definition LinearAlgebra.cpp:1852
std
Namespace for the standard library.