target_points.java

package com.elphel.imagej.calibration;
/**
** -----------------------------------------------------------------------------**
** target_points.java
**
** Measures focus sharpnes at orthogonal directions, differenct colors,
** Displays results for manual focusing/image plane tilting.
** NOTE: Requires special targets !
**
** Copyright (C) 2010 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**  
**  target_points.java is free software: you can redistribute it and/or modify
**  it under the terms of the GNU General Public License as published by
**  the Free Software Foundation, either version 3 of the License, or
**  (at your option) any later version.
**
**  This program is distributed in the hope that it will be useful,
**  but WITHOUT ANY WARRANTY; without even the implied warranty of
**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**  GNU General Public License for more details.
**
**  You should have received a copy of the GNU General Public License
**  along with this program.  If not, see <http://www.gnu.org/licenses/>.
** -----------------------------------------------------------------------------**
**
*/

import ij.*;
import ij.process.*;
import ij.gui.*;

import java.awt.*;
import java.awt.event.*;

import ij.plugin.frame.*;

import java.util.List;

import com.elphel.imagej.jp4.JP46_Reader_camera;

import java.util.ArrayList;

import ij.text.*;

import java.lang.Integer;


public class target_points extends PlugInFrame implements ActionListener {
private static final long serialVersionUID = -3057496866952930812L;
JP46_Reader_camera jp4_instance;
//  MTF_Bayer MTF_Bayer_instance;
  Panel panel;
  static Frame instance;

 public static int DEBUG_LEVEL = 1;
 public static int MASTER_DEBUG_LEVEL = 1;

 public static int FFTSize=64;
 public static int FFTScanStep=8;
 
 public static int test_x=FFTSize;
 public static int test_y=FFTSize;
 public static int     displayWidth=800;
 public static int     displayHeight=600;
 public static float [][] input_bayer=null;
 public static float [][] convolved_bayer=null;
 public static float [][] normalized_convolved_bayer=null;
 
 public static float   [] target_kernel=null;
 public static int   [][] clusterMaps=null; 
 public static double [][][] targetCoordinates;

 public static double targetOuterDMin =38; // minimal outer diameter of the target image , in pixels
 public static double targetOuterDMax =47; // maximal outer diameter of the target image , in pixels
 public static int    numTargetRings  = 2;  // number of target black rings (notg counting center black circle)
 public static int    pixelsSubdivide =10;  // subdivide pixels by this number (each direction) when generating targets
 public static double deconvInvert =   0.1; // when FFT component is lass than this fraction of the maximal value, replace 1/z with Z
 public static double filteredRadius=   7.0;  //pix - search maximums after convolution in (2*filteredRadius+1) squares
 public static double backgroundRadius= 15.0; //pix - consider ring area between backgroundRadius and filteredRadius as reference
 public static double clusterThreshold= 1.0; //1.2; // tested with 0.8 - many extras, but filtered out
 public static int    clusterSize=      20;  // cluster size (will be expanded/shrank before finding centroid

 public static int    discrAngularFreq=     2 ;  // pixels on FFT image of tragets converted polar (the smaller, the less angular variations)

/// these parameters are dependent on targets, use debug mode and manula fft for 64x64 polar coordinates target areas
 public static int    discrRadialMinFreq=  7 ; // pixels on FFT image of targets converted polar (radial component)
 public static int    discrRadialMaxFreq=  9 ; // pixels on FFT image of targets converted polar (radial component)
 public static double discrThreshold=      0.1; // FFT energy fraction in selecter area should be > this threshold to pass the test
 public static double maxChromaticDistance= 10.0; // Maximal distance between the same target on different color copmponents


 public static double[][][] targets; // For each target {averageX, averageY, num_non_zero_components},{X1,Y1,Qulaity1},...,{X4,Y4,Qulaity4}
/**
discrRadialMinFreq=(size*(2* numTargetRings +1)/targetOuterDMax)-1
discrRadialMaxFreq=(size*(2* numTargetRings +1)/targetOuterDMin)+1
*/

 private ImagePlus     imp_src;
 public ImageProcessor ip_display;
 public ImagePlus      imp_display;
 public ImagePlus      imp_camera=null;
 
  Plot plotResult;
 

  public target_points() {
    super("target_points");
    if (IJ.versionLessThan("1.39t")) return;
    if (instance!=null) {
      instance.toFront();
      return;
    }
    instance = this;
    addKeyListener(IJ.getInstance());

    setLayout(new FlowLayout());
    panel = new Panel();
    addButton("Configure");
    addButton("Split Bayer");
    addButton("Create Target");
    addButton("Split&Convolve");
    add(panel);
    pack();
    GUI.center(this);
    setVisible(true);
    initHamming(FFTSize);
//    initDisplay();
    jp4_instance=       new JP46_Reader_camera();
  }
  void addButton(String label) {
    Button b = new Button(label);
    b.addActionListener(this);
    b.addKeyListener(IJ.getInstance());
    panel.add(b);
  }
  public void actionPerformed(ActionEvent e) {
    int i,j,ir2, size;
    String label = e.getActionCommand();
    if (label==null) return;
    if (label.equals("Configure")) {
      if (showDialog()) {
        initHamming(FFTSize);
      }
      return;
    }
    if (label.equals("Create Target")) {
      DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
      target_kernel=createTargetDialog ();
      return;
    }
    imp_src = WindowManager.getCurrentImage();
    String newTitle= imp_src.getTitle();
    Rectangle r=new Rectangle(imp_src.getWidth(),imp_src.getHeight());

    if (label.equals("Split Bayer")) {
      DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
      input_bayer=splitBayer (imp_src);
      showBayers(input_bayer, r.width>>1, r.height>>1, newTitle);
      return;
    }
    if (label.equals("Split&Convolve")) {
      DEBUG_LEVEL=MASTER_DEBUG_LEVEL;
      input_bayer=splitBayer (imp_src);
      if (DEBUG_LEVEL>5) showBayers(input_bayer, r.width>>1, r.height>>1, newTitle);
       
      if (target_kernel==null) {
        IJ.showMessage("Error","Target kernel does not exist, please generate one");
        return;
      }
      size=(int) Math.sqrt(target_kernel.length);
//target_kernel
/* Convolve Bayer slices with prepared target template  */
      convolved_bayer=new float[input_bayer.length][];
      for (i=0; i<input_bayer.length; i++) {
        IJ.showStatus("Convolving Bayer "+i);
/* Double in convolution works twice faster than float!*/
        convolved_bayer[i]=doubleConvolveWithTarget(input_bayer[i], target_kernel, r.width>>1, r.height>>1, size);

      }
      if (DEBUG_LEVEL>2) showBayers(convolved_bayer, r.width>>1, r.height>>1, newTitle+"_"+deconvInvert);
//
//    filteredRadius   =(int) gd.getNextNumber();
//    backgroundRadius =(int) gd.getNextNumber();
/* normalize convolved Bayer slices   */
/* prepare pixel mask for the normalization (ring) */
      int filtSize= (int) filteredRadius;
      boolean [][] mask=new boolean[2*filtSize+1][2*filtSize+1];
      for (i=0;i<(2*filtSize+1);i++) for (j=0;j<(2*filtSize+1);j++) {
        ir2=(i-filtSize)*(i-filtSize) + (j-filtSize)*(j-filtSize);
        mask[i][j]=(ir2>(filteredRadius*filteredRadius)) && (ir2 < (backgroundRadius*backgroundRadius));
      }
// public static float [][] normalized_convolved_bayer=null;
      normalized_convolved_bayer=new float[input_bayer.length][];
      for (i=0; i<input_bayer.length; i++) {
        IJ.showStatus("Normalizing Bayer "+i);
        normalized_convolved_bayer[i]=normalizeAtRing(convolved_bayer[i], r.width>>1, r.height>>1, mask);
      }
      if (DEBUG_LEVEL>1) showBayers(normalized_convolved_bayer, r.width>>1, r.height>>1, newTitle+"_"+deconvInvert+"_normalized");
// public static int   [][] clusterMaps=null; 
//clusterThreshold
      clusterMaps= new int[4][];
      targetCoordinates=new double[4][][];
      for (i=0; i<input_bayer.length; i++) {
        IJ.showStatus("Clusterizing Bayer "+i);
//        clusterMaps[i]=clusteriseTargets(normalized_convolved_bayer[i], r.width>>1, r.height>>1,clusterThreshold,clusterSize);
        targetCoordinates[i]=clusteriseTargets(normalized_convolved_bayer[i], r.width>>1, r.height>>1,clusterThreshold,clusterSize);
      }
      float[][]clusterPixels=new float[4][input_bayer[0].length];
      if (DEBUG_LEVEL>2) {
//         float[][]clusterPixels=new float[4][input_bayer[0].length];
         for (i=0; i<clusterPixels.length; i++) {
          for (j=0; j<clusterPixels[0].length; j++) clusterPixels[i][j]=0.0f;
          for (j=0;j<targetCoordinates[i].length;j++) {
             clusterPixels[i][((int)(Math.round(targetCoordinates[i][j][1])*(r.width>>1))) +((int)Math.round(targetCoordinates[i][j][0]))]=1.0f;
          }
        }
        showBayers(clusterPixels, r.width>>1, r.height>>1, newTitle+"_"+deconvInvert+"_clusters");
      }
      float [] rectTarget;
      ImageProcessor ip_dbg;
      ImagePlus      imp_dbg;
      double          likely;
      double [] likelyness;
      int numGoodTargets;
      double [][] goodTargets; /// x/y/above Threshold
      double r0=(targetOuterDMin+targetOuterDMax)/(2* numTargetRings +1)/8; /// copmpensate for the center circle twice wider than rings
/**
discrRadialMinFreq=(size*(2* numTargetRings +1)/targetOuterDMax)-1
discrRadialMaxFreq=(size*(2* numTargetRings +1)/targetOuterDMin)+1
*/
      for (i=0; i<clusterPixels.length; i++) {
        likelyness= new double[targetCoordinates[i].length];
        numGoodTargets=0;
        for (j=0;j<targetCoordinates[i].length;j++) {
          rectTarget=circle2DoubleRect (input_bayer[i], r.width>>1, r.height>>1, size,targetCoordinates[i][j][0],targetCoordinates[i][j][1], r0);
          likely=likelyTarget(rectTarget, size,discrAngularFreq,discrRadialMinFreq,discrRadialMaxFreq);
          likelyness[j]=likely;
          if (DEBUG_LEVEL>2) {
   System.out.println("Cluster="+j+" x="+targetCoordinates[i][j][0]+" y="+targetCoordinates[i][j][1]+" likely="+likely);

          }
          if (DEBUG_LEVEL>3) {
            if (i==0) {  // just to reduce debug clutter
              ip_dbg=new FloatProcessor(size,size);
              ip_dbg.setPixels(rectTarget);
              ip_dbg.resetMinAndMax();
              imp_dbg=  new ImagePlus(newTitle+"_rect_"+j, ip_dbg);
              imp_dbg.show();
            }
          }
          if (likely >=discrThreshold) numGoodTargets++;
        }
        goodTargets=new double[numGoodTargets][3];
        numGoodTargets=0;
        for (j=0;j<targetCoordinates[i].length;j++) {
          if (likelyness[j]>=discrThreshold) {
            goodTargets[numGoodTargets][2]=likelyness[j]/discrThreshold;
            goodTargets[numGoodTargets][0]=2*targetCoordinates[i][j][0]+(i&1);      // convert to full image, use Bayer shift
            goodTargets[numGoodTargets][1]=2*targetCoordinates[i][j][1]+((i>>1)&1); // convert to full image, use Bayer shift
            numGoodTargets++;
          } 
        }
        targetCoordinates[i]=goodTargets;
      }

      if (DEBUG_LEVEL>1) {
         for (i=0; i<clusterPixels.length; i++) {
          for (j=0; j<clusterPixels[0].length; j++) clusterPixels[i][j]=0.0f;
          for (j=0;j<targetCoordinates[i].length;j++) {
             clusterPixels[i][((int)(Math.round((targetCoordinates[i][j][1] - ((i>>1)&1))/2) *(r.width>>1))) +
                              ((int)(Math.round((targetCoordinates[i][j][0] - ( i    &1))/2)))]=1.0f;
   System.out.println("Bayer="+i+" Target="+(j+1)+" x="+targetCoordinates[i][j][0]+" y="+targetCoordinates[i][j][1]+" quality="+targetCoordinates[i][j][2]);
          }
        }
        showBayers(clusterPixels, r.width>>1, r.height>>1, newTitle+"_"+deconvInvert+"_clusters");
      }
/* TODO: Verify that all Bayer components have the same targets (build composite table) */
      targets= combineTargets(targetCoordinates, maxChromaticDistance);
      showTargetsLocationTable(targets,  newTitle, 2,  (DEBUG_LEVEL>1));
      return;
    }
  }

/* Combine target locations from 4 Bayer components */
  double [][][] combineTargets(double[][][] targetCoord, ///[bayer][number][x,y,q>1]
                               double maxDistance) {     /// maximal distance between the same target in different Bayer components

   int [][] targetNumbers=new int[targetCoord.length][];
   int i,i1,j,k,n,l,maxLen, bNum, numTargets;
   double d2=maxDistance*maxDistance;
   double [][][] targets;
   double avX,avY;
   maxLen=0;
   bNum=-1;
   for (i=0;i<targetCoord.length;i++) {
      l=targetCoord[i].length;
      targetNumbers[i]=new int [l];
      for (j=0;j<l;j++) targetNumbers[i][j]=0;
      if (maxLen<l){
        maxLen=l;
        bNum=i;
      }
   }
 /* assign target number according to the component that has most of the targets (does not mean others do not have that this one is missing */
   for (j=0;j<maxLen;j++) targetNumbers[bNum][j]=j+1;
   numTargets=maxLen; // may increase later
 /* compare all other color components with the coordinates in the seslected one (not too many to bother with good guess) */
   for (i=0;i<targetNumbers.length;i++) if (i!=bNum) {
     for (j=0;j<targetNumbers[i].length;j++) if (targetNumbers[i][j]==0){
        for (k=0;k<targetNumbers[bNum].length;k++) {
          if (((targetCoord[bNum][k][0]-targetCoord[i][j][0])*(targetCoord[bNum][k][0]-targetCoord[i][j][0])+
               (targetCoord[bNum][k][1]-targetCoord[i][j][1])*(targetCoord[bNum][k][1]-targetCoord[i][j][1]))<=d2) {
              targetNumbers[i][j]=targetNumbers[bNum][k];
            break;
          }
        }
     }
   }
/* See if any targets are missing, add them */
   for (i=0;i<targetNumbers.length;i++) if (i!=bNum) {
     for (j=0;j<targetNumbers[i].length;j++) if (targetNumbers[i][j]==0){
       numTargets++;
       targetNumbers[i][j]=numTargets;
       for (i1=i+1;i1<targetNumbers.length;i1++) if (i1!=bNum) {
         for (k=0;k<targetNumbers[i1].length;k++) if (targetNumbers[i1][k]==0) {
           targetNumbers[i1][k]=numTargets;
         }
       }
     }
   }
   if (DEBUG_LEVEL>2) {
     System.out.println("numTargets="+numTargets);
     for (i=0;i<targetCoord.length;i++) for (j=0;j<targetCoord[i].length;j++) {
       System.out.println("["+targetNumbers[i][j]+"] "+i+":"+j+" "+targetCoord[i][j][0]+","+targetCoord[i][j][1]+" :"+targetCoord[i][j][2]);
     }
   }




   targets = new double [numTargets][targetNumbers.length+1][3];
   for (i=0;i<numTargets;i++) for (j=0;j<targets[i].length;j++) for (k=0;k<3;k++) targets [i][j][k]=0.0;
   for (i=0;i<targetNumbers.length;i++) for (j=0;j<targetNumbers[i].length;j++) if (targetNumbers[i][j]!=0){ // should be always non-zero
     k=targetNumbers[i][j]-1;
     targets[k][i+1][0]=targetCoord[i][j][0]; // x
     targets[k][i+1][1]=targetCoord[i][j][1]; // y
     targets[k][i+1][2]=targetCoord[i][j][2]; // quality >1.0
   }
/* Calculate average values*/
   for (i=0;i<numTargets;i++) {
     avX=0.0;
     avY=0.0;
     n=0;
     for (j=1;j< targets[i].length; j++) if (targets[i][j][2]>0){
       avX+=targets[i][j][0];
       avY+=targets[i][j][1];
       n++;
     }
     if (n>0) {
       targets[i][0][0]=avX/n;
       targets[i][0][1]=avY/n;
       targets[i][0][2]=n;
     }
   }
   if (DEBUG_LEVEL>2) {
     System.out.println("targets");
     for (i=0;i<targets.length;i++) {
       System.out.println(i+" | "+targets[i][0][0]+","+targets[i][0][1]+" :"+targets[i][0][2]+
                            " | "+targets[i][1][0]+","+targets[i][1][1]+" :"+targets[i][1][2]+
                            " | "+targets[i][2][0]+","+targets[i][2][1]+" :"+targets[i][2][2]+
                            " | "+targets[i][3][0]+","+targets[i][3][1]+" :"+targets[i][3][2]+
                            " | "+targets[i][4][0]+","+targets[i][4][1]+" :"+targets[i][4][2]);
     }
   }

   return targets;
  }

  public void showTargetsLocationTable(double[][][] targets, String title, int precision, boolean showQuality) {
    int i,n;
    
    String header="#\tX\tY";
    for (i=1;i<targets[0].length;i++) header+="\tdX"+i+"\tdY"+i+(showQuality?("\tQ"+i):"");
    StringBuffer sb = new StringBuffer();
    for (n=0;n<targets.length;n++) {
      sb.append((n+1)+
                  "\t"+IJ.d2s(targets[n][0][0],precision)+  // Average X 
                  "\t"+IJ.d2s(targets[n][0][1],precision)); // Average Y 
      for (i=1;i<targets[0].length;i++) {
        if (targets[n][i][2]>0) {
          sb.append(  "\t"+(((targets[n][i][0]-targets[n][0][0])>0)?"+":"")+IJ.d2s(targets[n][i][0]-targets[n][0][0],precision)+  // X
                      "\t"+(((targets[n][i][1]-targets[n][0][1])>0)?"+":"")+IJ.d2s(targets[n][i][1]-targets[n][0][1],precision)); // Y
        } else {
          sb.append(  "\t---\t---");
        }
        if (showQuality) sb.append("\t"+((targets[n][0][2]>0)?IJ.d2s(targets[n][i][2],precision):"---"));
      }
      sb.append(  "\n");
    }
    new TextWindow(title+"_Target_Locations_Table", header, sb.toString(), showQuality?900:700,500);
  }



/* determines if it was likely a target of concentric circles, after convertion to polar coordinates expect nearly horizontal b/w bands */
double likelyTarget(float[] pixels, // pixel array
                          int size, // image size (should be square for FFT
                          int hor,  // horizontal selection area (half width)
                          int vertMin,
                          int vertMax) // vertical selection area (half height > half width for horizointal bands)
 {

      ImageProcessor ip;
      FHT            fht;
      double[][][]   fft;
      double         s1,s2,e;
      int i,j;
      ip=new FloatProcessor(size,size);
      ip.setPixels(pixels);
      fht =  new FHT(ip);
// Swapping quadrants, so the center will be 0,0
      fht.swapQuadrants();
// get to frequency
      fht.transform();
// Convert from FHT to complex FFT
      fft= FHT2FFTHalf (fht);
      s1=0;
      s2=0;
     
      for (i=0;i<(size/2+1);i++) for (j=0;j<size;j++) {
       if ((i>0) || (j>0)) { // skip DC
         e=fft[i][j][0]*fft[i][j][0]+fft[i][j][1]*fft[i][j][1];
         if ((i>=vertMin) && (i<=vertMax) && ((j<=hor) || (j>=size-hor))) s1+=e;
         else s2+=e;
       }
      }
      return s1/(s1+s2); // fraction inside selected area, use as likelyhood of the needed target
 }

float [] circle2DoubleRect (float [] pixels, int width, int height, int size, double x0, double y0, double r0) {
  float [] outPixels=new float[size*size];
  int x,y;
  double a,r;
  int px,py;
  for (y=0;y<size;y++) for (x=0;x<size;x++) {
    if ((y>(size>>1)) ||  ((y==(size>>1)) && (x>=(size>>1)))){
      r=y-(size>>1); // +0.5?
      a=size-x-1;
    } else {
      r=(size>>1)-y; // -0.5?
      a=size+x;
    }
    r+=r0; /// to match periodic pattern on both sides of zero (center circle is twice wider)
    a*=Math.PI/size;
    px=((int)Math.round(x0+r*Math.cos(a)) + width ) % width;
    py=((int)Math.round(y0+r*Math.sin(a)) + height) % height;
    outPixels[y*size+x]=pixels[py*width+px];    
  }
  return outPixels;
}

  double [][] clusteriseTargets(float [] pixels,int width, int height, double threshold, int clusterSize) {
    if ((width*height) != pixels.length) {
      IJ.showMessage("Error in  clasteriseTargets","pixels.length ("+pixels.length+") does not match width ("+width+") x height ("+height+") = "+(width*height));
      return null;
    }
    
    int x,y,i,j;
    Integer Index, NewIndex, NextIndex;
    int clusterNumber=1;
    int []clusterMap=new int[width*height];
    List <Integer> pixelList=new ArrayList<Integer>(100);
    int [] dirs={1,-width+1,-width,-width-1,-1,+width-1,width,width+1};
    int listIndex;
    float f;
    boolean first;
    double cx,cy,cm; // for centroid calculation;
    int ix,iy;
    Double [] cxy=null;
    for (i=0;i<clusterMap.length;i++) clusterMap[i]=0; /// 0 - unused, -1 - "do not use"
    List <Double[]> Centroids=new ArrayList<Double[]>(100);
    for (y=0;y<height;y++) for (x=0;x<width;x++) {
      if ((pixels[y*width+x]>=threshold) && (clusterMap[y*width+x]==0)) {
/// mark all connected pixels above the threshold
         Index=y*width+x;
         pixelList.clear();
         pixelList.add (Index);
         clusterMap[Index]=clusterNumber;
         listIndex=0;
         while (listIndex<pixelList.size() ) {
           Index=pixelList.get(listIndex++);
           for (i=0;i<dirs.length;i++) {
             NewIndex=Index+dirs[i];
             if ((NewIndex>=0) && (NewIndex<clusterMap.length) && (clusterMap[NewIndex]==0) && (pixels[NewIndex]>=threshold)) {
               pixelList.add (NewIndex);
               clusterMap[NewIndex]=clusterNumber;
             }
           }
         } //  while (!pixelList.isEmpty() )
    if (DEBUG_LEVEL>2) {
   System.out.println("Cluster="+clusterNumber+", n="+pixelList.size()+" x="+x+" y="+y);
    }
	 if (clusterSize>0) { // 0 - leave as is
           if (pixelList.size()>clusterSize) { // shrink
             while (pixelList.size()>clusterSize) {
               i=0;
               f=pixels[pixelList.get(i)];
               for (j=1;j<pixelList.size();j++) if (pixels[pixelList.get(j)]<f){
                 i=j;
                 f=pixels[pixelList.get(j)];
               }
               clusterMap[pixelList.get(i)]=-1; // Do not use looking for the next cluster
               pixelList.remove(i);
             }
           } else if (pixelList.size()<clusterSize) { // expand
             while (pixelList.size()<clusterSize) {
               first=true;
               f=0.0f;
               NextIndex=0;
               for (j=0;j<pixelList.size();j++) {
                 Index= pixelList.get(j);
                 for (i=0;i<dirs.length;i++){
                   NewIndex=Index+dirs[i];
                   if ((NewIndex>=0) && (NewIndex<clusterMap.length) && (clusterMap[NewIndex]==0) && (first || (pixels[NewIndex]>f))) {
                     f=pixels[NewIndex];
                     NextIndex=NewIndex;
                     first=false;
                   }
                 }
               }
               pixelList.add (NextIndex);
               clusterMap[NextIndex]=clusterNumber;
             }
           }
/* calculate centroid */
           cx=0.0; cy=0.0; cm=0.0;
           for (i=0;i<pixelList.size();i++) {
             j=pixelList.get(i);
             iy=j/width;
             ix=j-width*iy;
//   System.out.println("j="+j+" x="+ix+" y="+iy);

             f=pixels[j];
             cm+=f;
             cx+=f*ix;
             cy+=f*iy;
           }
           cx/=cm;
           cy/=cm;
           cxy=new Double[2];
           cxy[0]=cx;
           cxy[1]=cy;
           Centroids.add(cxy);
//   System.out.println("New cluster size="+pixelList.size()+" x="+cx+" y="+cy);

         }
         clusterNumber++;
      }
    }
    double[][] coordList=new double[Centroids.size()][2];
    for (i=0;i<coordList.length;i++) {
      coordList[i][0]=Centroids.get(i)[0];
      coordList[i][1]=Centroids.get(i)[1];
    }
//    Double[][] coordList= (Double[][]) Centroids.toArray();
  if (DEBUG_LEVEL>2) {
    for (i=0;i<coordList.length;i++) System.out.println(i+": x="+coordList[i][0]+" y="+coordList[i][1]);
  }
//    return clusterMap;
    return coordList;
  }
////    System.out.println("measureTargets(), N="+N);

/* Normalize pixels values as ratios of difference to average in the surrounding ring to variation in the ring*/
/* TODO: don't roll over, limit */
/// BUG: Seems something wrong - if convolution kernel had DC component - generated all "1.0"

  float [] normalizeAtRing(float [] pixels, int width, int height, boolean[][] mask ) {
    if ((width*height) != pixels.length) {
      IJ.showMessage("Error in  normalizeAtRing","pixels.length ("+pixels.length+") does not match width ("+width+") x height ("+height+") = "+(width*height));
      return null;
    }
    int i,j,x,y,x1,y1,x2,y2, pre_x,pre_y;
    int nFiltPix=0;
    float [] result=new float [width*height];
    double s,s2,d, mean,sigma, meang,sigmag;
    double min=0.0;
    double max=0.0;
    boolean first;
    int ir=(mask.length-1)>>1;
    for (i=0;i<mask.length; i++) for (j=0;j<mask[0].length;j++) if (mask[i][j]) nFiltPix++;

    s= 0.0;
    s2=0.0;
    for (y=0;y<height;y++) for (x=0;x<width;x++) {
      d=pixels[y*width+x];
      s+=d;
      s2+=d*d;
    }
    meang= s/(width+height);
    sigmag=Math.sqrt(s2/(width+height)-meang*meang);

    for (y=0;y<height;y++) for (x=0;x<width;x++) {
      s= 0.0;
      s2=0.0;
      pre_y=y+ir+height; // preparing for "%", making sure it will be positive 
      pre_x=x+ir+width; // preparing for "%", making sure it will be positive 
      first=true;
      for (y1=0;y1 < mask.length; y1++) {
        y2=(pre_y-y1)%height;
        for (x1=0;x1<mask[0].length;x1++) if (mask[y1][x1]) {
          x2=(pre_x-x1)%width;
          d=pixels[y2*width+x2];
          if (first) {
            min=d;
            max=d;
            first=false;
          }
          if (d>max) max=d;
          if (d<min) min=d;
          s+=d;
          s2+=d*d;
        }
      }
      mean= s/nFiltPix;
      sigma=Math.sqrt(s2/nFiltPix-mean*mean);
      sigma=Math.sqrt(sigma*sigmag); // average with image-global sigma
      
//      mean=max;
//      sigma=max-min;
      d=pixels[y*width+x];
      if (sigma>0) { // should always be so
        result[y*width+x]= (float) ((d-mean)/sigma);
      } else {
        result[y*width+x]=1.0f; // any number?
      }
    }
    return result;
  }

/* Convolve image (one Bayer slice) with the inverted target kernel
    Center should be at size/2, size/2 - will convolve only (size-1)*(size-1) */
/**Which is faster - double or float? Double i TWICE faster!*/
  float [] doubleConvolveWithTarget(float [] pixels, float [] kernel_full, int width, int height, int size) {
    int hsize=size/2;
    double [] kernel;
    
    if ((width*height) != pixels.length) {
      IJ.showMessage("Error","pixels.length ("+pixels.length+") does not match width ("+width+") x height ("+height+") = "+(width*height));
      return null;
    }
    if ((size*size) != kernel_full.length) {
      IJ.showMessage("Error","kernel.length ("+kernel_full.length+") does not match size ("+size+") ^2  = "+(size*size));
      return null;
    }
    int i,j;
    double d; // is float faster than double? or opposite (then it makes sesne to convert everything to double first
/* if kernel has even dimensions - ignore first (0) row and first (0) column */

    if ((size & 1)!=0) {
      kernel= new double[size*size];
      for (i=0;i<kernel.length;i++) kernel[i]=kernel_full[i];
    } else {
      size-=1;
      hsize-=1;
      d=0.0f;
      kernel= new double[size*size];
      for (i=0;i<size;i++) for (j=0;j<size;j++) {
        kernel[i*size+j]=kernel_full[(i+1)*(size+1)+(j+1)];
        d+=kernel[i*size+j];
      }
      d/=size*size;
//    System.out.println("Subtracting average value ("+d+") from the convolution kernel");
      for (i=0;i<kernel.length;i++) kernel[i]-=d;
    }
    double [] dPixels=new double[pixels.length];
    for (i=0;i<pixels.length;i++) dPixels[i]=pixels[i];
    
    if (DEBUG_LEVEL>10) IJ.showMessage("Debug doubleConvolveWithTarget","pixels.length="+pixels.length+"\nwidth="+width+"\nheight="+height+"\nkernel.length="+kernel.length+"\nsize="+size);
    float [] result=new float [width*height]; /* this is still float - one conversion on tghe output*/
    int x,y,x1,y1, x2, y2, pre_y,pre_x;
//    double d;
    boolean yMiddle=false;
    int index_kernel, index_source;
    for (y=0;y<height;y++) {


/**/
      yMiddle= (y>=hsize) && (y<(height-hsize));
      if (yMiddle) { // calculate faster when no need to care about borders
        for (x=hsize;x<width-hsize;x++) {
          d=0;
          index_kernel=0;
          index_source=(y+hsize)*width+x+hsize;
          for (y1=0;y1<size;y1++) {
            for (x1=0;x1<size;x1++) {
//   if (index_source<0)   System.out.println("index_source="+index_source+" index_kernel="+index_kernel+" x="+x+" y="+y+" x1="+x1+" y1="+y1);
              d+=dPixels[index_source--]*kernel[index_kernel++]; ///out of bounds: -834
            }
            index_source-=width-size; 
          }
          result[y*width+x]= (float) d;
        }
      }
/**/
/* now finish beginnings/ends of the middle lines and process complete first/last lines*/
      pre_y=y+hsize+height; // preparing for "%", making sure it will be positive 
      for (x=0;x<width;x++) if ((x<hsize) || (x>=(width-hsize)) || !yMiddle){
        d=0;
        pre_x=x+hsize+width; // preparing for "%", making sure it will be positive 
        for (y1=0;y1<size;y1++) {
          y2=(pre_y-y1)%height;
          for (x1=0;x1<size;x1++) {
            x2=(pre_x-x1)%width;
            d+=dPixels[y2*width+x2]*kernel[y1*size+x1];
          }
          result[y*width+x]= (float) d;
        }
      }
    }
    return result;
  } 
  





  public void processWindowEvent(WindowEvent e) {
    super.processWindowEvent(e);
    if (e.getID()==WindowEvent.WINDOW_CLOSING) {
      instance = null;	
    }
  }

  public boolean showDialog() {
    int i;
    GenericDialog gd = new GenericDialog("Target Points parameters");
    gd.addStringField ("Filename prefix:                   ", jp4_instance.getTitle(), 20);
    gd.addNumericField("FFT_Size:                          ", FFTSize, 0);
//    gd.addNumericField("Target minimal outer diameter (pix)", targetOuterDMin, 2);
//    gd.addNumericField("Target maximal outer diameter (pix)", targetOuterDMax, 2);
//    gd.addNumericField("Number of target rings             ", numTargetRings, 0);
//    gd.addNumericField("Subdivide pixels for target generation ", pixelsSubdivide, 0);

    gd.addNumericField("Filtered radius (pix)             ", filteredRadius,   2); //3;  //pix - search maximums after convolution in (2*filteredRadius+1) squares
    gd.addNumericField("Background radius (pix)           ", backgroundRadius, 3); //25; //pix - consider ring area between backgroundRadius and filteredRadius as reference
    gd.addNumericField("Cluster threshold                 ", clusterThreshold, 3); //1.5
    gd.addNumericField("Cluster size (pix)                ", clusterSize, 0); //20

    gd.addNumericField("Target discriminator angular freq.   ", discrAngularFreq,  0); //2 ;  // pixels on FFT image of tragets converted polar (the smaller, the less angular variations)
    gd.addNumericField("Target discriminator radial min freq ", discrRadialMinFreq, 0); //10 ; // pixels on FFT image of tragets converted polar (radial component)
    gd.addNumericField("Target discriminator radial max freq ", discrRadialMaxFreq, 0); //10 ; // pixels on FFT image of tragets converted polar (radial component)
    gd.addNumericField("Target discriminator threshold    ", discrThreshold, 3); //0.3; // FFT energy fraction in selecter area should be > this threshold to pass the test
    gd.addNumericField("Max chromatic aberration (pix)    ", maxChromaticDistance, 1); //10.0; // Maximal distance between the same target on different color copmponents


    gd.addNumericField("Debug Level:                       ", MASTER_DEBUG_LEVEL, 0);


    gd.showDialog();
    if (gd.wasCanceled()) return false;
    jp4_instance.setTitle(gd.getNextString());
    FFTSize=1;
    for (i=             (int) gd.getNextNumber(); i >1; i>>=1) FFTSize <<=1; /* make FFTSize to be power of 2*/
//    targetOuterDMin =         gd.getNextNumber(); // minimal outer diameter of the target image , in pixels
//    targetOuterDMax =         gd.getNextNumber(); // maximal outer diameter of the target image , in pixels
//    numTargetRings  =   (int) gd.getNextNumber();  // number of target black rings (notg counting center black circle)
//    pixelsSubdivide  =   (int) gd.getNextNumber();  // Subdivide pixels for target generation

    filteredRadius   =      gd.getNextNumber();
    backgroundRadius =      gd.getNextNumber();
    clusterThreshold=       gd.getNextNumber();
    clusterSize=      (int) gd.getNextNumber();

    discrAngularFreq=    (int) gd.getNextNumber();
    discrRadialMinFreq=  (int) gd.getNextNumber();
    discrRadialMaxFreq=  (int) gd.getNextNumber();
    discrThreshold=            gd.getNextNumber();
    maxChromaticDistance=      gd.getNextNumber();
    MASTER_DEBUG_LEVEL= (int) gd.getNextNumber();
    return true;
  }

  public float []createTargetDialog() {
    int i;
    GenericDialog gd = new GenericDialog("Target template parameters");
    gd.addNumericField("FFT_Size:                          ", FFTSize, 0);
    gd.addNumericField("Target minimal outer diameter (pix)", targetOuterDMin, 2);
    gd.addNumericField("Target maximal outer diameter (pix)", targetOuterDMax, 2);
    gd.addNumericField("Number of target rings             ", numTargetRings, 0);
    gd.addNumericField("Subdivide pixels for target generation ", pixelsSubdivide, 0);
    gd.addNumericField("Invert deconvolution if less than",   deconvInvert, 3);

//    gd.addNumericField("Debug Level:                       ", MASTER_DEBUG_LEVEL, 0);


    gd.showDialog();
    if (gd.wasCanceled()) return null;
    FFTSize=1;
    for (i=             (int) gd.getNextNumber(); i >1; i>>=1) FFTSize <<=1; /* make FFTSize to be power of 2*/
    targetOuterDMin =         gd.getNextNumber(); // minimal outer diameter of the target image , in pixels
    targetOuterDMax =         gd.getNextNumber(); // maximal outer diameter of the target image , in pixels
    numTargetRings  =   (int) gd.getNextNumber();  // number of target black rings (notg counting center black circle)
    pixelsSubdivide  =  (int) gd.getNextNumber();  // Subdivide pixels for target generation
    deconvInvert=             gd.getNextNumber();  //0.05; // when FFT component is lass than this fraction of the maximal value, replace 1/z with Z
//    MASTER_DEBUG_LEVEL= (int) gd.getNextNumber();
    return createTarget(FFTSize,pixelsSubdivide,targetOuterDMin,targetOuterDMax,numTargetRings,deconvInvert);
  }

  public float [] createTarget(int size, int subdiv, double DMin, double DMax, int nRings, double deconvInvert) {
   ImageProcessor ip_target;
   FHT fht_target;
   double[][][] fft_target;
   int hsizeP1= (size>>1)+1;
   double [][] dpixels=new double [hsizeP1][hsizeP1];
   double [] rMinIn=   new double[nRings+1];
   double [] rMaxIn=   new double[nRings+1];
   double [] rMinOut=  new double[nRings+1];
   double [] rMaxOut=  new double[nRings+1];
   int i,j,i1,j1,n;
   double x,y,r,ks,ke;
   double subFraction=1.0/(subdiv*subdiv);
   double DCLevel=0.0;
   double a,k,r2,k2;
   if (DMin>DMax) {
     x=DMin;
     DMin=DMax;
     DMax=x;
   }
   for (n=0;n<=nRings;n++) {
     rMinIn[n]= DMin*(n*2  )/(2*(2*nRings+1));
     rMinOut[n]=DMin*(n*2+1)/(2*(2*nRings+1));
     rMaxIn[n]= DMax*(n*2  )/(2*(2*nRings+1));
     rMaxOut[n]=DMax*(n*2+1)/(2*(2*nRings+1));
   }

   for (i=0;i<hsizeP1; i++) for (j=0;j<hsizeP1; j++) {
     dpixels[i][j]=0.0;
     for (i1=0;i1<subdiv; i1++) for (j1=0;j1<subdiv; j1++) {
       x=j+0.1*j1;
       y=i+0.1*i1;
       r=Math.sqrt(x*x+y*y);
       for (n=0;n<=nRings;n++) if ((rMinIn[n] <= r)&& (rMaxOut[n]>r)){
         if (rMaxIn[n]>r)   ke=(r-rMinIn[n])/(rMaxIn[n]-rMinIn[n]);
         else ke=1.0;
         if (rMinOut[n]<=r) ks=(r-rMinOut[n])/(rMaxOut[n]-rMinOut[n]);
         else ks=0.0;
///         dpixels[i][j]+=subFraction*(ke-ks);
         dpixels[i][j]-=subFraction*(ke-ks);
       }
     }
     r=dpixels[i][j];
// some piuxels will appear once, some - twice, most - four times
     if ((i>0) &&(i<(size>>1))) r*=2.0; 
     if ((j>0) &&(j<(size>>1))) r*=2.0;
     DCLevel+=r;
   }
   DCLevel/=(size*size);
   for (i=0;i<hsizeP1; i++) for (j=0;j<hsizeP1; j++)  dpixels[i][j]-=DCLevel;

   ip_target = new FloatProcessor(FFTSize,FFTSize);
   for (i=0;i<size; i++) for (j=0;j<size; j++) {
//     ip_target.putPixelValue(j,i, (float) dpixels[(i>=hsizeP1)?(size-i):i][(j>=hsizeP1)?(size-j):j]);
     ip_target.putPixelValue(j ^ (size>>1),i ^ (size>>1), (float) dpixels[(i>=hsizeP1)?(size-i):i][(j>=hsizeP1)?(size-j):j]);
   }
   ip_target.resetMinAndMax();
   if (DEBUG_LEVEL>5) {
     ImagePlus imp_target=  new ImagePlus("Target_direct_"+deconvInvert, ip_target);
     imp_target.show();
   }

   fht_target =  new FHT(ip_target);
// Swapping quadrants, so the center will be 0,0
   fht_target.swapQuadrants();
// get to frequency
   fht_target.transform();
   float [] fht_target_pixels=(float []) fht_target.getPixels();
   
   if (DEBUG_LEVEL>5) {
     ImageProcessor ip_fht_target = new FloatProcessor(size,size);
     ip_fht_target.setPixels(fht_target_pixels);
     ip_fht_target.resetMinAndMax();
     ImagePlus imp_fht_target= new ImagePlus("FHT_"+deconvInvert, ip_fht_target);
     imp_fht_target.show();
   }


// Convert from FHT to complex FFT
   fft_target= FHT2FFTHalf (fht_target);
/* */

/// deconvInvert
/// Now tricky thing. Invert Z for large values, but make them Z - for small ones. So it will be a mixture of correlation and deconvolution
// here the targets are round, but what will the the correct way fo assymmetrical ones?


/// First - find maximal value
   
//   double[][][] fft_target;
   double fft_max=0;
   for (i=0;i<fft_target.length; i++) for (j=0;j<fft_target[0].length;j++) {
     r2=fft_target[i][j][0]*fft_target[i][j][0]+fft_target[i][j][1]*fft_target[i][j][1];
     if (r2>fft_max) fft_max=r2;
   }
   k=Math.sqrt(fft_max)*deconvInvert;
   k2=k*k;

   for (i=0;i<fft_target.length; i++) for (j=0;j<fft_target[0].length;j++) {
     r=Math.sqrt(fft_target[i][j][0]*fft_target[i][j][0]+fft_target[i][j][1]*fft_target[i][j][1]);
     a=-Math.atan2(fft_target[i][j][1],fft_target[i][j][0]); /// will be zero for these targets)
     r=r/(r*r+k2);
     fft_target[i][j][0]=r*Math.cos(a);
     fft_target[i][j][1]=r*Math.sin(a);
   }   

// Convert fft array back to fht array
/**/

   fht_target_pixels= FFTHalf2FHT (fft_target);
// set fht_target pixels with new values 
   fht_target.setPixels (fht_target_pixels);
/// optionally show the result
   if (DEBUG_LEVEL>5) {
     ImageProcessor ip_fht_target1 = new FloatProcessor(size,size);
     ip_fht_target1.setPixels(fht_target_pixels);
     ip_fht_target1.resetMinAndMax();
     ImagePlus imp_fht_target1= new ImagePlus("Inverted_FHT_"+deconvInvert, ip_fht_target1);
     imp_fht_target1.show();
   }
/// transform 

   fht_target.inverseTransform();

   fht_target.swapQuadrants();


   fht_target.resetMinAndMax();
//   ImagePlus imp= new ImagePlus(title, ip_fht);
   if (DEBUG_LEVEL>1) {
     ImagePlus imp_target_inverted= new ImagePlus("Inverted_"+deconvInvert, fht_target);
     imp_target_inverted.show();
   }
//   return direct_target;
   return (float[])fht_target.getPixels();
  }




 /* ignore ROI, use whole image */
  public float[][] splitBayer (ImagePlus imp) {
    ImageProcessor ip=imp.getProcessor();
    Rectangle r=new Rectangle(imp.getWidth(),imp.getHeight());
    float [] pixels;
    pixels=(float[])ip.getPixels();    
    if (DEBUG_LEVEL>10) IJ.showMessage("splitBayer","r.width="+r.width+
                              "\nr.height="+r.height+
                              "\nlength="+pixels.length);
    float [][] bayer_pixels=new float[4][pixels.length>>2];
    int x,y,base,base_b,bv;
    int half_height=r.height>>1;
    int half_width=r.width>>1;
    for (y=0; y<half_height; y++) for (bv=0;bv<2;bv++){
      base=r.width*((y<<1)+bv);
      base_b=half_width*y;
      if (bv==0) for (x=0; x<half_width; x++) {
        bayer_pixels[0][base_b]=  pixels[base++];
        bayer_pixels[1][base_b++]=pixels[base++];
      } else  for (x=0; x<half_width; x++) {
        bayer_pixels[2][base_b]=  pixels[base++];
        bayer_pixels[3][base_b++]=pixels[base++];
      }
    }
    return bayer_pixels;
  }

  public void showBayers(float[][] bayer_pixels, int width, int height, String title) {
    int i;
   if (DEBUG_LEVEL>10) IJ.showMessage("showBayers","width="+width+
                              "\nheight="+height+
                              "\nlength="+bayer_pixels[0].length);

    ImageProcessor[] ip= new ImageProcessor[4];
    ImagePlus[]      imp=new ImagePlus[4];
    for (i=0;i<4;i++) {
      ip[i]=new FloatProcessor(width,height);
      ip[i].setPixels(bayer_pixels[i]);
      ip[i].resetMinAndMax();
      imp[i]=  new ImagePlus(title+"_"+i, ip[i]);
      imp[i].show();
    }
  }



  public float[] initHamming(int size) {
    float [] hamming =new float [size*size];
    float [] hamming_line=new float [size];
    int i,j;
    for (i=0; i<size; i++) hamming_line[i]=(float) (0.54-0.46*Math.cos((i*2.0*Math.PI)/size));
    for (i=0; i<size; i++) for (j=0; j<size; j++){
       hamming[size*i+j]=hamming_line[i]*hamming_line[j];
    }
    return hamming;
  }


  public void initDisplay() {
    if ((imp_display==null) || (imp_display.getWidth()!=displayWidth) || (imp_display.getHeight()!=displayHeight)) {
      if (imp_display!=null)     imp_display.close();
      ip_display=   new ColorProcessor (displayWidth,displayHeight);
      imp_display=  new ImagePlus("Target Points", ip_display);
      imp_display.show();

    }
  }


/* converts FHT results (frequency space) to complex numbers of [FFTSize/2+1][FFTSize] */
private double[][][] FHT2FFTHalf (FHT fht) {
   float[] fht_pixels=(float[])fht.getPixels();
   double[][][] fftHalf=new double[(FFTSize>>1)+1][FFTSize][2];
   int row1,row2,col1,col2;

   for (row1=0;row1<=(FFTSize>>1);row1++) {
     row2=(FFTSize-row1) %FFTSize;
     for (col1=0;col1<FFTSize;col1++) {
       col2=(FFTSize-col1) %FFTSize;
       fftHalf[row1][col1][0]=   0.5*(fht_pixels[row1*FFTSize+col1] + fht_pixels[row2*FFTSize+col2]);
       fftHalf[row1][col1][1]=   0.5*(fht_pixels[row2*FFTSize+col2] - fht_pixels[row1*FFTSize+col1]);
     }
   }
   return fftHalf;
}

/* converts FFT arrays of complex numbers of [FFTSize/2+1][FFTSize] to FHT arrays */
private float[] FFTHalf2FHT (double [][][] fft) {
   float[] fht_pixels=new float [FFTSize*FFTSize];
   int row1,row2,col1,col2;
   for (row1=0;row1<=(FFTSize>>1);row1++) {
     row2=(FFTSize-row1) %FFTSize;
     for (col1=0;col1 < FFTSize;col1++) {
       col2=(FFTSize-col1) %FFTSize;
/* out of bounds */
       fht_pixels[row1*FFTSize+col1]=(float)(fft[row1][col1][0]+fft[row1][col1][1]);
       fht_pixels[row2*FFTSize+col2]=(float)(fft[row1][col1][0]-fft[row1][col1][1]);
     }
   }
   return fht_pixels;
}
	/**
 * Main method for debugging.
 *
 * For debugging, it is convenient to have a method that starts ImageJ, loads an
 * image and calls the plugin, e.g. after setting breakpoints.
 * Grabbed from https://github.com/imagej/minimal-ij1-plugin
 * @param args unused
 */
public static void main(String[] args) {
	// set the plugins.dir property to make the plugin appear in the Plugins menu
	Class<?> clazz = Aberration_Calibration.class;
	String url = clazz.getResource("/" + clazz.getName().replace('.', '/') + ".class").toString();
	String pluginsDir = url.substring(5, url.length() - clazz.getName().length() - 6);
	System.setProperty("plugins.dir", pluginsDir);
	// start ImageJ
	new ImageJ();
	// run the plugin
	IJ.runPlugIn(clazz.getName(), "");
}
}