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



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 loosly 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 indeces 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 thes 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 interation 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 solutins that have average resiudals highter 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


  //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 transfrom from its sample-line
  //sytem to the stiched cube sample-line system Comibine 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;
    }
  }

  //translation from 8 to 8
  trans[7].theta = 0.0;
  trans[7].dx = 0.0;
  trans[7].dy = 0.0;

  //now lets find the extents of the stiched 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
  //transfrom it now into the stiched 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 covert 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 scan
    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();
}
