Select Git revision
-
acpaquette authoredacpaquette authored
cam2map.cpp 21.73 KiB
#include "cam2map.h"
#include "Camera.h"
#include "CubeAttribute.h"
#include "IException.h"
#include "IString.h"
#include "ProjectionFactory.h"
#include "PushFrameCameraDetectorMap.h"
#include "Pvl.h"
#include "Target.h"
#include "TProjection.h"
using namespace std;
namespace Isis {
// Global variables
void bandChange(const int band);
Cube *icube;
Camera *incam;
void cam2map(UserInterface &ui, Pvl *log) {
// Open the input cube
Cube icube;
CubeAttributeInput inAtt = ui.GetInputAttribute("FROM");
if (inAtt.bands().size() != 0) {
icube.setVirtualBands(inAtt.bands());
}
icube.open(ui.GetFileName("FROM"));
// Get the map projection file provided by the user
Pvl userMap;
userMap.read(ui.GetFileName("MAP"));
PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse);
cam2map(&icube, userMap, userGrp, ui, log);
}
void cam2map(Cube *icube, Pvl &userMap, PvlGroup &userGrp, UserInterface &ui, Pvl *log) {
ProcessRubberSheet p;
cam2map(icube, userMap, userGrp, p, ui, log);
}
void cam2map(Cube *icube, Pvl &userMap, PvlGroup &userGrp, ProcessRubberSheet &p,
UserInterface &ui, Pvl *log){
// Get the camera from the input cube
p.SetInputCube(icube);
incam = icube->camera();
// Make sure it is not the sky
if (incam->target()->isSky()) {
QString msg = "The image [" + ui.GetFileName("FROM") +
"] is targeting the sky, use skymap instead.";
throw IException(IException::User, msg, _FILEINFO_);
}
// Get the mapping grop
Pvl camMap;
incam->BasicMapping(camMap);
PvlGroup &camGrp = camMap.findGroup("Mapping");
// Make the target info match the user mapfile
double minlat, maxlat, minlon, maxlon;
incam->GroundRange(minlat, maxlat, minlon, maxlon, userMap);
camGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace);
camGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace);
camGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace); camGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace);
// We want to delete the keywords we just added if the user wants the range
// out of the mapfile, otherwise they will replace any keywords not in the
// mapfile
// Use the updated label to create the output projection
int samples, lines;
TProjection *outmap = NULL;
bool trim = ui.GetBoolean("TRIM");
bool occlusion = ui.GetBoolean("OCCLUSION");
// Make sure the target name of the input cube and map file match.
if (userGrp.hasKeyword("TargetName") && !icube->group("Instrument").findKeyword("TargetName").isNull()) {
if (!PvlKeyword::stringEqual(incam->target()->name(), userGrp.findKeyword("TargetName")[0])) {
QString msg = "The TargetName: [" + incam->target()->name() + "] of the input cube: [" + icube->fileName() +
"] does not match the TargetName: [" + userGrp.findKeyword("TargetName")[0] + "] of the map file: [" +
ui.GetFileName("MAP") + "].";
throw IException(IException::User, msg, _FILEINFO_);
}
}
if ( !ui.GetBoolean("MATCHMAP") ) {
if (ui.GetString("DEFAULTRANGE") == "MAP") {
camGrp.deleteKeyword("MinimumLatitude");
camGrp.deleteKeyword("MaximumLatitude");
camGrp.deleteKeyword("MinimumLongitude");
camGrp.deleteKeyword("MaximumLongitude");
}
// Otherwise, remove the keywords from the map file so the camera keywords
// will be propogated correctly
else {
while(userGrp.hasKeyword("MinimumLatitude")) {
userGrp.deleteKeyword("MinimumLatitude");
}
while(userGrp.hasKeyword("MinimumLongitude")) {
userGrp.deleteKeyword("MinimumLongitude");
}
while(userGrp.hasKeyword("MaximumLatitude")) {
userGrp.deleteKeyword("MaximumLatitude");
}
while(userGrp.hasKeyword("MaximumLongitude")) {
userGrp.deleteKeyword("MaximumLongitude");
}
}
// If the user decided to enter a ground range then override
if ( ui.WasEntered("MINLON") ) {
userGrp.addKeyword(PvlKeyword("MinimumLongitude",
toString(ui.GetDouble("MINLON"))), Pvl::Replace);
}
if ( ui.WasEntered("MAXLON") ) {
userGrp.addKeyword(PvlKeyword("MaximumLongitude",
toString(ui.GetDouble("MAXLON"))), Pvl::Replace);
}
if ( ui.WasEntered("MINLAT") ) {
userGrp.addKeyword(PvlKeyword("MinimumLatitude",
toString(ui.GetDouble("MINLAT"))), Pvl::Replace);
}
if ( ui.WasEntered("MAXLAT") ) {
userGrp.addKeyword(PvlKeyword("MaximumLatitude",
toString(ui.GetDouble("MAXLAT"))), Pvl::Replace);
}
// If they want the res. from the mapfile, delete it from the camera so
// nothing gets overriden if (ui.GetString("PIXRES") == "MAP") {
camGrp.deleteKeyword("PixelResolution");
}
// Otherwise, delete any resolution keywords from the mapfile so the camera
// info is propogated over
else if (ui.GetString("PIXRES") == "CAMERA") {
if (userGrp.hasKeyword("Scale")) {
userGrp.deleteKeyword("Scale");
}
if (userGrp.hasKeyword("PixelResolution")) {
userGrp.deleteKeyword("PixelResolution");
}
}
// Copy any defaults that are not in the user map from the camera map file
for (int k = 0; k < camGrp.keywords(); k++) {
if (!userGrp.hasKeyword(camGrp[k].name())) {
userGrp += camGrp[k];
}
}
// If the user decided to enter a resolution then override
if (ui.WasEntered("PIXRES")) {
if (ui.GetString("PIXRES") == "MPP") {
userGrp.addKeyword(PvlKeyword("PixelResolution",
toString(ui.GetDouble("RESOLUTION"))),
Pvl::Replace);
if (userGrp.hasKeyword("Scale")) {
userGrp.deleteKeyword("Scale");
}
}
else if (ui.GetString("PIXRES") == "PPD") {
userGrp.addKeyword(PvlKeyword("Scale",
toString(ui.GetDouble("RESOLUTION"))),
Pvl::Replace);
if (userGrp.hasKeyword("PixelResolution")) {
userGrp.deleteKeyword("PixelResolution");
}
}
}
// See if the user want us to handle the longitude seam
if ( (ui.GetString("DEFAULTRANGE") == "CAMERA" || ui.GetString("DEFAULTRANGE") == "MINIMIZE") ) {
if (incam->IntersectsLongitudeDomain(userMap)) {
if (ui.GetString("LONSEAM") == "AUTO") {
if ((int) userGrp["LongitudeDomain"] == 360) {
userGrp.addKeyword(PvlKeyword("LongitudeDomain", "180"), Pvl::Replace);
if (incam->IntersectsLongitudeDomain(userMap)) {
// Its looks like a global image so switch back to the users preference
userGrp.addKeyword(PvlKeyword("LongitudeDomain", "360"), Pvl::Replace);
}
}
else {
userGrp.addKeyword(PvlKeyword("LongitudeDomain", "360"), Pvl::Replace);
if (incam->IntersectsLongitudeDomain(userMap)) {
// Its looks like a global image so switch back to the
// users preference
userGrp.addKeyword(PvlKeyword("LongitudeDomain", "180"), Pvl::Replace);
}
}
// Make the target info match the new longitude domain
double minlat, maxlat, minlon, maxlon;
incam->GroundRange(minlat, maxlat, minlon, maxlon, userMap);
if (!ui.WasEntered("MINLAT")) {
userGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace);
}
if (!ui.WasEntered("MAXLAT")) {
userGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace); }
if (!ui.WasEntered("MINLON")) {
userGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace);
}
if (!ui.WasEntered("MAXLON")) {
userGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace);
}
}
else if (ui.GetString("LONSEAM") == "ERROR") {
QString msg = "The image [" + ui.GetFileName("FROM") + "] crosses the " +
"longitude seam";
throw IException(IException::User, msg, _FILEINFO_);
}
}
}
// Determine the image size
if (ui.GetString("DEFAULTRANGE") == "MINIMIZE") {
outmap = (TProjection *) ProjectionFactory::CreateForCube(userMap, samples, lines, *incam);
trim = false;
}
else {//if (ui.GetString("DEFAULTRANGE") == "CAMERA" || DEFAULTRANGE = MAP) {
outmap = (TProjection *) ProjectionFactory::CreateForCube(userMap, samples, lines, false);
}
// else {
// outmap = (TProjection *) ProjectionFactory::CreateForCube(userMap, samples, lines, false);
// }
}
else { // MATCHMAP=TRUE
//does this have any affect on anything?
camGrp.deleteKeyword("MinimumLatitude");
camGrp.deleteKeyword("MaximumLatitude");
camGrp.deleteKeyword("MinimumLongitude");
camGrp.deleteKeyword("MaximumLongitude");
camGrp.deleteKeyword("PixelResolution");
// Copy any defaults that are not in the user map from the camera map file
outmap = (TProjection *) ProjectionFactory::CreateForCube(userMap,
samples,
lines,
true);//change usrmap to camgrp?
}
// Output the mapping group used to the Gui session log
PvlGroup cleanMapping = outmap->Mapping();
// Allocate the output cube and add the mapping labels
QString fname = ui.GetFileName("TO");
Isis::CubeAttributeOutput &atts = ui.GetOutputAttribute("TO");
Cube *ocube = p.SetOutputCube(fname, atts, samples, lines, icube->bandCount());
ocube->putGroup(cleanMapping);
// Set up the interpolator
Interpolator *interp = NULL;
if (ui.GetString("INTERP") == "NEARESTNEIGHBOR") {
interp = new Interpolator(Interpolator::NearestNeighborType);
}
else if (ui.GetString("INTERP") == "BILINEAR") {
interp = new Interpolator(Interpolator::BiLinearType);
}
else if (ui.GetString("INTERP") == "CUBICCONVOLUTION") {
interp = new Interpolator(Interpolator::CubicConvolutionType);
}
// See if we need to deal with band dependent camera models
if (!incam->IsBandIndependent()) {
p.BandChange(bandChange);
}
// See if center of input image projects. If it does, force tile
// containing this center to be processed in ProcessRubberSheet.
// TODO: WEIRD ... why is this needed ... Talk to Tracie ... JAA??
double centerSamp = icube->sampleCount() / 2.;
double centerLine = icube->lineCount() / 2.;
if (incam->SetImage(centerSamp, centerLine)) {
if (outmap->SetUniversalGround(incam->UniversalLatitude(),
incam->UniversalLongitude())) {
p.ForceTile(outmap->WorldX(), outmap->WorldY());
}
}
// Create an alpha cube group for the output cube
if (!ocube->hasGroup("AlphaCube")) {
PvlGroup alpha("AlphaCube");
alpha += PvlKeyword("AlphaSamples", toString(icube->sampleCount()));
alpha += PvlKeyword("AlphaLines", toString(icube->lineCount()));
alpha += PvlKeyword("AlphaStartingSample", toString(0.5));
alpha += PvlKeyword("AlphaStartingLine", toString(0.5));
alpha += PvlKeyword("AlphaEndingSample", toString(icube->sampleCount() + 0.5));
alpha += PvlKeyword("AlphaEndingLine", toString(icube->lineCount() + 0.5));
alpha += PvlKeyword("BetaSamples", toString(icube->sampleCount()));
alpha += PvlKeyword("BetaLines", toString(icube->lineCount()));
ocube->putGroup(alpha);
}
// We will need a transform class
Transform *transform = 0;
// Okay we need to decide how to apply the rubbersheeting for the transform
// Does the user want to define how it is done?
if (ui.GetString("WARPALGORITHM") == "FORWARDPATCH") {
transform = new cam2mapForward(icube->sampleCount(),
icube->lineCount(),
incam,
samples,
lines,
outmap,
trim);
int patchSize = ui.GetInteger("PATCHSIZE");
if (patchSize <= 1) {
patchSize = 3; // Make the patchsize reasonable
}
p.setPatchParameters(1, 1, patchSize, patchSize, patchSize-1, patchSize-1);
p.processPatchTransform(*transform, *interp);
}
else if (ui.GetString("WARPALGORITHM") == "REVERSEPATCH") {
transform = new cam2mapReverse(icube->sampleCount(),
icube->lineCount(), incam, samples,lines,
outmap, trim, occlusion);
int patchSize = ui.GetInteger("PATCHSIZE");
int minPatchSize = 4;
if (patchSize < minPatchSize) {
patchSize = minPatchSize;
}
p.SetTiling(patchSize, patchSize);
p.StartProcess(*transform, *interp);
}
// The user didn't want to override the program smarts.
// Handle framing cameras. Always process using the backward
// driven system (tfile).
else if (incam->GetCameraType() == Camera::Framing) {
transform = new cam2mapReverse(icube->sampleCount(),
icube->lineCount(), incam, samples,lines, outmap, trim, occlusion);
p.SetTiling(4, 4);
p.StartProcess(*transform, *interp);
}
// The user didn't want to override the program smarts.
// Handle linescan cameras. Always process using the forward
// driven patch option. Faster and we get better orthorectification
//
// TODO: For now use the default patch size. Need to modify
// to determine patch size based on 1) if the limb is in the file
// or 2) if the DTM is much coarser than the image
else if (incam->GetCameraType() == Camera::LineScan) {
transform = new cam2mapForward(icube->sampleCount(),
icube->lineCount(), incam, samples,lines,
outmap, trim);
p.processPatchTransform(*transform, *interp);
}
// The user didn't want to override the program smarts.
// Handle pushframe cameras. Always process using the forward driven patch
// option. It is much faster than the tfile method. We will need to
// determine patch sizes based on the size of the push frame.
//
// TODO: What if the user has run crop, enlarge, or shrink on the push
// frame cube. Things probably won't work unless they do it just right
// TODO: What about the THEMIS VIS Camera. Will tall narrow (128x4) patches
// work okay?
else if (incam->GetCameraType() == Camera::PushFrame) {
transform = new cam2mapForward(icube->sampleCount(),
icube->lineCount(), incam, samples,lines,
outmap, trim);
// Get the frame height
PushFrameCameraDetectorMap *dmap = (PushFrameCameraDetectorMap *) incam->DetectorMap();
int frameSize = dmap->frameletHeight() / dmap->LineScaleFactor();
// Check for even/odd cube to determine starting line
PvlGroup &instGrp = icube->label()->findGroup("Instrument", Pvl::Traverse);
int startLine = 1;
// Get the alpha cube group in case they cropped the image
AlphaCube acube(*icube);
double betaLine = acube.AlphaLine(1.0);
if (fabs(betaLine - 1.0) > 0.0000000001) {
if (fabs(betaLine - (int) betaLine) > 0.00001) {
string msg = "Input file is a pushframe camera cropped at a ";
msg += "fractional pixel. Can not project";
throw IException(IException::User, msg, _FILEINFO_);
}
int offset = (((int) (betaLine + 0.5)) - 1) % frameSize;
startLine -= offset;
}
if (((QString)instGrp["Framelets"]).toUpper() == "EVEN") {
startLine += frameSize;
}
p.setPatchParameters(1, startLine, 5, frameSize,
4, frameSize * 2);
p.processPatchTransform(*transform, *interp);
}
// The user didn't want to override the program smarts. The other camera
// types have not be analyized. This includes Radar and Point. Continue to
// use the reverse geom option with the default tiling hints
else {
transform = new cam2mapReverse(icube->sampleCount(), icube->lineCount(), incam, samples,lines,
outmap, trim, occlusion);
int tileStart, tileEnd;
incam->GetGeometricTilingHint(tileStart, tileEnd);
p.SetTiling(tileStart, tileEnd);
p.StartProcess(*transform, *interp);
}
// Wrap up the warping process
p.EndProcess();
// add mapping to print.prt
if(log) {
log->addGroup(cleanMapping);
}
// Cleanup
delete outmap;
delete transform;
delete interp;
}
// Transform object constructor
cam2mapForward::cam2mapForward(const int inputSamples, const int inputLines,
Camera *incam, const int outputSamples,
const int outputLines, TProjection *outmap,
bool trim) {
p_inputSamples = inputSamples;
p_inputLines = inputLines;
p_incam = incam;
p_outputSamples = outputSamples;
p_outputLines = outputLines;
p_outmap = outmap;
p_trim = trim;
}
// Transform method mapping input line/samps to lat/lons to output line/samps
bool cam2mapForward::Xform(double &outSample, double &outLine,
const double inSample, const double inLine) {
// See if the input image coordinate converts to a lat/lon
if (!p_incam->SetImage(inSample,inLine)) return false;
// Does that ground coordinate work in the map projection
double lat = p_incam->UniversalLatitude();
double lon = p_incam->UniversalLongitude();
if (!p_outmap->SetUniversalGround(lat,lon)) return false;
// See if we should trim
if ((p_trim) && (p_outmap->HasGroundRange())) {
if (p_outmap->Latitude() < p_outmap->MinimumLatitude()) return false;
if (p_outmap->Latitude() > p_outmap->MaximumLatitude()) return false;
if (p_outmap->Longitude() < p_outmap->MinimumLongitude()) return false;
if (p_outmap->Longitude() > p_outmap->MaximumLongitude()) return false;
}
// Get the output sample/line coordinate
outSample = p_outmap->WorldX();
outLine = p_outmap->WorldY();
// Make sure the point is inside the output image
if (outSample < 0.5) return false;
if (outLine < 0.5) return false;
if (outSample > p_outputSamples + 0.5) return false;
if (outLine > p_outputLines + 0.5) return false;
// Everything is good return true;
}
int cam2mapForward::OutputSamples() const {
return p_outputSamples;
}
int cam2mapForward::OutputLines() const {
return p_outputLines;
}
// Transform object constructor
cam2mapReverse::cam2mapReverse(const int inputSamples, const int inputLines,
Camera *incam, const int outputSamples,
const int outputLines, TProjection *outmap,
bool trim, bool occlusion) {
p_inputSamples = inputSamples;
p_inputLines = inputLines;
p_incam = incam;
p_outputSamples = outputSamples;
p_outputLines = outputLines;
p_outmap = outmap;
p_trim = trim;
p_occlusion = occlusion;
}
// Transform method mapping output line/samps to lat/lons to input line/samps
bool cam2mapReverse::Xform(double &inSample, double &inLine,
const double outSample, const double outLine) {
// See if the output image coordinate converts to lat/lon
if (!p_outmap->SetWorld(outSample, outLine)) return false;
// See if we should trim
if ((p_trim) && (p_outmap->HasGroundRange())) {
if (p_outmap->Latitude() < p_outmap->MinimumLatitude()) return false;
if (p_outmap->Latitude() > p_outmap->MaximumLatitude()) return false;
if (p_outmap->Longitude() < p_outmap->MinimumLongitude()) return false;
if (p_outmap->Longitude() > p_outmap->MaximumLongitude()) return false;
}
// Get the universal lat/lon and see if it can be converted to input line/samp
double lat = p_outmap->UniversalLatitude();
double lon = p_outmap->UniversalLongitude();
if (!p_incam->SetUniversalGround(lat, lon)) return false;
// Make sure the point is inside the input image
if (p_incam->Sample() < 0.5) return false;
if (p_incam->Line() < 0.5) return false;
if (p_incam->Sample() > p_inputSamples + 0.5) return false;
if (p_incam->Line() > p_inputLines + 0.5) return false;
// Everything is good
inSample = p_incam->Sample();
inLine = p_incam->Line();
// Good to ground one last time to check for occlusion
p_incam->SetImage(inSample, inLine);
if (p_occlusion){
if (abs(lat - p_incam->UniversalLatitude()) > 0.00001 || abs(lon - p_incam->UniversalLongitude()) > 0.00001) {
return false;
}
}
return true;
}
int cam2mapReverse::OutputSamples() const {
return p_outputSamples;
}
int cam2mapReverse::OutputLines() const {
return p_outputLines;
}
void bandChange(const int band) {
incam->SetBand(band);
}
}