CentroidApolloPan.cpp

//Isis.h and IsisDebug.h if needed


//C++ standard libraries if needed
#include <math.h>
#include <time.h>

//QT libraries if needed if needed


//third party libraries if needed
#include <gsl/gsl_rng.h>
#include "GSLUtility.h"

//Isis Headers if needed
#include "CentroidApolloPan.h"


using namespace std;

namespace Isis {

  CentroidApolloPan::CentroidApolloPan(double pixelSizeMicrons)
  {
    if( pixelSizeMicrons > 0)
      this->m_pixelSize = pixelSizeMicrons;
    else
      this->m_pixelSize = 5.0;  //if a negative or zero value is passed set the default
  }

  
  CentroidApolloPan::~CentroidApolloPan(){};

  bool CentroidApolloPan::setPixelSize(double microns)
  {
    if( microns <= 0)return false;
      this->m_pixelSize = microns;
    return true;
  }

  int CentroidApolloPan::selectAdaptive(Chip *inputChip,Chip *selectionChip) {
    /*
    Given a range of DN this function creates a biniary chip for all continuous pixels that have the 
    DN within the specified range using the center pixel of the chip as the seed value
    input:
      m_minDN,m_maxDN  set using the this->SetDNRange(..)
      inputChip  chip centered around some centroidable feature (dark or light continous block of 
      pixels)
    
    output:
      selectionChip  binary chip of selected and unselected pixels

    */
    
    //check the sizes of the chips and make the selectionChip match the inputChip
    int lines, samples;

    int i,j,k,l,step;

    double minDN,maxDN,temp;

    vector< double> dn;
    
    lines = inputChip->Lines();
    samples = inputChip->Samples();
    //printf("DEBUG: samples: %d  lines: %d\n",samples,lines);

    if (lines <= 0 || samples <= 0)
      return 0;  //abort if the input chips isn't 2D

    //check DN's around the edge, we want this centroid to stop before the edge of the chip
    minDN = this->getMinDN();
    maxDN = this->getMaxDN();

    //build an array of all the DN around the border of the chip
    for (i=1;i<=lines;i++) {
      if (i == 1 || i == lines) {
        step = 1;
      }
      else {
        step = samples-1;
      }
      for(j=1;j<=samples;j+=step) {
        dn.push_back(inputChip->GetValue(j,i));
      }
    }

    //find the 90 percentile of dn
    j = dn.size();
    for (i=0;;i++) {
      //find the largest ellement of dn
      l=0;
      temp = dn[0];
      for (k=1;k<int(dn.size());k++) {
        if (dn[k] > temp) {
          temp = dn[k];
          l =k;
        }
      }
      //erase the largest ellement of dn or quit
      if (i < int(j/10)) dn.erase(dn.begin()+l);
      else break;
    }

    if( temp > minDN) this->setDNRange(temp,maxDN);  //set updated DN range

    //printf("DEBUG: minDN: %lf\n",temp); getchar();

    this->select(inputChip,selectionChip);

    this->setDNRange(minDN,maxDN); //restore original DN range   

    return 1;
  }  


  int CentroidApolloPan::elipticalReduction(Chip *selectionChip, double percent_selected, 
                                            double play, int patience_limit)
  {
    int i,j,k,l,
        emptySets,  //number of consecutive failures to find a better ellipse
        randPt,    //random potential center point selection
        samples,  //samples in selectionChip
        lines;    //lines in selectionChip

    std::vector< std::vector<int> > pts;      //vector of pixels on the border of the slection
    std::vector< std::vector<double> >  cen_pts;    //vector of possible center pixels

    std::vector<double> pt(2);  //individual pixel coordinate

    Ellipse ellNew,ellBest;  //new and best ellipses

    double  a = 60.0*5.0/m_pixelSize,  //approximate semiMajor axis length
            b = 60.0*5.0/m_pixelSize,  //approximate semiMinor axis length
            dpt[2];        //double 2D point
      
    //printf("DEBUG: a: %lf b: %lf pixel_size: %lf\n",a,b,m_pixelSize);    
    gsl_rng *randGen = gsl_rng_alloc(gsl_rng_taus);
    gsl_rng_set(randGen, time(NULL));  //using the current time as a seed value
    
    samples = selectionChip->Samples();
    lines   = selectionChip->Lines();
    //printf("DEBUG selectionchip samples: %d lines: %d\n",samples,lines);

    //populating a vector of all possible integral center points
    k = int(ceil(double(samples) - a));
    l = int(ceil(double(lines  ) - b));
    j=int(floor(b+1.0));
    i=int(floor(a+1.0));

    cen_pts.clear();
    
    /*
     * inorder to concentrate the search near the center of the chip we'll varry the step size for 
     * the center search nodes quadratically from 1*5.0/m_pixelSize at the center to 
     * 5.0*5.0/m_pixelSize at the edge
     */
    double searchStep = 5.0*5.0/m_pixelSize;
    double aSearch,centerSample,centerLine,delta,cSearch;
    cSearch = searchStep/5.0;
    centerSample = double((k+i)/2.0);
    centerLine   = double((l+j)/2.0);
    delta = (centerSample-double(i))*(centerSample-double(i)) + (centerLine  -double(l))*(centerLine  -double(l));  //square of the maximum distance from center of the search space
   
    aSearch = (searchStep-cSearch)/delta;

    for (pt[0]=floor(a+1.0); pt[0]<k; pt[0]+=searchStep) {
      for (pt[1]=floor(b+1.0); pt[1]<l; pt[1]+=searchStep) {
        cen_pts.push_back(pt);
        delta = (centerSample-double(pt[0]))*(centerSample-double(pt[0])) + (centerLine-double(pt[1]))*(centerLine-double(pt[1]));  //square of the distance from center of the search space
        searchStep = aSearch*delta + cSearch;
        if (searchStep > l - pt[1] && l - pt[1] > 1e-4)  //don't jump outside the search area
          searchStep = l - pt[1];
      }
      if (searchStep > k - pt[0] && k - pt[0] > 1e-4)  //don't jump outside the search area
        searchStep = k - pt[0];
    }

    //STEP 1: finding points along the boundary of the selection
    selectionEdge(selectionChip,&pts);
    //printf("DEBUG: edgePoints: %d\n",pts.size());

    if (pts.size() ==0) return 0;
  
    emptySets = 0;
    ellBest.area = 0.0;  //initialize the current best ellipse to zero area

    while (emptySets < patience_limit && cen_pts.size()>0) {
      //printf("DEBUGA\n");
      //STEP 2: choosing a new random hypothesis point
      randPt = gsl_rng_uniform_int( randGen  ,cen_pts.size() );
      dpt[0] = cen_pts[randPt][0];
      dpt[1] = cen_pts[randPt][1];

      cen_pts.erase( cen_pts.begin() + randPt);  //erasing this ellement ensures it will never be chosen again
      //printf("DEBUGB\n");
      //STEP 3: Now define the ellipse
      ellipseFromCenterAxesAngle(&ellNew, dpt[0], dpt[1], a, b, 0.0);
      //printf("Debug a: %lf b: %lf theta: %lf center: %lf %lf\n",ellNew.semiMajor,ellNew.semiMinor,ellNew.majorAxis[0],ellNew.cen[0],ellNew.cen[1]);getchar();
      //printf("DEBUGC\n");
      //STEP 4: Do a least squares generization.  Any points within a distance of play pixels to the edge of the ellipse are included in the gernalization. The distance check is repeated for every iteration so the ellipse can effectively grow to include more points.
      if (!bestFitEllipse(&ellNew, &pts, play, 50)) {
        emptySets++;
        continue;  //if optimization fails go on to the next hypothesised ellipse
      }
      //printf("DEBUGD\n");
      //STEP 5: If the generization is successfully check to see if the area is at least as great as the current best, and that the ellipse is contained in the chip
      if (ellNew.area < ellBest.area) {
        emptySets++;
        continue;
      }
      //printf("DEBUGE\n");
      if (!ellipseInChip(&ellNew,selectionChip)) {  //if the ellipse is entirely contained in the chip
        emptySets++;
        continue;
      }
      //printf("DEBUGF\n");
      //STEP 6: if the area is great enough check that the percent sellected is at least percent_selected.
      if (elipsePercentSelected(selectionChip,&ellNew) < percent_selected) {
        emptySets++;
        continue;
      }

      //STEP 7 saving the new Best ellipse so far
      //printf("DEBUG: new best area: %lf\n",ellNew.area);
      emptySets=0;
      ellBest.area = ellNew.area;
      ellBest.A[0] = ellNew.A[0];
      ellBest.A[1] = ellNew.A[1];
      ellBest.A[2] = ellNew.A[2];

      ellBest.cen[0] = ellNew.cen[0];
      ellBest.cen[1] = ellNew.cen[1];

      ellBest.semiMajor = ellNew.semiMajor;
      ellBest.semiMinor = ellNew.semiMinor;
      ellBest.majorAxis[0] = ellNew.majorAxis[0];
      ellBest.majorAxis[1] = ellNew.majorAxis[1];
      ellBest.minorAxis[0] = ellNew.minorAxis[0];
      ellBest.minorAxis[1] = ellNew.minorAxis[1];
    }

    if( ellBest.area == 0)
      return 0;  //no ellipse meeting the selection criteria was found


    //go through and unselect the points outside the trimming ellipse
    for (i=1;i<=samples;i++) {
      for (j=1;j<=lines;j++) {
        dpt[0] = double(i) ;
        dpt[1] = double(j) ;
              
        if (!pointInEllipse(&ellBest,dpt,play)) {  //if the point isn't within play pixles of being within the elipse
          if( selectionChip->GetValue(i,j) == 1)
            selectionChip->SetValue(i,j,3.0);
        }
      }
    }
    return 1;

  }
  
}