diff --git a/isis/src/apollo/apps/apollopanstitcher/apollopanstitcher.cpp b/isis/src/apollo/apps/apollopanstitcher/apollopanstitcher.cpp new file mode 100644 index 0000000000..97b5809a13 --- /dev/null +++ b/isis/src/apollo/apps/apollopanstitcher/apollopanstitcher.cpp @@ -0,0 +1,787 @@ +/* +Program uses the fiducial marks along the edges of Apollo Panoramic Images to stich eight sub-scans into a single conitinuous image + +*/ +#include "apollopanstitcher.h" + +//C++ standard libraries if needed +#include +#include + +//QT libraries if needed if needed +#include + +//third party libraries if needed +#include + +//Isis Headers if needed +#include "AutoReg.h" +#include "AutoRegFactory.h" +#include "Centroid.h" +#include "CentroidApolloPan.h" +#include "Chip.h" +#include "Cube.h" +#include "CubeAttribute.h" +#include "Interpolator.h" +#include "ProcessMosaic.h" +#include "ProcessRubberSheet.h" +#include "Pvl.h" +#include "PvlGroup.h" +#include "PvlKeyword.h" +#include "PvlObject.h" +#include "Table.h" +#include "Trans2d3p.h" + +/** + * + * @internal + * @History 2011-11-02 Orrin Thomas - Original Version. + * + * @History 2017-06-29 Christopher Combs - Moved initialization of trans[7] to before + * the loop that uses it. Fixes #4948. + * + */ + +using namespace std; +using namespace Isis; + +//constants for 5 micro resolution scans +#define SCALE 10.0 //scale used for down sizingS patternS an search chips +#define SEARCHh 1400.0 //number of lines (in 5-micron-pixels) in search space for the first + // fiducial +#define SEARCHc 350.0 //number of samples per edge(in 5-micron-pixels) in each sub-search area +#define AVERs 5286.0 //average smaples (in 5-micron-pixels) between fiducials +#define AVERl 23459.0 //average diference (in 5-micron-pixels) between the top and bottom + // fiducials +#define TRANS_N 28520.0 //nomimal dx between scan lines, scani x + 28532 = + // (approx) scani+1 x + // also the size of the search area for the first fiducial + +double R_MOON[3]; + +typedef struct TRANS2d3p { //2d 3 parameter transformation + double theta; + double dx; + double dy; + double averR; //average residual + double maxR; //max resiudal + double limit[2]; //where to stop using this transformation, will eventually be a point in the + //coordinate system of the new cube, a vertical line through this point will + //be the dividing line +}TRANS2d3p; + +namespace Isis { + void apolloPanStitcher(UserInterface &ui) { + Chip patternS,searchS; //scaled pattern and search chips + Cube *panC[8], //panoramic image + outputC, //outputC + fidC; //Fiducial image + int regStatus, + scanS, //samples in a sub scan + scanL, //lines in a sub scan + refL, //number of lines in the patternS + refS; //number of samples in the patternS + + Pvl pvl; + + bool foundFirst; + + int i,j,k; + + double sampleFrom, sampleTo; //limits for rubber sheeting and mosaicing the sub-scans + double play; //how close a pixel must be to be considered on the ellipse + + //vector of possible solutions to a transformation-all transformations resulting from any + // possible combination of minimal sets of data are calculated to find the best solution + // in an outlier resistant way + std::vector solV; + TRANS2d3p sol; //individual solution to be add to the above vector + + TRANS2d3p trans[8]; //final solutions for the transformations-one for each sub-scan + + QString fileName,tempIString; + QString fileBaseName, tempString; + + double l = 1, s = 1, sample, line, temp,//line and sample coordinates for looping through the panC + ct, //"cosine of theta" used to prevent evaluted the trig function thousands of times + st; //"sine of theta" used to prevent evaluted the trig function thousands of times + + std::vector > scanFid; //fiducials of Scan i + std::vector > conFid; //fiducials of Scan i + std::vector fid(2); //individual fiducial measurement + std::vector pairFid(4); //pair of conjugate fiducial coordinates + + //read the image resolutions and scale the constants acordingly + double resolution = ui.GetDouble("MICRONS"), //pixel size in microns + scale = SCALE *5.0/resolution, //reduction scale for quicker autoregistrations + searchHeight = SEARCHh*5.0/resolution, //number of lines (in 5-micron-pixels) in + // search space for the first fiducial + searchCellSize = SEARCHc*5.0/resolution, //height/width of search chips block + averageSamples = AVERs *5.0/resolution, //scaled smaples between fiducials + nominalTrans = TRANS_N*5.0/resolution, //nominal sample direction overlap between + // successive sub scans + averageLines = AVERl *5.0/resolution; //scaled average diference between the top + // and bottom fiducials + + if (15.0 / resolution < 1.5) { + play = 1.5; + } + else { + play = 15.0 / resolution; + } + + //copy the patternS chip (the entire ApolloPanFiducialMark.cub) + FileName fiducialFileName("$apollo15/calibration/ApolloPanFiducialMark.cub"); + fidC.open(fiducialFileName.expanded(),"r"); + if (!fidC.isOpen()) { + string msg = "Unable to open the fiducial patternS cube: ApolloPanFiducialMark.cub\n"; + throw IException(IException::User, msg, _FILEINFO_); + } + refL = fidC.lineCount(); + refS = fidC.sampleCount(); + //scaled pattern chip for fast matching + patternS.SetSize((int)floor((refS-2)/SCALE), (int)floor((refL-2)/SCALE)); + patternS.TackCube((refS-1)/2, (refL-1)/2); + patternS.Load(fidC, 0, SCALE); + + //make sure the eight "FROM" cubes exist, and can can be opened for reading + fileBaseName = ui.GetString("FILE_BASE"); + for (i = 1; i <= 8; i++) { + panC[i-1] = new Cube; + fileName = fileBaseName + "-000" + toString(i) + ".cub"; + panC[i-1]->open(fileName, "r"); + if (!panC[i-1]->isOpen()) { + string msg = "Unable to open input cube: " + IString(fileName) + "\n"; + throw IException(IException::User, msg, _FILEINFO_); + } + } + + //parameters for maximum correlation autoregestration + // see: ///usgs/pkgs/isis3nightly2011-09-21/isis/doc/documents/patternSMatch/patternSMatch.html#DistanceTolerance + FileName fiducialPvl("$ISISROOT/appdata/templates/apollo/PanFiducialFinder.def"); + pvl.read(fiducialPvl.expanded()); //read in the autoreg parameters + AutoReg *arS = AutoRegFactory::Create(pvl); + + *arS->PatternChip() = patternS; //patternS chip is constant + + //set up a centroid measurer + //Centroid centroid; + CentroidApolloPan centroid(resolution); + Chip inputChip,selectionChip; + if (panC[0]->pixelType() == 1) { //UnsignedByte + centroid.setDNRange(12, 1e99); //8 bit bright target + } + else { + centroid.setDNRange(3500, 1e99); //16 bit bright target + } + inputChip.SetSize(int(ceil(200 * 5.0 / resolution)), int(ceil(200 * 5.0 / resolution))); + + //find conjugate finducials and calcualte transformations + for (i = 0; i < 7; i++) { + //for each scan segment 0-6 + //1. find the probable conjugate fiducials (those with the smallest sample coordinates) + //2. use the nominal transformation to find their matches in the next scan + //3. robustly calculate the transformation (that is in a manner resistent to erroneous + // measurements) + + //Step 1: find the probable conjugate fiducials (those with the smallest sample coordinates) + scanS = panC[i+1]->sampleCount(); + + scanFid.clear(); + + //search sizes are constanst + searchS.SetSize(int(searchCellSize/scale),int(searchCellSize/scale)); + //now start searching along a horizontal line for the first fiducial mark + foundFirst = false; + for (l = searchCellSize / 2; + l < searchHeight + searchCellSize / 2.0 && !foundFirst; + l = l + (searchCellSize - 125 * 5.0 / resolution) ) { + for (s = searchCellSize / 2; + s < (nominalTrans + searchCellSize) / 2.0 && !foundFirst; + s = s + (searchCellSize - 125 * 5.0 / resolution) ) { + searchS.TackCube(s, l); + searchS.Load(*panC[i], 0, scale); + *arS->SearchChip() = searchS; + regStatus = arS->Register(); + if (regStatus == AutoReg::SuccessPixel) { + inputChip.TackCube(arS->CubeSample(), arS->CubeLine()); + inputChip.Load(*panC[i], 0, 1); + inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); + //continuous dynamic range selection + centroid.selectAdaptive(&inputChip, &selectionChip); + //selectionChip.Write("firstselectionTemp.cub"); + //elliptical trimming/smoothing + if (centroid.elipticalReduction(&selectionChip, 95, play, 2000)) { + //center of mass to reduce selection to a single measure + centroid.centerOfMass(&selectionChip, &sample, &line); + inputChip.SetChipPosition(sample, line); + fid[0] = inputChip.CubeSample(); + fid[1] = inputChip.CubeLine(); + scanFid.push_back(fid); + foundFirst = true; //once the first fiducial is found stop + } + } + } + } + if (s >= averageLines + searchCellSize / 2.0) { + QString msg = "Unable to locate a fiducial mark in the input cube [" + panC[i]->fileName() + + "]. Check FROM and MICRONS parameters."; + throw IException(IException::Io, msg, _FILEINFO_); + return; + } + + + for (s = scanFid[0][0], l = scanFid[0][1]; s + nominalTrans < scanS; s+= averageSamples / 2.0) { + //look for the bottom fiducial + searchS.TackCube(s, l + averageLines); + searchS.Load(*panC[i], 0, scale); + *arS->SearchChip() = searchS; + regStatus = arS->Register(); + if (regStatus == AutoReg::SuccessPixel) { + inputChip.TackCube( arS->CubeSample(), arS->CubeLine()); + inputChip.Load(*panC[i], 0, 1); + inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); + } + else { + inputChip.TackCube(s, l + averageLines); + inputChip.Load(*panC[i], 0, 1); + inputChip.SetCubePosition(s, l + averageLines); + } + //continuous dynamic range selection + centroid.selectAdaptive(&inputChip, &selectionChip); + //elliptical trimming/smoothing... if this fails move on + if (centroid.elipticalReduction(&selectionChip, 95, play, 2000) != 0) { + //center of mass to reduce selection to a single measure + centroid.centerOfMass(&selectionChip,&sample,&line); + inputChip.SetChipPosition(sample,line); + fid[0] = inputChip.CubeSample(); + fid[1] = inputChip.CubeLine(); + scanFid.push_back(fid); + } + + //look for the top fiducial + //is this is the first time through the loop? + if (s == scanFid[0][0]) { + continue; //then the top fiducial was already found + } + searchS.TackCube(s, l); + searchS.Load(*panC[i], 0, scale); + *arS->SearchChip() = searchS; + regStatus = arS->Register(); + if (regStatus == AutoReg::SuccessPixel) { + inputChip.TackCube(arS->CubeSample(), arS->CubeLine()); + inputChip.Load(*panC[i], 0, 1); + inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); + } + else { + inputChip.TackCube(s, l + averageLines); + inputChip.Load(*panC[i], 0, 1); + inputChip.SetCubePosition(s, l + averageLines); + } + //continuous dynamic range selection + centroid.selectAdaptive(&inputChip, &selectionChip); + //elliptical trimming/smoothing... if this fails move on + if (centroid.elipticalReduction(&selectionChip, 95, play, 2000) != 0) { + //center of mass to reduce selection to a single measure + centroid.centerOfMass(&selectionChip, &sample, &line); + inputChip.SetChipPosition(sample, line); + //s and l are refined here to help follow any trends in the image/film + s=fid[0] = inputChip.CubeSample(); + l=fid[1] = inputChip.CubeLine(); + scanFid.push_back(fid); + } + } + + // STEP 2. use the nominal transformation to find their matches in the next scan + conFid.clear(); + for (j = scanFid.size() - 1; j >= 0; j--) { + searchS.TackCube(scanFid[j][0] + nominalTrans,scanFid[j][1]); + searchS.Load(*panC[i+1], 0, scale); + *arS->SearchChip() = searchS; + regStatus = arS->Register(); + //if autoreg is successful a smaller window can be used for the following steps + if (regStatus == AutoReg::SuccessPixel) { + inputChip.SetSize(int(200 * 5.0 / resolution), int(150 * 5.0 / resolution)); + inputChip.TackCube(arS->CubeSample(), arS->CubeLine()); + inputChip.Load(*panC[i+1], 0, 1); + inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); + } + else { //if autoreg is unsuccessful, a larger window will be used + inputChip.SetSize(int(ceil(200 * 5.0 / resolution)), int(ceil(200 * 5.0 / resolution))); + inputChip.TackCube(scanFid[j][0] + nominalTrans, scanFid[j][1]); + inputChip.Load(*panC[i+1], 0, 1); + inputChip.SetCubePosition(scanFid[j][0] + nominalTrans, scanFid[j][1]); + } + //continuous dynamic range selection + centroid.selectAdaptive(&inputChip, &selectionChip); + //elliptical trimming/smoothing... if this fails move on + if (centroid.elipticalReduction(&selectionChip, 95, play, 2000) != 0) { + //center of mass to reduce selection to a single measure + centroid.centerOfMass(&selectionChip, &sample, &line); + inputChip.SetChipPosition(sample, line); + //now that the complete conjugate pair has been located save it + pairFid[2] = inputChip.CubeSample(); + pairFid[3] = inputChip.CubeLine(); + pairFid[0] = scanFid[j][0]; + pairFid[1] = scanFid[j][1]; + conFid.push_back(pairFid); + } + } + //STEP 3. robustly calculate the transformation (that is in a manner resistent to erroneous + // measurements) + //calculate the transformation from i+1 to i + //RANSAC2 algorithim + //Find two point solutions and then attempt to generalize to more points + //this helps to identify measurement blunders by seeing if certain points are excluded + //consistently--or equivalently if certain solutions keep repeating themselves + //the limitation in this case is that there will never be more than 6 points, if the number + //of points drops to <=three blunder detection becomes impossible + + /*This is a Gauss-Helemert non-linear least squares adjustment + + Matrices + a partials wrt unknowns + b partials wrt to measured quantities + w constant portion of linearized equations evaluate using estimates of unknown parameters + p p covariance matrix of measured quantities + m b*p*transpose(b) which is the propogated covariance of the design equations + delta is vector of corrections to estimated unknowns + v residual vector + + linearized math model: + a*delta + b*v = w + + normal equation: + transpose(a)*inverse(m)*a*delta = transpose(a)*inverse(m)*w + + Solution: + delta = minverse(transpose(a)*inverse(m)*a)*transpose(a)*inverse(m)*w + iterate until corrections in delta are insignificant + + In this case to keep all resdiuals in pixel units and weight each observation identically + p is modeled as the identity matrix thus the solution can be built using a, b, m, and w + submatrices (indicated with a dot suffix) as follows + + Also in this particular case m is constant (proof given below) and thus un-necessary + + Normal equation: + sum(transpose(adot)*adot)*delta = sum(transpose(adot)*wdot) or ata*delta = atw + + solution: + delta = inverse(ata)*atw iterate until corrections in delta are insignificant + */ + + //clear out previous data + solV.clear(); + //matrix names loosely following the naming convention of Uotilia see large comment above + double wdot[2], adot[2][3], ata[6], atf[3]; + + int l, m, n, o, ni; //additonal indices for matrix multipilcations, etc + + for (j = 0; j < int(conFid.size()) - 1; j++) { + for (k = j + 1; k < int(conFid.size()); k++) { + //Starting by finding the two point solution: using the following close form algorithim + sol.theta = atan2(conFid[j][3] - conFid[k][3], conFid[j][2] - conFid[k][2]) - + atan2(conFid[j][1] - conFid[k][1], conFid[j][0] - conFid[k][0]); + sol.dx = conFid[j][2] - conFid[j][0] * cos(sol.theta) + conFid[j][1] * sin(sol.theta); + sol.dy = conFid[j][3] - conFid[j][0] * sin(sol.theta) - conFid[j][1] * cos(sol.theta); + + //generalize to include all points that are in general agreement with the solution + //NOTE: this is a Gauss-Helmert least squares optimization + //iteration limit is imposed here + // apon successful convergance the code breaks out of this for loop + for (ni = 0; ni < 50; ni++) { + //initailize normal equation + for (m = 0; m < 3; m++) + atf[m]=ata[m]=0.0; + for (; m < 6; m++) + ata[m] = 0.0; + + ct = cos(sol.theta); //cos and sin computed only once + st = sin(sol.theta); + + for (m = 0; m < int(conFid.size()); m++) { //for each point + wdot[0] = -(conFid[m][0] * ct - conFid[m][1] * st + sol.dx - conFid[m][2]); + wdot[1] = -(conFid[m][0] * st + conFid[m][1] * ct + sol.dy - conFid[m][3]); + + // the following calculation of mdot matrix is left in the code for ease of + // understanding, the matrix is constant (as shown below) + // and thus neither caculated nor stored in memory + //partials wrt measured quantities + //bdot[0][0] = ct; //partial of equation 1 wrt to x scani + //bdot[1][0] = st; //partial of equation 2 wrt to x scani + + //bdot[0][1] = -st; //partial of equation 1 wrt to y scani + //bdot[1][1] = ct; //partial of equation 2 wrt to y scani + + //bdot[0][2] = -1.0; //partial of equation 1 wrt to x scani+1 + //bdot[1][2] = 0.0; //partial of equation 2 wrt to x scani+1 + + //bdot[0][3] = 0.0; //partial of equation 1 wrt to y scani+1 + //bdot[1][3] = -1.0; //partial of equation 2 wrt to y scani+1 + + //covariance_matrix assumed identity--that keeps are weighted residuals in pixel + // units for easy interpretation + //mdot = bdot*covariance_matrix*transpose(bdot); + + //mdot[0] = bdot[0][0]*bodt[0][0] + + // bdot[0][1]*bodt[0][1] + + // bdot[0][2]*bodt[0][2] + + // bdot[0][3]*bodt[0][3] + // = ct*ct + st*st + 1 + 0 = 2 + //mdot[1] = bdot[0][0]*bodt[1][0] + + // bdot[0][1]*bodt[1][1] + + // bdot[0][2]*bodt[1][2] + + // bdot[0][3]*bodt[1][3] + // = ct*st - ct*st + 0 + 0 = 0 + //mdot[2] = bdot[1][0]*bodt[1][0] + + // bdot[1][1]*bodt[1][1] + + // bdot[1][2]*bodt[1][2] + + // bdot[1][3]*bodt[1][3] + // = st*st + ct*ct + 0 + 1 = 2 + + //mdot inverse + //mdot[0] = 0.5; + //mdot[1] = 0.0; + //mdot[2] = 0.5; + + if (sqrt(0.5 * (wdot[0] * wdot[0] + wdot[1] * wdot[1])) > 3.0) { + continue; //if the R^2 residual is greater than 3.0 pixels go on to the next point + } + + //partials wrt unknowns + adot[0][0] = -conFid[m][0]*st - conFid[m][1]*ct; //wrt to theta + adot[1][0] = conFid[m][0]*ct - conFid[m][1]*st; + + adot[0][1] = 1.0; //wrt dx + adot[1][1] = 0.0; + + adot[0][2] = 0.0; //wrt dy + adot[1][2] = 1.0; + + //build the normal equations + //add transpose(adot)*adot to ata + //note: becuase m is constant, and every weight is equal, it is irrelevant and ignored + for (l = 0; l < 3; l++) { + for (n = 0; n <= l; n++) { + for (o = 0; o < 2; o++) { + ata[isymp(l,n)] += adot[o][l] * adot[o][n]; + } + } + } + + //add transpose(adot)*wdot to atf + //note: becuase m is constant, and every weight is equal, it is irrelevant and ignored + atf[0] += adot[0][0] * wdot[0] + adot[1][0] * wdot[1]; + atf[1] += adot[0][1] * wdot[0] + adot[1][1] * wdot[1]; + atf[2] += adot[0][2] * wdot[0] + adot[1][2] * wdot[1]; + } + + //solve normal equations + //calculation is done in place as delta is returned in atf + if (choleski_solve(ata, atf, 3, 2) != 1) { + continue; //if the solution fails move on to the next two point solution + } + + //check for nonsense numbers + if (atf[2] != atf[2]) continue; + + //add corections + sol.theta += atf[0]; + sol.dx += atf[1]; + sol.dy += atf[2]; + + + if (fabs(atf[0]) < 1e-10 && fabs(atf[0]) < 1e-5 && fabs(atf[0]) < 1e-5) { + // solution converged + // find residual stats + // Note: residuals are caluclated for all points-including any that might have been + // excluded from the solution above + for (m = 0, sol.maxR = 0.0, sol.averR = 0.0, sol.limit[0] = 0.0, sol.limit[1] = 0.0; + m sol.maxR) { + sol.maxR = temp; + } + + sol.averR += temp; + + sol.limit[0] += conFid[m][0]; + sol.limit[1] += conFid[m][1]; + } + sol.averR /= double(conFid.size()); + + // save the solution + solV.push_back(sol); + break; // break out of iteration counting loop + } + } // iteration counting for loop + } + } // end of two point solutions loops + if (solV.size() < 1) { + string msg = "ERROR: unable to to calculate transformation coefficients. Check FILE_BASE and" + " MICRON parameters.\n"; + throw IException(IException::Io, msg, _FILEINFO_); + return; + } + + // now that all the two point solutions have been found go through and choose the best + // first lets thin by maximum residual + temp = solV[0].maxR; // find the lowest max residual + for (j = 1; j < int(solV.size()); j++) { + if (solV[j].maxR < temp) { + temp = solV[j].maxR; + } + } + // toss out any solution that has a maximum residual higher than ceil(lowest max residual) + temp = ceil(temp); + for (j = int(solV.size())-1; j >= 0; j--) { + if (solV[j].maxR > temp) { + solV.erase(solV.begin() + j); + } + } + + // then by average residual + temp = solV[0].averR; // find the lowest average residual + for (j = 1; j < int(solV.size()); j++) { + if (solV[j].averR < temp) { + temp = solV[j].averR; + } + } + // toss out any solutions that have average resiudals higher than ceil(lowest average residual) + temp = ceil(temp); + for (j = int(solV.size()) - 1; j >= 0; j--) { + if (solV[j].averR > temp) { + solV.erase(solV.begin() + j); + } + } + + // finally choose the transformation with the lowest theta angle as the final filter + k = 0; + for (j = 1; j < int(solV.size()); j++) { + if (fabs(solV[j].theta) < fabs(solV[k].theta)) { + k = j; + } + } + + // save solV[k] + trans[i].theta = solV[k].theta; + trans[i].dx = solV[k].dx; + trans[i].dy = solV[k].dy; + trans[i].averR = solV[k].averR; // average residual + trans[i].maxR = solV[k].maxR; // max resiudal + trans[i].limit[0] = solV[k].limit[0] / double(conFid.size()); + trans[i].limit[1] = solV[k].limit[1] / double(conFid.size()); + }// end of scan loop + + // translation from 8 to 8 + trans[7].theta = 0.0; + trans[7].dx = 0.0; + trans[7].dy = 0.0; + + // We now have seven transformations that convert from scan i to scan i+1 for i = 0 to 7 + // as scan 8 is the left most scan and has a nominally identity transform from its sample-line + // system to the stitched cube sample-line system. Combine transformations so that each transforms + // is from scan i to scan 8 + for (i = 0; i < 7; i++) { + for (j = i + 1; j < 8; j++) { + // Make i equal to the combined transformation of i and j + trans[i].theta = atan2( cos(trans[j].theta)*sin(trans[i].theta) + + sin(trans[j].theta)*cos(trans[i].theta), + cos(trans[j].theta)*cos(trans[i].theta) - + sin(trans[j].theta)*sin(trans[i].theta)); + + temp = cos(trans[j].theta)*trans[i].dx - sin(trans[j].theta)*trans[i].dy + + trans[j].dx; + trans[i].dy = sin(trans[j].theta)*trans[i].dx + cos(trans[j].theta)*trans[i].dy + + trans[j].dy; + + trans[i].dx = temp; + } + } + + + + // Now lets find the extents of the stitched image + double minS = 1, + maxS = panC[7]->sampleCount(), + minL = 1, + maxL = panC[7]->lineCount(); + + for (i = 0; i < 7; i++) { + scanS = panC[i]->sampleCount(); + scanL = panC[i]->lineCount(); + + // Convert the four corner points to the scan 8 domain and determine the greatest extents + temp = cos(trans[i].theta) - sin(trans[i].theta) + trans[i].dx; + if (temp < minS) minS = temp; + if (temp > maxS) maxS = temp; + temp = cos(trans[i].theta)*scanS - sin(trans[i].theta) + trans[i].dx; + if (temp < minS) minS = temp; + if (temp > maxS) maxS = temp; + temp = cos(trans[i].theta)*scanS - sin(trans[i].theta)*scanL + trans[i].dx; + if (temp < minS) minS = temp; + if (temp > maxS) maxS = temp; + temp = cos(trans[i].theta) - sin(trans[i].theta)*scanL + trans[i].dx; + if (temp < minS) minS = temp; + if (temp > maxS) maxS = temp; + + temp = sin(trans[i].theta) + cos(trans[i].theta) + trans[i].dy; + if (temp < minL) minL = temp; + if (temp > maxL) maxL = temp; + temp = sin(trans[i].theta)*scanS + cos(trans[i].theta) + trans[i].dy; + if (temp < minL) minL = temp; + if (temp > maxL) maxL = temp; + temp = sin(trans[i].theta)*scanS + cos(trans[i].theta)*scanL + trans[i].dy; + if (temp < minL) minL = temp; + if (temp > maxL) maxL = temp; + temp = sin(trans[i].theta) + cos(trans[i].theta)*scanL + trans[i].dy; + if (temp < minL) minL = temp; + if (temp > maxL) maxL = temp; + } + + // Update the transformations to make the minimum line = 1 + for (i = 0; i < 8; i++) { + trans[i].dy += 1- minL; + } + + maxL += ceil(1 - minL); + minL = 1; + + // Update the transformations to make the minimum sample = 1 + for (i = 0; i < 8; i++) { + trans[i].dx += 1- minS; + } + + maxS += ceil(1 - minS); + minS = 1; + + // trans[i].limit was previous calculated as the center of mass of the conjugate fiducials in the + // scani system + // Transform it now into the stitched cube system to be used as a limit between scans + for (i = 0; i < 7; i++) { + temp = trans[i].limit[0]*cos(trans[i].theta) - + trans[i].limit[1]*sin(trans[i].theta) + trans[i].dx; + trans[i].limit[1] = trans[i].limit[0]*sin(trans[i].theta) + + trans[i].limit[1]*cos(trans[i].theta) + trans[i].dy; + trans[i].limit[0] = floor(temp); + } + + + + // Finally shift and invert the transformations so that they convert from the rubber sheeted + // sub-cubes back to the sub-scans. The result of this will be transformations that will convert + // from a set of coordinate systems that differ only by an integral transformation back to + // sub-scans, each sub-scan will be rubber-sheeted into one of these systems, and the limit[0] + // values record where they belong in the stiched image/system. + for (i = 0; i < 8; i++) { + if (i == 7) { + sampleFrom = 1.0; + } + else { + sampleFrom = trans[i].limit[0]; + } + //a shift so that the transform puts sampleFrom at sample 1 of the rubber sheeted sub cubes + trans[i].dx += -sampleFrom+1; + trans[i].theta = -trans[i].theta; + temp = -cos(trans[i].theta) * trans[i].dx + sin(trans[i].theta) * trans[i].dy; + trans[i].dy = -sin(trans[i].theta) * trans[i].dx - cos(trans[i].theta) * trans[i].dy; + trans[i].dx = temp; + } + + // Make the final maxes integral + maxS = ceil(maxS); + maxL = ceil(maxL); + + // Attributes of input and ouput cubes to process class + CubeAttributeOutput att; + CubeAttributeInput attI; + + attI.setAttributes("some.cub+1"); //will only be processing one band from the input + + //make output attributes match the input + att.setFileFormat( panC[0]->format() ); + att.setByteOrder( panC[0]->byteOrder() ); + att.setPixelType( panC[0]->pixelType() ); + if (panC[0]->labelsAttached()) { + att.setLabelAttachment(AttachedLabel); + } + else { + att.setLabelAttachment(DetachedLabel); + } + + //define an output cube + outputC.setDimensions(int(maxS), int(maxL), 1); + outputC.setPixelType(panC[0]->pixelType()); //set pixel type + tempString = ui.GetFileName("TO"); + outputC.create(tempString); + outputC.close(); //closing the output cube so that it can be opened by the mosaic process + ProcessMosaic mosaic; + mosaic.SetOutputCube("TO", ui); + mosaic.SetBandBinMatch(false); + + + //transform and mosaic the content from each scan + for (i = 0; i < 8; i++) { //for each scan + FileName tempFile = FileName::createTempFile("$temporary/tempscan.cub"); + + scanS = panC[i]->sampleCount(); + scanL = panC[i]->lineCount(); + + //define the sample range + if (i == 0) { + sampleTo = maxS; + } + else { + sampleTo = trans[i-1].limit[0]; + } + if (i == 7) { + sampleFrom = 1.0; + } + else { + sampleFrom = trans[i].limit[0]; + } + + Interpolator bilinearInt(Interpolator::BiLinearType); + ProcessRubberSheet rubberS; + + //use ProcessRubber sheet to create the sub-cube for this scanz + Trans2d3p transform(trans[i].theta, trans[i].dx, trans[i].dy, + int(sampleTo - sampleFrom), + int(maxL)); + rubberS.SetInputCube(panC[i]); + rubberS.SetOutputCube(tempFile.expanded(), + att, + transform.OutputSamples(), + transform.OutputLines(), + 1); + rubberS.Progress()->SetText( + QObject::tr("Transforming Scan%1: ").arg(i+1).toStdString().c_str()); + rubberS.StartProcess(transform, bilinearInt); + //end process closes the cubes and deletes the dynamically allocated memory (panC[i]) + rubberS.EndProcess(); + + //use process mosaic to add the sub cube to the stiched cube + mosaic.SetInputCube(tempFile.expanded(), + attI, + 1, + 1, + 1, + transform.OutputSamples(), + transform.OutputLines(), + 1); + rubberS.Progress()->SetText(QObject::tr("Mosaicing Scan%1: ").arg(i+1).toStdString().c_str()); + mosaic.StartProcess(int(sampleFrom), 1, 1); + + //clear input cube + mosaic.ClearInputCubes(); + + QFile::remove(tempFile.expanded()); //delete temporary cube if it exists + } + mosaic.EndProcess(); + } +} diff --git a/isis/src/apollo/apps/apollopanstitcher/apollopanstitcher.h b/isis/src/apollo/apps/apollopanstitcher/apollopanstitcher.h new file mode 100644 index 0000000000..4e6bb91107 --- /dev/null +++ b/isis/src/apollo/apps/apollopanstitcher/apollopanstitcher.h @@ -0,0 +1,11 @@ +#ifndef apollopanstitcher_h +#define apollopanstitcher_h + +#include "Pvl.h" +#include "UserInterface.h" + +namespace Isis{ + extern void apolloPanStitcher(UserInterface &ui); +} + +#endif diff --git a/isis/src/apollo/apps/apollopanstitcher/main.cpp b/isis/src/apollo/apps/apollopanstitcher/main.cpp index 65c7c077bd..cbf3b429b3 100644 --- a/isis/src/apollo/apps/apollopanstitcher/main.cpp +++ b/isis/src/apollo/apps/apollopanstitcher/main.cpp @@ -1,794 +1,14 @@ -/* -Program uses the fiducial marks along the edges of Apollo Panoramic Images to stich eight sub-scans into a single conitinuous image - -*/ - -//Isis.h and IsisDebug.h if needed #include "Isis.h" -//C++ standard libraries if needed -#include -#include - -//QT libraries if needed if needed -#include - -//third party libraries if needed -#include - -//Isis Headers if needed -#include "AutoReg.h" -#include "AutoRegFactory.h" -#include "Centroid.h" -#include "CentroidApolloPan.h" -#include "Chip.h" -#include "Cube.h" -#include "CubeAttribute.h" -#include "Interpolator.h" -#include "ProcessMosaic.h" -#include "ProcessRubberSheet.h" -#include "Pvl.h" -#include "PvlGroup.h" -#include "PvlKeyword.h" -#include "PvlObject.h" -#include "Table.h" -#include "Trans2d3p.h" -#include "UserInterface.h" - - - - -/** - * - * @internal - * @History 2011-11-02 Orrin Thomas - Original Version. - * - * @History 2017-06-29 Christopher Combs - Moved initialization of trans[7] to before - * the loop that uses it. Fixes #4948. - * - */ - +#include "apollopanstitcher.h" +#include "Application.h" +#include "Pvl.h" using namespace std; using namespace Isis; -//constants for 5 micro resolution scans -#define SCALE 10.0 //scale used for down sizingS patternS an search chips -#define SEARCHh 1400.0 //number of lines (in 5-micron-pixels) in search space for the first - // fiducial -#define SEARCHc 350.0 //number of samples per edge(in 5-micron-pixels) in each sub-search area -#define AVERs 5286.0 //average smaples (in 5-micron-pixels) between fiducials -#define AVERl 23459.0 //average diference (in 5-micron-pixels) between the top and bottom - // fiducials -#define TRANS_N 28520.0 //nomimal dx between scan lines, scani x + 28532 = - // (approx) scani+1 x - // also the size of the search area for the first fiducial - -double R_MOON[3]; - -typedef struct TRANS2d3p { //2d 3 parameter transformation - double theta; - double dx; - double dy; - double averR; //average residual - double maxR; //max resiudal - double limit[2]; //where to stop using this transformation, will eventually be a point in the - //coordinate system of the new cube, a vertical line through this point will - //be the dividing line -}TRANS2d3p; - void IsisMain() { - Chip patternS,searchS; //scaled pattern and search chips - Cube *panC[8], //panoramic image - outputC, //outputC - fidC; //Fiducial image - int regStatus, - scanS, //samples in a sub scan - scanL, //lines in a sub scan - refL, //number of lines in the patternS - refS; //number of samples in the patternS - - Pvl pvl; - - bool foundFirst; - - int i,j,k; - - double sampleFrom, sampleTo; //limits for rubber sheeting and mosaicing the sub-scans - double play; //how close a pixel must be to be considered on the ellipse - - //vector of possible solutions to a transformation-all transformations resulting from any - // possible combination of minimal sets of data are calculated to find the best solution - // in an outlier resistant way - std::vector solV; - TRANS2d3p sol; //individual solution to be add to the above vector - - TRANS2d3p trans[8]; //final solutions for the transformations-one for each sub-scan - - QString fileName,tempIString; - QString fileBaseName, tempString; UserInterface &ui = Application::GetUserInterface(); - - double l = 1, s = 1, sample, line, temp,//line and sample coordinates for looping through the panC - ct, //"cosine of theta" used to prevent evaluted the trig function thousands of times - st; //"sine of theta" used to prevent evaluted the trig function thousands of times - - std::vector > scanFid; //fiducials of Scan i - std::vector > conFid; //fiducials of Scan i - std::vector fid(2); //individual fiducial measurement - std::vector pairFid(4); //pair of conjugate fiducial coordinates - - //read the image resolutions and scale the constants acordingly - double resolution = ui.GetDouble("MICRONS"), //pixel size in microns - scale = SCALE *5.0/resolution, //reduction scale for quicker autoregistrations - searchHeight = SEARCHh*5.0/resolution, //number of lines (in 5-micron-pixels) in - // search space for the first fiducial - searchCellSize = SEARCHc*5.0/resolution, //height/width of search chips block - averageSamples = AVERs *5.0/resolution, //scaled smaples between fiducials - nominalTrans = TRANS_N*5.0/resolution, //nominal sample direction overlap between - // successive sub scans - averageLines = AVERl *5.0/resolution; //scaled average diference between the top - // and bottom fiducials - - if (15.0 / resolution < 1.5) { - play = 1.5; - } - else { - play = 15.0 / resolution; - } - - //copy the patternS chip (the entire ApolloPanFiducialMark.cub) - FileName fiducialFileName("$apollo15/calibration/ApolloPanFiducialMark.cub"); - fidC.open(fiducialFileName.expanded(),"r"); - if (!fidC.isOpen()) { - string msg = "Unable to open the fiducial patternS cube: ApolloPanFiducialMark.cub\n"; - throw IException(IException::User, msg, _FILEINFO_); - } - refL = fidC.lineCount(); - refS = fidC.sampleCount(); - //scaled pattern chip for fast matching - patternS.SetSize((int)floor((refS-2)/SCALE), (int)floor((refL-2)/SCALE)); - patternS.TackCube((refS-1)/2, (refL-1)/2); - patternS.Load(fidC, 0, SCALE); - - //make sure the eight "FROM" cubes exist, and can can be opened for reading - fileBaseName = ui.GetString("FILE_BASE"); - for (i = 1; i <= 8; i++) { - panC[i-1] = new Cube; - fileName = fileBaseName + "-000" + toString(i) + ".cub"; - panC[i-1]->open(fileName, "r"); - if (!panC[i-1]->isOpen()) { - string msg = "Unable to open input cube: " + IString(fileName) + "\n"; - throw IException(IException::User, msg, _FILEINFO_); - } - } - - //parameters for maximum correlation autoregestration - // see: ///usgs/pkgs/isis3nightly2011-09-21/isis/doc/documents/patternSMatch/patternSMatch.html#DistanceTolerance - FileName fiducialPvl("$ISISROOT/appdata/templates/apollo/PanFiducialFinder.def"); - pvl.read(fiducialPvl.expanded()); //read in the autoreg parameters - AutoReg *arS = AutoRegFactory::Create(pvl); - - *arS->PatternChip() = patternS; //patternS chip is constant - - //set up a centroid measurer - //Centroid centroid; - CentroidApolloPan centroid(resolution); - Chip inputChip,selectionChip; - if (panC[0]->pixelType() == 1) { //UnsignedByte - centroid.setDNRange(12, 1e99); //8 bit bright target - } - else { - centroid.setDNRange(3500, 1e99); //16 bit bright target - } - inputChip.SetSize(int(ceil(200 * 5.0 / resolution)), int(ceil(200 * 5.0 / resolution))); - - //find conjugate finducials and calcualte transformations - for (i = 0; i < 7; i++) { - //for each scan segment 0-6 - //1. find the probable conjugate fiducials (those with the smallest sample coordinates) - //2. use the nominal transformation to find their matches in the next scan - //3. robustly calculate the transformation (that is in a manner resistent to erroneous - // measurements) - - //Step 1: find the probable conjugate fiducials (those with the smallest sample coordinates) - scanS = panC[i+1]->sampleCount(); - - scanFid.clear(); - - //search sizes are constanst - searchS.SetSize(int(searchCellSize/scale),int(searchCellSize/scale)); - //now start searching along a horizontal line for the first fiducial mark - foundFirst = false; - for (l = searchCellSize / 2; - l < searchHeight + searchCellSize / 2.0 && !foundFirst; - l = l + (searchCellSize - 125 * 5.0 / resolution) ) { - for (s = searchCellSize / 2; - s < (nominalTrans + searchCellSize) / 2.0 && !foundFirst; - s = s + (searchCellSize - 125 * 5.0 / resolution) ) { - searchS.TackCube(s, l); - searchS.Load(*panC[i], 0, scale); - *arS->SearchChip() = searchS; - regStatus = arS->Register(); - if (regStatus == AutoReg::SuccessPixel) { - inputChip.TackCube(arS->CubeSample(), arS->CubeLine()); - inputChip.Load(*panC[i], 0, 1); - inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); - //continuous dynamic range selection - centroid.selectAdaptive(&inputChip, &selectionChip); - //selectionChip.Write("firstselectionTemp.cub"); - //elliptical trimming/smoothing - if (centroid.elipticalReduction(&selectionChip, 95, play, 2000)) { - //center of mass to reduce selection to a single measure - centroid.centerOfMass(&selectionChip, &sample, &line); - inputChip.SetChipPosition(sample, line); - fid[0] = inputChip.CubeSample(); - fid[1] = inputChip.CubeLine(); - scanFid.push_back(fid); - foundFirst = true; //once the first fiducial is found stop - } - } - } - } - if (s >= averageLines + searchCellSize / 2.0) { - QString msg = "Unable to locate a fiducial mark in the input cube [" + panC[i]->fileName() - + "]. Check FROM and MICRONS parameters."; - throw IException(IException::Io, msg, _FILEINFO_); - return; - } - - - for (s = scanFid[0][0], l = scanFid[0][1]; s + nominalTrans < scanS; s+= averageSamples / 2.0) { - //look for the bottom fiducial - searchS.TackCube(s, l + averageLines); - searchS.Load(*panC[i], 0, scale); - *arS->SearchChip() = searchS; - regStatus = arS->Register(); - if (regStatus == AutoReg::SuccessPixel) { - inputChip.TackCube( arS->CubeSample(), arS->CubeLine()); - inputChip.Load(*panC[i], 0, 1); - inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); - } - else { - inputChip.TackCube(s, l + averageLines); - inputChip.Load(*panC[i], 0, 1); - inputChip.SetCubePosition(s, l + averageLines); - } - //continuous dynamic range selection - centroid.selectAdaptive(&inputChip, &selectionChip); - //elliptical trimming/smoothing... if this fails move on - if (centroid.elipticalReduction(&selectionChip, 95, play, 2000) != 0) { - //center of mass to reduce selection to a single measure - centroid.centerOfMass(&selectionChip,&sample,&line); - inputChip.SetChipPosition(sample,line); - fid[0] = inputChip.CubeSample(); - fid[1] = inputChip.CubeLine(); - scanFid.push_back(fid); - } - - //look for the top fiducial - //is this is the first time through the loop? - if (s == scanFid[0][0]) { - continue; //then the top fiducial was already found - } - searchS.TackCube(s, l); - searchS.Load(*panC[i], 0, scale); - *arS->SearchChip() = searchS; - regStatus = arS->Register(); - if (regStatus == AutoReg::SuccessPixel) { - inputChip.TackCube(arS->CubeSample(), arS->CubeLine()); - inputChip.Load(*panC[i], 0, 1); - inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); - } - else { - inputChip.TackCube(s, l + averageLines); - inputChip.Load(*panC[i], 0, 1); - inputChip.SetCubePosition(s, l + averageLines); - } - //continuous dynamic range selection - centroid.selectAdaptive(&inputChip, &selectionChip); - //elliptical trimming/smoothing... if this fails move on - if (centroid.elipticalReduction(&selectionChip, 95, play, 2000) != 0) { - //center of mass to reduce selection to a single measure - centroid.centerOfMass(&selectionChip, &sample, &line); - inputChip.SetChipPosition(sample, line); - //s and l are refined here to help follow any trends in the image/film - s=fid[0] = inputChip.CubeSample(); - l=fid[1] = inputChip.CubeLine(); - scanFid.push_back(fid); - } - } - - // STEP 2. use the nominal transformation to find their matches in the next scan - conFid.clear(); - for (j = scanFid.size() - 1; j >= 0; j--) { - searchS.TackCube(scanFid[j][0] + nominalTrans,scanFid[j][1]); - searchS.Load(*panC[i+1], 0, scale); - *arS->SearchChip() = searchS; - regStatus = arS->Register(); - //if autoreg is successful a smaller window can be used for the following steps - if (regStatus == AutoReg::SuccessPixel) { - inputChip.SetSize(int(200 * 5.0 / resolution), int(150 * 5.0 / resolution)); - inputChip.TackCube(arS->CubeSample(), arS->CubeLine()); - inputChip.Load(*panC[i+1], 0, 1); - inputChip.SetCubePosition(arS->CubeSample(), arS->CubeLine()); - } - else { //if autoreg is unsuccessful, a larger window will be used - inputChip.SetSize(int(ceil(200 * 5.0 / resolution)), int(ceil(200 * 5.0 / resolution))); - inputChip.TackCube(scanFid[j][0] + nominalTrans, scanFid[j][1]); - inputChip.Load(*panC[i+1], 0, 1); - inputChip.SetCubePosition(scanFid[j][0] + nominalTrans, scanFid[j][1]); - } - //continuous dynamic range selection - centroid.selectAdaptive(&inputChip, &selectionChip); - //elliptical trimming/smoothing... if this fails move on - if (centroid.elipticalReduction(&selectionChip, 95, play, 2000) != 0) { - //center of mass to reduce selection to a single measure - centroid.centerOfMass(&selectionChip, &sample, &line); - inputChip.SetChipPosition(sample, line); - //now that the complete conjugate pair has been located save it - pairFid[2] = inputChip.CubeSample(); - pairFid[3] = inputChip.CubeLine(); - pairFid[0] = scanFid[j][0]; - pairFid[1] = scanFid[j][1]; - conFid.push_back(pairFid); - } - } - //STEP 3. robustly calculate the transformation (that is in a manner resistent to erroneous - // measurements) - //calculate the transformation from i+1 to i - //RANSAC2 algorithim - //Find two point solutions and then attempt to generalize to more points - //this helps to identify measurement blunders by seeing if certain points are excluded - //consistently--or equivalently if certain solutions keep repeating themselves - //the limitation in this case is that there will never be more than 6 points, if the number - //of points drops to <=three blunder detection becomes impossible - - /*This is a Gauss-Helemert non-linear least squares adjustment - - Matrices - a partials wrt unknowns - b partials wrt to measured quantities - w constant portion of linearized equations evaluate using estimates of unknown parameters - p p covariance matrix of measured quantities - m b*p*transpose(b) which is the propogated covariance of the design equations - delta is vector of corrections to estimated unknowns - v residual vector - - linearized math model: - a*delta + b*v = w - - normal equation: - transpose(a)*inverse(m)*a*delta = transpose(a)*inverse(m)*w - - Solution: - delta = minverse(transpose(a)*inverse(m)*a)*transpose(a)*inverse(m)*w - iterate until corrections in delta are insignificant - - In this case to keep all resdiuals in pixel units and weight each observation identically - p is modeled as the identity matrix thus the solution can be built using a, b, m, and w - submatrices (indicated with a dot suffix) as follows - - Also in this particular case m is constant (proof given below) and thus un-necessary - - Normal equation: - sum(transpose(adot)*adot)*delta = sum(transpose(adot)*wdot) or ata*delta = atw - - solution: - delta = inverse(ata)*atw iterate until corrections in delta are insignificant - */ - - //clear out previous data - solV.clear(); - //matrix names loosely following the naming convention of Uotilia see large comment above - double wdot[2], adot[2][3], ata[6], atf[3]; - - int l, m, n, o, ni; //additonal indices for matrix multipilcations, etc - - for (j = 0; j < int(conFid.size()) - 1; j++) { - for (k = j + 1; k < int(conFid.size()); k++) { - //Starting by finding the two point solution: using the following close form algorithim - sol.theta = atan2(conFid[j][3] - conFid[k][3], conFid[j][2] - conFid[k][2]) - - atan2(conFid[j][1] - conFid[k][1], conFid[j][0] - conFid[k][0]); - sol.dx = conFid[j][2] - conFid[j][0] * cos(sol.theta) + conFid[j][1] * sin(sol.theta); - sol.dy = conFid[j][3] - conFid[j][0] * sin(sol.theta) - conFid[j][1] * cos(sol.theta); - - //generalize to include all points that are in general agreement with the solution - //NOTE: this is a Gauss-Helmert least squares optimization - //iteration limit is imposed here - // apon successful convergance the code breaks out of this for loop - for (ni = 0; ni < 50; ni++) { - //initailize normal equation - for (m = 0; m < 3; m++) - atf[m]=ata[m]=0.0; - for (; m < 6; m++) - ata[m] = 0.0; - - ct = cos(sol.theta); //cos and sin computed only once - st = sin(sol.theta); - - for (m = 0; m < int(conFid.size()); m++) { //for each point - wdot[0] = -(conFid[m][0] * ct - conFid[m][1] * st + sol.dx - conFid[m][2]); - wdot[1] = -(conFid[m][0] * st + conFid[m][1] * ct + sol.dy - conFid[m][3]); - - // the following calculation of mdot matrix is left in the code for ease of - // understanding, the matrix is constant (as shown below) - // and thus neither caculated nor stored in memory - //partials wrt measured quantities - //bdot[0][0] = ct; //partial of equation 1 wrt to x scani - //bdot[1][0] = st; //partial of equation 2 wrt to x scani - - //bdot[0][1] = -st; //partial of equation 1 wrt to y scani - //bdot[1][1] = ct; //partial of equation 2 wrt to y scani - - //bdot[0][2] = -1.0; //partial of equation 1 wrt to x scani+1 - //bdot[1][2] = 0.0; //partial of equation 2 wrt to x scani+1 - - //bdot[0][3] = 0.0; //partial of equation 1 wrt to y scani+1 - //bdot[1][3] = -1.0; //partial of equation 2 wrt to y scani+1 - - //covariance_matrix assumed identity--that keeps are weighted residuals in pixel - // units for easy interpretation - //mdot = bdot*covariance_matrix*transpose(bdot); - - //mdot[0] = bdot[0][0]*bodt[0][0] + - // bdot[0][1]*bodt[0][1] + - // bdot[0][2]*bodt[0][2] + - // bdot[0][3]*bodt[0][3] - // = ct*ct + st*st + 1 + 0 = 2 - //mdot[1] = bdot[0][0]*bodt[1][0] + - // bdot[0][1]*bodt[1][1] + - // bdot[0][2]*bodt[1][2] + - // bdot[0][3]*bodt[1][3] - // = ct*st - ct*st + 0 + 0 = 0 - //mdot[2] = bdot[1][0]*bodt[1][0] + - // bdot[1][1]*bodt[1][1] + - // bdot[1][2]*bodt[1][2] + - // bdot[1][3]*bodt[1][3] - // = st*st + ct*ct + 0 + 1 = 2 - - //mdot inverse - //mdot[0] = 0.5; - //mdot[1] = 0.0; - //mdot[2] = 0.5; - - if (sqrt(0.5 * (wdot[0] * wdot[0] + wdot[1] * wdot[1])) > 3.0) { - continue; //if the R^2 residual is greater than 3.0 pixels go on to the next point - } - - //partials wrt unknowns - adot[0][0] = -conFid[m][0]*st - conFid[m][1]*ct; //wrt to theta - adot[1][0] = conFid[m][0]*ct - conFid[m][1]*st; - - adot[0][1] = 1.0; //wrt dx - adot[1][1] = 0.0; - - adot[0][2] = 0.0; //wrt dy - adot[1][2] = 1.0; - - //build the normal equations - //add transpose(adot)*adot to ata - //note: becuase m is constant, and every weight is equal, it is irrelevant and ignored - for (l = 0; l < 3; l++) { - for (n = 0; n <= l; n++) { - for (o = 0; o < 2; o++) { - ata[isymp(l,n)] += adot[o][l] * adot[o][n]; - } - } - } - - //add transpose(adot)*wdot to atf - //note: becuase m is constant, and every weight is equal, it is irrelevant and ignored - atf[0] += adot[0][0] * wdot[0] + adot[1][0] * wdot[1]; - atf[1] += adot[0][1] * wdot[0] + adot[1][1] * wdot[1]; - atf[2] += adot[0][2] * wdot[0] + adot[1][2] * wdot[1]; - } - - //solve normal equations - //calculation is done in place as delta is returned in atf - if (choleski_solve(ata, atf, 3, 2) != 1) { - continue; //if the solution fails move on to the next two point solution - } - - //check for nonsense numbers - if (atf[2] != atf[2]) continue; - - //add corections - sol.theta += atf[0]; - sol.dx += atf[1]; - sol.dy += atf[2]; - - - if (fabs(atf[0]) < 1e-10 && fabs(atf[0]) < 1e-5 && fabs(atf[0]) < 1e-5) { - // solution converged - // find residual stats - // Note: residuals are caluclated for all points-including any that might have been - // excluded from the solution above - for (m = 0, sol.maxR = 0.0, sol.averR = 0.0, sol.limit[0] = 0.0, sol.limit[1] = 0.0; - m sol.maxR) { - sol.maxR = temp; - } - - sol.averR += temp; - - sol.limit[0] += conFid[m][0]; - sol.limit[1] += conFid[m][1]; - } - sol.averR /= double(conFid.size()); - - // save the solution - solV.push_back(sol); - break; // break out of iteration counting loop - } - } // iteration counting for loop - } - } // end of two point solutions loops - if (solV.size() < 1) { - string msg = "ERROR: unable to to calculate transformation coefficients. Check FILE_BASE and" - " MICRON parameters.\n"; - throw IException(IException::Io, msg, _FILEINFO_); - return; - } - - // now that all the two point solutions have been found go through and choose the best - // first lets thin by maximum residual - temp = solV[0].maxR; // find the lowest max residual - for (j = 1; j < int(solV.size()); j++) { - if (solV[j].maxR < temp) { - temp = solV[j].maxR; - } - } - // toss out any solution that has a maximum residual higher than ceil(lowest max residual) - temp = ceil(temp); - for (j = int(solV.size())-1; j >= 0; j--) { - if (solV[j].maxR > temp) { - solV.erase(solV.begin() + j); - } - } - - // then by average residual - temp = solV[0].averR; // find the lowest average residual - for (j = 1; j < int(solV.size()); j++) { - if (solV[j].averR < temp) { - temp = solV[j].averR; - } - } - // toss out any solutions that have average resiudals higher than ceil(lowest average residual) - temp = ceil(temp); - for (j = int(solV.size()) - 1; j >= 0; j--) { - if (solV[j].averR > temp) { - solV.erase(solV.begin() + j); - } - } - - // finally choose the transformation with the lowest theta angle as the final filter - k = 0; - for (j = 1; j < int(solV.size()); j++) { - if (fabs(solV[j].theta) < fabs(solV[k].theta)) { - k = j; - } - } - - // save solV[k] - trans[i].theta = solV[k].theta; - trans[i].dx = solV[k].dx; - trans[i].dy = solV[k].dy; - trans[i].averR = solV[k].averR; // average residual - trans[i].maxR = solV[k].maxR; // max resiudal - trans[i].limit[0] = solV[k].limit[0] / double(conFid.size()); - trans[i].limit[1] = solV[k].limit[1] / double(conFid.size()); - }// end of scan loop - - // translation from 8 to 8 - trans[7].theta = 0.0; - trans[7].dx = 0.0; - trans[7].dy = 0.0; - - // We now have seven transformations that convert from scan i to scan i+1 for i = 0 to 7 - // as scan 8 is the left most scan and has a nominally identity transform from its sample-line - // system to the stitched cube sample-line system. Combine transformations so that each transforms - // is from scan i to scan 8 - for (i = 0; i < 7; i++) { - for (j = i + 1; j < 8; j++) { - // Make i equal to the combined transformation of i and j - trans[i].theta = atan2( cos(trans[j].theta)*sin(trans[i].theta) + - sin(trans[j].theta)*cos(trans[i].theta), - cos(trans[j].theta)*cos(trans[i].theta) - - sin(trans[j].theta)*sin(trans[i].theta)); - - temp = cos(trans[j].theta)*trans[i].dx - sin(trans[j].theta)*trans[i].dy - + trans[j].dx; - trans[i].dy = sin(trans[j].theta)*trans[i].dx + cos(trans[j].theta)*trans[i].dy - + trans[j].dy; - - trans[i].dx = temp; - } - } - - - - // Now lets find the extents of the stitched image - double minS = 1, - maxS = panC[7]->sampleCount(), - minL = 1, - maxL = panC[7]->lineCount(); - - for (i = 0; i < 7; i++) { - scanS = panC[i]->sampleCount(); - scanL = panC[i]->lineCount(); - - // Convert the four corner points to the scan 8 domain and determine the greatest extents - temp = cos(trans[i].theta) - sin(trans[i].theta) + trans[i].dx; - if (temp < minS) minS = temp; - if (temp > maxS) maxS = temp; - temp = cos(trans[i].theta)*scanS - sin(trans[i].theta) + trans[i].dx; - if (temp < minS) minS = temp; - if (temp > maxS) maxS = temp; - temp = cos(trans[i].theta)*scanS - sin(trans[i].theta)*scanL + trans[i].dx; - if (temp < minS) minS = temp; - if (temp > maxS) maxS = temp; - temp = cos(trans[i].theta) - sin(trans[i].theta)*scanL + trans[i].dx; - if (temp < minS) minS = temp; - if (temp > maxS) maxS = temp; - - temp = sin(trans[i].theta) + cos(trans[i].theta) + trans[i].dy; - if (temp < minL) minL = temp; - if (temp > maxL) maxL = temp; - temp = sin(trans[i].theta)*scanS + cos(trans[i].theta) + trans[i].dy; - if (temp < minL) minL = temp; - if (temp > maxL) maxL = temp; - temp = sin(trans[i].theta)*scanS + cos(trans[i].theta)*scanL + trans[i].dy; - if (temp < minL) minL = temp; - if (temp > maxL) maxL = temp; - temp = sin(trans[i].theta) + cos(trans[i].theta)*scanL + trans[i].dy; - if (temp < minL) minL = temp; - if (temp > maxL) maxL = temp; - } - - // Update the transformations to make the minimum line = 1 - for (i = 0; i < 8; i++) { - trans[i].dy += 1- minL; - } - - maxL += ceil(1 - minL); - minL = 1; - - // Update the transformations to make the minimum sample = 1 - for (i = 0; i < 8; i++) { - trans[i].dx += 1- minS; - } - - maxS += ceil(1 - minS); - minS = 1; - - // trans[i].limit was previous calculated as the center of mass of the conjugate fiducials in the - // scani system - // Transform it now into the stitched cube system to be used as a limit between scans - for (i = 0; i < 7; i++) { - temp = trans[i].limit[0]*cos(trans[i].theta) - - trans[i].limit[1]*sin(trans[i].theta) + trans[i].dx; - trans[i].limit[1] = trans[i].limit[0]*sin(trans[i].theta) + - trans[i].limit[1]*cos(trans[i].theta) + trans[i].dy; - trans[i].limit[0] = floor(temp); - } - - - - // Finally shift and invert the transformations so that they convert from the rubber sheeted - // sub-cubes back to the sub-scans. The result of this will be transformations that will convert - // from a set of coordinate systems that differ only by an integral transformation back to - // sub-scans, each sub-scan will be rubber-sheeted into one of these systems, and the limit[0] - // values record where they belong in the stiched image/system. - for (i = 0; i < 8; i++) { - if (i == 7) { - sampleFrom = 1.0; - } - else { - sampleFrom = trans[i].limit[0]; - } - //a shift so that the transform puts sampleFrom at sample 1 of the rubber sheeted sub cubes - trans[i].dx += -sampleFrom+1; - trans[i].theta = -trans[i].theta; - temp = -cos(trans[i].theta) * trans[i].dx + sin(trans[i].theta) * trans[i].dy; - trans[i].dy = -sin(trans[i].theta) * trans[i].dx - cos(trans[i].theta) * trans[i].dy; - trans[i].dx = temp; - } - - // Make the final maxes integral - maxS = ceil(maxS); - maxL = ceil(maxL); - - // Attributes of input and ouput cubes to process class - CubeAttributeOutput att; - CubeAttributeInput attI; - - attI.setAttributes("some.cub+1"); //will only be processing one band from the input - - //make output attributes match the input - att.setFileFormat( panC[0]->format() ); - att.setByteOrder( panC[0]->byteOrder() ); - att.setPixelType( panC[0]->pixelType() ); - if (panC[0]->labelsAttached()) { - att.setLabelAttachment(AttachedLabel); - } - else { - att.setLabelAttachment(DetachedLabel); - } - - //define an output cube - outputC.setDimensions(int(maxS), int(maxL), 1); - outputC.setPixelType(panC[0]->pixelType()); //set pixel type - tempString = ui.GetFileName("TO"); - outputC.create(tempString); - outputC.close(); //closing the output cube so that it can be opened by the mosaic process - ProcessMosaic mosaic; - mosaic.SetOutputCube("TO"); - mosaic.SetBandBinMatch(false); - - - //transform and mosaic the content from each scan - for (i = 0; i < 8; i++) { //for each scan - FileName tempFile = FileName::createTempFile("$temporary/tempscan.cub"); - - scanS = panC[i]->sampleCount(); - scanL = panC[i]->lineCount(); - - //define the sample range - if (i == 0) { - sampleTo = maxS; - } - else { - sampleTo = trans[i-1].limit[0]; - } - if (i == 7) { - sampleFrom = 1.0; - } - else { - sampleFrom = trans[i].limit[0]; - } - - Interpolator bilinearInt(Interpolator::BiLinearType); - ProcessRubberSheet rubberS; - - //use ProcessRubber sheet to create the sub-cube for this scanz - Trans2d3p transform(trans[i].theta, trans[i].dx, trans[i].dy, - int(sampleTo - sampleFrom), - int(maxL)); - rubberS.SetInputCube(panC[i]); - rubberS.SetOutputCube(tempFile.expanded(), - att, - transform.OutputSamples(), - transform.OutputLines(), - 1); - rubberS.Progress()->SetText( - QObject::tr("Transforming Scan%1: ").arg(i+1).toStdString().c_str()); - rubberS.StartProcess(transform, bilinearInt); - //end process closes the cubes and deletes the dynamically allocated memory (panC[i]) - rubberS.EndProcess(); - - //use process mosaic to add the sub cube to the stiched cube - mosaic.SetInputCube(tempFile.expanded(), - attI, - 1, - 1, - 1, - transform.OutputSamples(), - transform.OutputLines(), - 1); - rubberS.Progress()->SetText(QObject::tr("Mosaicing Scan%1: ").arg(i+1).toStdString().c_str()); - mosaic.StartProcess(int(sampleFrom), 1, 1); - - //clear input cube - mosaic.ClearInputCubes(); - - QFile::remove(tempFile.expanded()); //delete temporary cube if it exists - } - mosaic.EndProcess(); + apolloPanStitcher(ui); } diff --git a/isis/src/apollo/apps/apollopanstitcher/tsts/Makefile b/isis/src/apollo/apps/apollopanstitcher/tsts/Makefile deleted file mode 100644 index 46d84c74c2..0000000000 --- a/isis/src/apollo/apps/apollopanstitcher/tsts/Makefile +++ /dev/null @@ -1,4 +0,0 @@ -BLANKS = "%-6s" -LENGTH = "%-40s" - -include $(ISISROOT)/make/isismake.tststree diff --git a/isis/src/apollo/apps/apollopanstitcher/tsts/default/Makefile b/isis/src/apollo/apps/apollopanstitcher/tsts/default/Makefile deleted file mode 100644 index c583a7b014..0000000000 --- a/isis/src/apollo/apps/apollopanstitcher/tsts/default/Makefile +++ /dev/null @@ -1,24 +0,0 @@ -APPNAME = apollopanstitcher -#This test excercises the stiching of 8 apollopanoramic subscans into a single image. -# The detection of the fiducial marks is complex and involves a RANSAC algorithm. -# For this reason the final stiched cubes are expected to vary slightly. To address -# this the median and number of valid pixels of the delta cube are used to test -# for equivalence of the output. - -include $(ISISROOT)/make/isismake.tsts - -commands: - $(APPNAME) TO=$(OUTPUT)/reduced8.cub \ - FILE_BASE = $(INPUT)/AS15_P_0177R10 \ - MICRONS = 50 > /dev/null; - #difference the two cubes - fx F1=$(OUTPUT)/reduced8.cub F2=$(INPUT)/reduced8.cub TO=$(OUTPUT)/fx.cub EQUATION="F1-F2" > /dev/null; - #use hist to find the median difference of the cubes and total valid pixels - hist FROM=$(OUTPUT)/fx.cub To=$(OUTPUT)/hist.txt NBINS=100000000 > /dev/null; - #form the histogram output grep the median and remove digits beyond the hundreths - cat $(OUTPUT)/hist.txt | grep "Median" \ - | sed 's/\([0-9][0-9]*\.[0-9][0-9]\)\([0-9][0-9]*\)/\1/g' > $(OUTPUT)/CubeStats.txt; - #form the histogram output grep the total valid Pixels - cat $(OUTPUT)/hist.txt | grep "Valid Pixels:" >> $(OUTPUT)/CubeStats.txt; - #remove the intermediate files - $(RM) $(OUTPUT)/fx.cub $(OUTPUT)/hist.txt $(OUTPUT)/reduced8.cub diff --git a/isis/src/base/objs/ProcessMosaic/ProcessMosaic.cpp b/isis/src/base/objs/ProcessMosaic/ProcessMosaic.cpp index 2b9c2bf8f4..5bf4f09afd 100644 --- a/isis/src/base/objs/ProcessMosaic/ProcessMosaic.cpp +++ b/isis/src/base/objs/ProcessMosaic/ProcessMosaic.cpp @@ -707,6 +707,11 @@ namespace Isis { * @throws IException::Message */ Cube *ProcessMosaic::SetOutputCube(const QString &psParameter) { + return SetOutputCube(psParameter, Application::GetUserInterface()); + } + + Cube *ProcessMosaic::SetOutputCube(const QString &psParameter, UserInterface &ui) { + QString fname = ui.GetFileName(psParameter); // Make sure there is only one output cube if (OutputCubes.size() > 0) { @@ -718,7 +723,6 @@ namespace Isis { // (e.g., "TO") and the cube size from an input cube Cube *cube = new Cube; try { - QString fname = Application::GetUserInterface().GetFileName(psParameter); cube->open(fname, "rw"); } catch (IException &) { diff --git a/isis/src/base/objs/ProcessMosaic/ProcessMosaic.h b/isis/src/base/objs/ProcessMosaic/ProcessMosaic.h index ab13a25134..5837656874 100644 --- a/isis/src/base/objs/ProcessMosaic/ProcessMosaic.h +++ b/isis/src/base/objs/ProcessMosaic/ProcessMosaic.h @@ -269,6 +269,7 @@ namespace Isis { // ProcessMosaic objects and child objects, unless redifined in the // child class Isis::Cube *SetOutputCube(const QString &psParameter); + Isis::Cube *SetOutputCube(const QString &psParameter, UserInterface &ui); void SetBandBinMatch(bool enforceBandBinMatch); diff --git a/isis/tests/FunctionalTestsApolloPanStitcher.cpp b/isis/tests/FunctionalTestsApolloPanStitcher.cpp new file mode 100644 index 0000000000..9d590dd0e4 --- /dev/null +++ b/isis/tests/FunctionalTestsApolloPanStitcher.cpp @@ -0,0 +1,42 @@ +#include +#include + +#include "apollopanstitcher.h" +#include "Fixtures.h" +#include "Pvl.h" +#include "PvlGroup.h" +#include "TestUtilities.h" +#include "Histogram.h" +#include "Endian.h" +#include "PixelType.h" + +#include "gtest/gtest.h" +#include "gmock/gmock.h" + +#include "Fixtures.h" + +using namespace Isis; + +static QString APP_XML = FileName("$ISISROOT/bin/xml/apollopanstitcher.xml").expanded(); + +TEST_F(TempTestingFiles, FunctionalTestsApolloPanStitcherDefault) { + QVector args = {"file_base=$ISISTESTDATA/isis/src/apollo/apps/apollopanstitcher/tsts/default/input/AS15_P_0177R10", + "to=" + tempDir.path() + "/reduced8.cub", + "microns=50"}; + + UserInterface options(APP_XML, args); + try { + apolloPanStitcher(options); + } + catch (IException &e) { + FAIL() << "Unable to stitcher apollo images: " << e.toString().toStdString().c_str() << std::endl; + } + Cube outputCube(options.GetFileName("TO")); + + std::unique_ptr hist (outputCube.histogram()); + + EXPECT_DOUBLE_EQ(hist->Average(), 53214.457630315941); + EXPECT_DOUBLE_EQ(hist->Sum(), 3243279908182.748); + EXPECT_EQ(hist->ValidPixels(), 60947345); + EXPECT_DOUBLE_EQ(hist->StandardDeviation(), 20175.877734537076); +}