/**
  *  StarField instance contains the image array and related information
  *  Glen Cowan, RHUL Physics Department
  *  December 2003
  */

import java.util.*;
import java.math.*;

public class StarField {
   
// private data fields
    
  private int nx, ny;
  private int minVal;           // the threshold
  private double gain;
  private double mode;
  private double median;
  private Pixel[][] pixelArray;
  protected Vector clusterVec;
  
// constructors

  public StarField(){
    nx = 0;
    ny = 0;
    minVal = 0;
    gain = 1.0;
    mode = 0.0;
    median = 0.0;
    clusterVec = new Vector();
  }
  
  public StarField(int nx, int ny){
    this();
    this.nx = nx;
    this.ny = ny; 
    pixelArray = new Pixel[nx][ny];
    for (int j=0; j<ny; j++){
      for (int i=0; i<nx; i++){
        pixelArray[i][j] = new Pixel(i, j);
      }
    }
  }
    
// public methods
  
  public int getNx(){ return nx; }
  public int getNy(){ return ny; }
  public void setNx(int nx) { this.nx = nx; }
  public void setNy(int ny) { this.ny = ny; }
  public int getMinVal(){ return minVal; }
  public void setMinVal(int minVal){ this.minVal = minVal; }
  public double getGain(){ return gain; }
  public void setGain(double gain){ this.gain = gain; }
  public double getMode(){ return mode; }
  public double getMedian(){ return median; }
  public void setMode(double mode){ this.mode = mode; }
  public void setMedian(double median){ this.median = median; }
  public Pixel[][] getPixelArray(){ return pixelArray; }
  public void setPixelArray(Pixel[][] pixelArray){
    if ( this.pixelArray == null ){
      int nx = this.getNx();
      int ny = this.getNy();
      this.pixelArray = new Pixel[nx][ny];
    }
    System.arraycopy(pixelArray, 0, this.pixelArray, 0, pixelArray.length);
  }
  public Vector getClusterVec(){ return clusterVec; }
  public void setClusterVec(Vector clusterVec){ this.clusterVec = clusterVec; }
  
  public void fillClusterVec(){
  
    int clusterNum = -1;
    for (int j=0; j<ny; j++){
      for (int i=0; i<nx; i++){
        int val = pixelArray[i][j].getVal();
        int cNum = pixelArray[i][j].getCluster();
        boolean startCluster = val >= minVal & cNum < 0;
        if ( startCluster ){
          clusterNum++;               // start counting at 0 
          Cluster c = new Cluster(clusterNum);
          c.addPixel(pixelArray[i][j]);
          Stack s = new Stack();
          s.push(pixelArray[i][j]);
          boolean clusterFinished = s.empty();
          while ( ! clusterFinished ) {

// Consider neighbouring pixels; if over minVal push on stack and add to cluster

            Pixel p = (Pixel)s.pop();
            int pI = p.getI();
            int pJ = p.getJ();
            if (pI > 0) { considerPixel(pixelArray[pI-1][pJ], s, c, minVal); }
            if (pI < nx-1) { considerPixel(pixelArray[pI+1][pJ], s, c, minVal); }
            if (pJ > 0) { considerPixel(pixelArray[pI][pJ-1], s, c, minVal); }
            if (pJ < ny-1) { considerPixel(pixelArray[pI][pJ+1], s, c, minVal); }

// Should we consider the corner neighbours?  If not, perhaps rethink peak finding algo.
//
//          if (pI>0 && pJ>0) { considerPixel(pixelArray[pI-1][pJ-1], s, c, minVal); }
//          if (pI<nx-1 && pJ>0) { considerPixel(pixelArray[pI+1][pJ-1], s, c, minVal); }
//          if (pI>0 && pJ<ny-1) { considerPixel(pixelArray[pI-1][pJ+1], s, c, minVal); }
//          if (pI<nx-1 && pJ<ny-1) {considerPixel(pixelArray[pI+1][pJ+1], s, c, minVal);}

            clusterFinished = s.empty();

          }             // clusterFinished
          clusterVec.addElement(c);
        }               // start new cluster
      }                 // loop over x
    }                   // loop over y

//    System.out.println("number of clusters found = " + clusterVec.size());
//    System.out.println("cluster num pixels and total counts:");
//    for(int i=0; i<clusterVec.size(); i++){
//      Cluster c = (Cluster)clusterVec.elementAt(i);
//      System.out.println(c.getClusterNum() + "\t" + c.numPixels() + 
//                         "\t" + c.totalCounts());
//    }
 
  }                   // end of fillClusterVec
  
//-------------------------------------------------------------------------

   public void trimClusterVec(int minPixels){
      
//  Eliminate clusters that contain less than the minimum number of pixels.
//  Need to reassign cluster numbers for affected clusters and pixels.

    int numClusters = clusterVec.size();
    int n = 0;
    while ( n < clusterVec.size() ){
      Cluster c = (Cluster)clusterVec.elementAt(n);
      if ( c.numPixels() < minPixels ){
        for (int j=0; j<c.numPixels(); j++){
          Pixel p = (Pixel)c.getPixelVec().elementAt(j);
          p.setCluster(-2);          // tags pixels as former cluster members
        }
        clusterVec.removeElementAt(n);
      }
      else {
        n++;
      }
    }

// reassign the cluster numbers to correspond to the clusVec element number.

    for (int i = 0; i<clusterVec.size(); i++){
      Cluster c = (Cluster)clusterVec.elementAt(i);
      c.setClusterNum(i);
      for (int j=0; j<c.numPixels(); j++){
        Pixel p = (Pixel)c.getPixelVec().elementAt(j);
        p.setCluster(i);
      }
//      System.out.println(i + "  " + c.getClusterNum() + "  " + c.numPixels() );
    }
    System.out.println("number of clusters found = " + clusterVec.size());

  }          // end of trimClusterVec


//-------------------------------------------------------------------------

  public void findPeaks(){

// Within each cluster find maxima (peaks).
// Loop over pixels, for each consider (if not on edge of image) the 8
// nearest neighbours.  If none is greater than the central pixel, a peak is found.
// (Problem:  will find multiple peaks if e.g. several adjacent pixels
// at a local maximum have the same value.)

    for (int k=0; k<clusterVec.size(); k++){
      Cluster c = (Cluster)clusterVec.elementAt(k);
      for (int l=0; l<c.numPixels(); l++){
        Vector pv = c.getPixelVec();
        Pixel pix = (Pixel)pv.elementAt(l);
        int val = pix.getVal();
        int pI = pix.getI();
        int pJ = pix.getJ();
        int a[] = new int[8];
        for(int i=0; i<8; i++){ a[i] = 0; }
        if (pI>0)               { a[0] = pixelArray[pI-1][pJ].getVal(); }
        if (pI>0 && pJ>0)       { a[1] = pixelArray[pI-1][pJ-1].getVal(); }
        if (pJ>0)               { a[2] = pixelArray[pI][pJ-1].getVal(); }
        if (pI<nx-1 && pJ>0)    { a[3] = pixelArray[pI+1][pJ-1].getVal(); }
        if (pI<nx-1)            { a[4] = pixelArray[pI+1][pJ].getVal(); }
        if (pI<nx-1 && pJ<ny-1) { a[5] = pixelArray[pI+1][pJ+1].getVal(); }
        if (pJ<ny-1)            { a[6] = pixelArray[pI][pJ+1].getVal(); }
        if (pI>0 && pJ<ny-1)    { a[7] = pixelArray[pI-1][pJ+1].getVal(); }

        int iMax = -1;
        for (int i=0; i<8; i++){
          if ( a[i] > val ) { iMax = i; }
        }
        if ( iMax < 0 ) { c.addPeak(pix); }

      }             // loop over pixels in a cluster
    }               // loop over clusters
    
  }                 // end of findPeaks

//-------------------------------------------------------------------------

  public void checkPeaks(){

// debug routine for findPeaks

    System.out.println("cluster, peak, i, j...");
    for (int k=0; k<clusterVec.size(); k++){
      Cluster c = (Cluster)clusterVec.elementAt(k);
      Vector pkv = c.getPeakVec();
      int numPeaks = pkv.size();
      for (int l=0; l<numPeaks; l++){
        Pixel pk = (Pixel)pkv.elementAt(l);  // recall a peak is object of type Pixel
        int pI = pk.getI();
        int pJ = pk.getJ();
        System.out.println(k + "  " + l + "  " + pI + "  " + pJ);
      }             // loop over peaks in a cluster
    }               // loop over clusters
    
  }                 // end of checkPeaks

// private methods -----------------------------------------------------------

  private void considerPixel(Pixel n, Stack s, Cluster c, int minVal){

// if the neighbour doesn't already belong to a cluster and is over
// threshold then push onto stack and add to cluster

    if ( n.getCluster() < 0 & n.getVal() >= minVal ) {
      s.push(n);
      int cNum = c.getClusterNum();
      n.setCluster(cNum);
      c.addPixel(n);
    }
  }               // end of considerPixel 
    
}                 // end of class StarField

//-------------------------------------------------------------------------

class Cluster {

// A cluster is a contiguous set of pixels with values of at least
// a minimum threshold.  A peak is a (possibly local) maximum.

  private Vector pixelVec;
  private int clusterNum;
  private Vector peakVec;
  private Star star;

// Constructors

  public Cluster(){
    pixelVec = new Vector();
    clusterNum = -1;
    peakVec = new Vector();
    Star s = new Star(-1, -1, -1);
    setStar(s);
  }

  public Cluster(Vector pixelVec){
    this();
    this.pixelVec = pixelVec;
    if ( pixelVec.size() > 0 ) { 
      clusterNum = ((Pixel)pixelVec.elementAt(0)).getCluster();
    }
    else {
      System.out.println("Cluster>> pixelVec has no entries");
      clusterNum = -1; 
    }
  }

  public Cluster(int clusterNum){
    this();
    pixelVec = new Vector();
    this.clusterNum = clusterNum;
  }

  public Cluster(Vector pixelVec, int clusterNum){
    this();
    this.pixelVec = pixelVec;
    this.clusterNum = clusterNum;
  }

// public methods

  public int numPixels(){ return pixelVec.size(); }
  public Vector getPixelVec(){ return pixelVec; }
  public int getClusterNum(){ return clusterNum; }
  public Vector getPeakVec(){ return peakVec; }

  public void addPixel(Pixel p){
    pixelVec.addElement(p);
    p.setCluster(clusterNum);       // recall clusterNum is private data of Pixel
  }
  
  public void addPeak(Pixel pix){
    peakVec.addElement(pix);
    if ( pix.getCluster() != clusterNum ) {
      System.out.println("Cluster.addPeak>> peak cluster number doesn't match");
      pix.setCluster(clusterNum);
    }
  }

  public void setClusterNum(int clusterNum){
    this.clusterNum = clusterNum;
  }
  
  public int totalCounts(){
    int counts = 0;
    for (int i = 0; i<pixelVec.size(); i++){
      counts += ((Pixel)pixelVec.elementAt(i)).getVal();
    }  
    return counts;
  }

  public void setStar (Star star){
    this.star = star;
  }

  public Star getStar(StarField sf){
    if (star.getX() < 0 || star.getY() < 0 || star.getVal() < 0) {
      return calculateStar(sf);
    }
    else {
      return star;
    }
  }

// for barycentre calculation, subtract threshold minVal.  For this need
// to pass StarField.

  public Star calculateStar(StarField sf){
    double xValSum = 0, yValSum = 0, valSum = 0;
    for (int i=0; i<pixelVec.size(); i++){
      int pixVal = (int)(((Pixel)pixelVec.elementAt(i)).getVal());
      double val = (double)(pixVal - sf.getMinVal());
      double x = (double)(((Pixel)pixelVec.elementAt(i)).getI());
      double y = (double)(((Pixel)pixelVec.elementAt(i)).getJ());
      xValSum += x*val;
      yValSum += y*val;
      valSum += val;
    }
    double muX = xValSum / valSum;
    double muY = yValSum / valSum;
    Star s = new Star(muX, muY, valSum);
    return s;
  }
  
}             // end of class Cluster

//-------------------------------------------------------------------------

class Pixel {

  private int i, j;             // the x and y indices
  private int val;              // in ADU
  private int cluster;          // negative if not member of cluster

  public Pixel(int i, int j){
    this.i = i;
    this.j = j;
    val = 0;
    cluster = -1;
  }

  public Pixel(int i, int j, int val){
    this(i, j);
    this.val = val;
  }

  public Pixel(int i, int j, int val, int cluster){
    this(i, j, val);
    this.cluster = cluster;
  }

  public int getI(){ return i; }
  public int getJ(){ return j; }
  public int getVal(){ return val; }
  public int getCluster(){ return cluster; }
  public void setCluster(int cluster){
    this.cluster = cluster;
  }

}         // end of class Pixel

//-------------------------------------------------------------------------

class Star implements Comparable {

//  A star's position and intensity are stored as doubles (in contrast to Pixel).

  private double val;
  private double x, y;
  protected double nStar;        // background subtracted value
  protected double sigmaNStar;   // Poisson standard deviation

  public Star(){
    this.val = 0;
    this.x = 0;
    this.y = 0;
    this.nStar = 0;
    this.sigmaNStar = 0;
  }

  public Star(double val){
    this();
    this.val = val;
  }

  public Star(double x, double y){
    this();
    this.x = x;
    this.y = y;
  }

  public Star(double x, double y, double val){
    this(x, y);
    this.val = val;
  }

  public double getVal(){ return val; }
  public void setVal(double val){ this.val = val; }
  public double getX(){ return x; }
  public void setX(double x){ this.x = x; }
  public double getY(){ return y; }
  public void setY(double y){ this.y = y; }
  public double getNStar(){ return nStar; }
  public void setNStar(double nStar){ this.nStar = nStar; }
  public double getSigmaNStar(){ return sigmaNStar; }
  public void setSigmaNStar(double sigmaNStar){ this.sigmaNStar = sigmaNStar; }

  public int compareTo (Object o){

// define comparison so that sorted order is descending

    Star s = (Star)o;
    double val1 = this.getVal();
    double val2 = s.getVal();
    int comp = 0;
    if ( val1 < val2 ) {
      comp = 1;
    } 
    else if ( val2 > val1 ){
      comp = -1;
    }
    else {
      comp = 0;
    }
    return comp;

  }          // end of compareTo

}         // end of class Star