/*
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 <math.h>
#include <stdio.h>

//QT libraries if needed if needed
#include <QFile>

//third party libraries if needed
#include <Ransac.h>

//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.
 *
 */



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 <TRANS2d3p> 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 <std::vector <double> > scanFid;  //fiducials of Scan i
  std::vector <std::vector <double> > conFid;  //fiducials of Scan i
  std::vector <double> fid(2);      //individual fiducial measurement
  std::vector <double> 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("$apollo15/templates/apolloPanFiducialFinder.pvl");
  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<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]);
              temp = sqrt(wdot[0]*wdot[0] + wdot[1]*wdot[1]);

              if (temp > 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();
}
