import java.awt.*;
import java.awt.event.*;
import java.io.*;
import ij.plugin.frame.*;
import ij.*;
import ij.process.*;
import ij.gui.*;

/* Author: Giovanni Di Domenico,  in ImageJ
 Contact: giovanni.di.domenico [at] unife.it
 First version: 2015/07/22
 Licence: Public Domain */

public class Focal_Spot  extends Dialog implements ActionListener
{  
  private int    nAngularStep   = 360;    // number of angular step
  private int    filterType     =   0;    // ramp filter
  private int    filterOrder    =   1;    // only for butterworth filter 
  private double axisOffset     = 32.5;   // (Ns+1)/2
  private double filterCutoff   =  0.5;   // frequency cutoff for filter
  private double holeDiameterMm = 10.0;   // hole diameter in mm
  private double pixelSizeMm    = 0.096;   // detector pixel size in mm  
  private static Dialog instance = null;
  
  public Focal_Spot()
  {
    super(new Frame(), "Focal Spot");
    if (instance != null) {
      instance.toFront();
      return;
    }
    instance = null;
    
    if (IJ.versionLessThan("1.21a")) {
      return;
    }
    
    doDialog();
  }
  
  /**
   * Build the dialog box
   *
   */
  private GridBagLayout  layout;
  private GridBagConstraints  constraint;
  private TextField  txtAngStep;      // number of Angula Step to sample the circle
  private TextField  txtAxisOffset;   // s-axis offset
  private TextField  txtCutoff;       // frequency cutoff
  private TextField  txtOrder;        // filter order (Butterworth)
  private TextField  txtHoleDiameter; // hole diameter
  private TextField  txtPixelSize;    // detector pixel size
  private Choice     choiceFilter;    // kind of filter 
  
  private List    lstReport;
  private Button   bnRunSession;
  private Button   bnClose;
  
  private void doDialog()
  {
    // Layout
    layout = new GridBagLayout();
    constraint = new GridBagConstraints();
    
    // input data
    txtAngStep      = new TextField("360",  9);
    txtAxisOffset   = new TextField("32.5",9); //1-4-2025 modified the default value in 32.5
    txtCutoff       = new TextField("0.5", 9);
    txtOrder        = new TextField("1",9);
    txtHoleDiameter = new TextField("10.0",9);
    txtPixelSize    = new TextField("0.096",9);
    
    // filter choice
    choiceFilter  = new Choice();
    choiceFilter.add("Ramp");
    choiceFilter.add("Hann");
    choiceFilter.add("Hamming");
    choiceFilter.add("Butterworth");
    choiceFilter.select(0);    
    
    // List report
    lstReport = new List(15);
    
    // Buttons
    bnClose  = new Button("Close");
    bnRunSession  = new Button("Run");
    
    // Panel buttons 
    Panel pnButtons = new Panel();
    pnButtons.setLayout(new FlowLayout(FlowLayout.CENTER, 8, 0));
    pnButtons.add(bnClose);
    pnButtons.add(bnRunSession);
    
    // Panel filtro 
    //Panel pnFiltro = new Panel();
    //pnFiltro.setLayout(new FlowLayout(FlowLayout.CENTER, 8, 0));
    //pnFiltro.add(new Label("Filter Type"));
    //pnFiltro.add(choiceFilter);
    //pnFiltro.add(new Label("Filter Order"));
    //pnFiltro.add(txtOrder);
    //pnFiltro.add(new Label("Filter Cutoff"));
    //pnFiltro.add(txtCutoff);
    
    // Panel main 
    Panel pnMain = new Panel();
    pnMain.setLayout(layout);
    addComponent(pnMain, 0, 0, 1, 1, 5, new Label("Angular Step"));
    addComponent(pnMain, 0, 1, 1, 1, 5, txtAngStep);
    addComponent(pnMain, 1, 0, 1, 1, 5, new Label("Axis Offset"));
    addComponent(pnMain, 1, 1, 1, 1, 5, txtAxisOffset);
    addComponent(pnMain, 2, 0, 1, 1, 5, new Label("Filter Type"));
    addComponent(pnMain, 2, 1, 1, 1, 5, choiceFilter);
    addComponent(pnMain, 3, 0, 1, 1, 5, new Label("Filter Cutoff"));
    addComponent(pnMain, 3, 1, 1, 1, 5, txtCutoff);
    addComponent(pnMain, 4, 0, 1, 1, 5, new Label("Filter Order"));
    addComponent(pnMain, 4, 1, 1, 1, 5, txtOrder);
    addComponent(pnMain, 5, 0, 1, 1, 5, new Label("Circular Hole Diameter [mm]"));
    addComponent(pnMain, 5, 1, 1, 1, 5, txtHoleDiameter);
    addComponent(pnMain, 6, 0, 1, 1, 5, new Label("Detector Pixel Size [mm]"));
    addComponent(pnMain, 6, 1, 1, 1, 5, txtPixelSize);

    // add Filter panel
    //addComponent(pnMain,  3, 0, 4, 1, 5, pnFiltro);
    // add list report
    addComponent(pnMain, 7, 0, 4, 1, 5, new Label("Report"));
    addComponent(pnMain, 8, 0, 4, 1, 5, lstReport);
    // add button panel
    addComponent(pnMain, 9, 0, 4, 1, 5, pnButtons);
    
    // Add Listeners
    bnClose.addActionListener(this);
    bnRunSession.addActionListener(this);
    
    // Building the main panel
    add(pnMain);
    pack();
    setResizable(false);
    GUI.center(this);
    setVisible(true);
  }
  
  /**
   * Add a component in a panel in the northeast of the cell.
   */
  final private void addComponent(
                                  final Panel pn,
                                  final int row, final int col, 
                                  final int width, final int height, 
                                  final int space,
                                  final Component comp){
    constraint.gridx = col;
    constraint.gridy = row;
    constraint.gridwidth = width;
    constraint.gridheight = height;
    constraint.anchor = GridBagConstraints.NORTHWEST;
    constraint.insets = new Insets(space, space, space, space);
    constraint.weightx = IJ.isMacintosh()?90:100;
    constraint.fill = constraint.HORIZONTAL;
    layout.setConstraints(comp, constraint);
    pn.add(comp);
  }
  
  /**
   * Implements the actionPerformed for the ActionListener.
   */
  public synchronized  void actionPerformed(ActionEvent e) 
  {
    if (e.getSource() == bnClose) {
      dispose();
      instance = null;
    }
    else if (e.getSource() == bnRunSession) {
      run();
    }    
    notify();
  }
  
  /**
   * Get a double value from a TextField between minimal and maximal values.
   */
  private double getDoubleValue(TextField text, double mini, double defaut, double maxi)
  {
    double d;
    try {
      d = (new Double(text.getText())).doubleValue();
      if (d < mini)  
        text.setText( "" + mini);
      if (d > maxi)  
        text.setText( "" + maxi);
    }
    catch (Exception e) {
      if (e instanceof NumberFormatException) 
        text.setText( "" + defaut);
    }
    d = (new Double(text.getText())).doubleValue();
    return d;
  } 
  
  /**
   * Get a integer value from a TextField between minimal and maximal values.
   */
  private int getIntegerValue(TextField text, int mini, int defaut, int maxi)
  {
    int d;
    try {
      d = (new Integer(text.getText())).intValue();
      if (d < mini)  text.setText( "" + mini);
      if (d > maxi)  text.setText( "" + maxi);
    }
    catch (Exception e) {
      if (e instanceof NumberFormatException) 
        text.setText( "" + defaut);
    }
    d = (new Integer(text.getText())).intValue();
    return d;
  }
  
  /**
   * Implementation of the method run of the interface class PlugInFilter.
   */
  private void run()
  {
    ImagePlus imp = WindowManager.getCurrentImage(); // image to analyze
    ImageProcessor ip = imp.getProcessor();          //
    ImageAccess ima = new ImageAccess(ip);           // 
    int nx = imp.getWidth();
    int ny = imp.getHeight();
    
    /*********************************************************************/
    nAngularStep = getIntegerValue(txtAngStep,0,360,1024);
    axisOffset   = getDoubleValue(txtAxisOffset,0.0,32.5,1024.0);//1-4-2025 modified the default value in 32.5
    filterCutoff = getDoubleValue(txtCutoff,0.0,0.5,0.5);
    filterOrder  = getIntegerValue(txtOrder,0, 1, 50);
    holeDiameterMm = getDoubleValue(txtHoleDiameter,0.0,10.0,100.0);
    pixelSizeMm = getDoubleValue(txtPixelSize,0.0,0.096,10.0);
    
    String selezione;
    // Tipo di filtro
    selezione = (String)choiceFilter.getSelectedItem();
    if( selezione == "Ramp" ) {
      filterType = 0;
    }
    else if ( selezione == "Hann" ) {
      filterType = 1;
    }
    else if ( selezione == "Hamming" ) {
      filterType = 2;
    }
    else if ( selezione == "Butterworth" ) {
      filterType = 3;
    }
    
    // Write the parameters list on tne report window
    lstReport.add("---- STARTING FOCAL SPOT Plugin ----");
    lstReport.add( "ANGULAR STEP  = " + nAngularStep);
    lstReport.add( "FILTER TYPE   = " + filterType);
    lstReport.add( "FILTER CUTOFF = " + filterCutoff);
    lstReport.add( "FILTER ORDER  = " + filterOrder);
    lstReport.add( "AXIS OFFSET   = " + axisOffset );
    lstReport.add( "HOLE DIAMETER [mm]  = " + holeDiameterMm);
    lstReport.add( "PIXEL SIZE    [mm]  = " + pixelSizeMm );
    /*************************************************************************/
    //
    //  check the roi selection
    //
    Roi roi = imp.getRoi();
    if (roi!=null && !roi.isArea()) roi = null;
    ImageProcessor mask =  roi!=null ? roi.getMask():null;
    Rectangle r = roi!=null ? roi.getBounds():new Rectangle(0,0,ip.getWidth(),ip.getHeight());
    //IJ.log("r.x: " + r.x + " r.y: " + r.y + " r.width: " + r.width + " r.height: " + r.height);
    //
    // calculate the average on the rectangular roi selection
    //
    double avg_fore = 0.0,threshold=0.0;
    int count = 0;
    for (int y=0; y<r.height; y++) {
      for (int x=0; x<r.width; x++) {
        count++;
        avg_fore += ip.getPixelValue(x+r.x, y+r.y);
      }
    }
    avg_fore /= count;
    //IJ.log("avg_fore: " + avg_fore);
    
    // 18-11-2019 
    // background estimation for image with background value
    // I use a roi  at top-left part of image (10x10)
    // 
    double avg_bgnd = 0.0;
    count = 0;
    for (int y=1; y<=10; y++) {
      for (int x=1; x<=10; x++) {
        count++;
        avg_bgnd += ip.getPixelValue(x, y);
      }
    }
    avg_bgnd /= count;
    
    // 18-11-2019
    // the threshold value is estimated respect to the background
    //
    threshold = avg_bgnd + (avg_fore - avg_bgnd)/2.0;
    //
    // build a horizontal selection large as the image 
    //
    double[] hori  = new double[nx];
    double[] xLeft = new double[r.height];
    double[] xRight= new double[r.height];
    double[] aa = new double[2];  
    for(int y = 0;y  < r.height;y++)
    {
      ima.getRow(r.y + y,hori);
      scanLine(hori,threshold,aa);
      //IJ.log("a : " + aa[0] + " b : " + aa[1]);
      xLeft[y]  = aa[0];
      xRight[y] = aa[1];
    }
    //
    //  build a vertical selection height as the image 
    //
    double[] vert = new double[ny];
    double[] yUp = new double[r.width];
    double[] yDown = new double[r.width];
    for(int x = 0;x  < r.width;x++)
    {
      ima.getColumn(r.x + x,vert);
      scanLine(vert,threshold,aa);
      //IJ.log("a : " + aa[0] + " b : " + aa[1]);
      yDown[x] = aa[0];
      yUp[x] = aa[1];
    }
    
    //
    // calculation of x_c, y_c and the radius along the horizontal and vertical direction 
    //
    double x_c = 0.0;
    double r_x = 0.0;
    for(int i = 0; i < r.height;i++)
    {
      x_c += (xRight[i]+xLeft[i])/2.0;
      r_x += (xRight[i]-xLeft[i])/2.0; 
    }
    x_c = Math.round(x_c/r.height);
    r_x = Math.round(r_x/r.height);
    
    double y_c = 0.0;
    double r_y = 0.0;
    for(int i = 0; i < r.width;i++)
    {
      y_c += (yUp[i]+yDown[i])/2.0;
      r_y += (yUp[i]-yDown[i])/2.0;
    }
    y_c = Math.round(y_c/r.width);
    r_y = Math.round(r_y/r.width);
    double magnification = 2.0*Math.max(r_x,r_y) * pixelSizeMm / holeDiameterMm;
    
    //
    // write the information about the hole center, radius and magnification
    //
    lstReport.add("---- OUTPUT FOCAL SPOT Plugin ----");
    lstReport.add("Hole Center X: " + Math.round(x_c));
    lstReport.add("Hole Center Y: " + Math.round(y_c));
    lstReport.add("Hole Radius X [pix]: " + Math.round(r_x));
    lstReport.add("Hole Radius Y [pix]: " + Math.round(r_y));
    lstReport.add("Magnification: " + magnification); 
    // 
    // sinogram calculation using the idea of radial reslice plugin (author: Julian Cooper). 
    // Starting from the (x_c,y_c) point the plugin calcculate the image profile
    // along nAngularStep lines covering 360° degrees.
    //
    double arcAngle = 360.0;
    double a = (Math.PI/180.0)*(arcAngle/nAngularStep);
    double XO,X1,YO,Y1;
    double ang;
	//
	// 1/04/2025 Modified the code to manage sinograms height less than 128
	//
	int radius = (int)Math.max(r_x,r_y);
	int sino_height = 64;
	if(radius > 64)
	{
		sino_height = 128;
		lstReport.add("Sinogram height changes: modify Axis Offset.");
	}
	int radius_plus = radius + sino_height/2;
	int radius_minus = radius - sino_height/2;
    lstReport.add("Radius:" + radius + ", radius_plus:" + radius_plus + ", radius_minus:" + radius_minus);
    
    ImageAccess sinogram = new ImageAccess(nAngularStep,radius_plus);
    double[] col = new double[radius_plus];
    XO = x_c; YO = y_c;
    //X1 = x_c + radius;Y1 = YO;
    for(int i = 0;i < nAngularStep;i++)
    {
      ang = i * a;
      X1 = (XO + radius_plus*Math.cos(ang));
      Y1 = (YO - radius_plus*Math.sin(ang));
      getLine(ip,XO,YO,X1,Y1,col);
      drawLine(XO, YO, X1, Y1, imp);
      sinogram.putColumn(i,col);
    }
	//lstReport.add("1st step Sinogram Creation done.");
    //
    // apply the 1st derivative filter along the vertical direction of sinogram
    //
    double[] colin = new double[radius_plus];
    double[] colout = new double[radius_plus];
    for(int ix = 0; ix < nAngularStep;ix++)
    {
      sinogram.getColumn(ix,colin);
      applyFirstDerivativeFilter(colin,colout);
      sinogram.putColumn(ix,colout);
    }
	//lstReport.add("2nd step Sinogram Creation done.");
	//ImagePlus sng0 = new ImagePlus("sinogram",sinogram.createFloatProcessor());
	//sng0.show();
    //sng0.updateAndDraw();
	ImageAccess sino = new ImageAccess(nAngularStep,sino_height);
    sinogram.getSubImage(radius_minus,radius_plus - 1,sino); // 1/04/2025 modified the final sinogram selection
    ImagePlus sng = new ImagePlus("sinogram",sino.createFloatProcessor());
    sng.show();
    sng.updateAndDraw();
	lstReport.add("3rd step Sinogram Visualization done.");
    // 18-11-2019
    // how to estimate the right offset?
    // I try to use the formula sino(phi,s) = sino(phi + pi,-s)
    // Starting from this formula I search the offset to minimize the
    // find_min_of(sum_phi_s(sino(phi,s+offset) - sino(phi + pi,-s+offset))^2)
    //
    double err_min = Double.MAX_VALUE;
    int offset_opt = 0,offset;
    for(int i_off=-5;i_off<=5;i_off++)
    {
      offset = i_off;
      double error = 0.0;
      int ncounts = 0;
	  int ns_element = sino_height/2 - Math.abs(i_off); 
      //IJ.log("i_off: " + i_off + ", ns_element: " + ns_element);
	  for(int iphi = 0; iphi < nAngularStep/2;iphi++)
      {
        for(int is = 0; is < ns_element;is++)
        {
		  // 	
		  // 1/04/2025 Modified code for the correct estimation of rotation axis offset
          // In a sinogram with height equal to NS and axis at center of FOV the corresponding points
          // are is=0 and is = NS-1, is=1 and is=NS-2, is=2 and is=NS-3 and so on,
          // if the axis offset is -1 the the corresponding points are
          // is=0 and is = NS-3, is=1 and is=Ns-4, and so on.
          // the correct rules is is - offset and NS-1 - 2*offset	    		  
		  int is_neg = (sino_height/2 - 1) - is + i_off;
          int is_pos = (sino_height/2 + 0) + is + i_off;
		  //if(iphi == 0)
		  //	IJ.log("is_neg: " + is_neg + ", is_pos: " + is_pos);
          if( (is_neg >= 0) && (is_neg < sino_height) && (is_pos >= 0) && (is_pos < sino_height))
          {
            double scarto = (sino.getPixel(iphi,is_neg)-sino.getPixel(iphi+nAngularStep/2,is_pos));
            error = error + scarto*scarto;
            ncounts++;
          }
        }
      }
      error = Math.sqrt(error/ncounts);
      IJ.log("Stack Image:" +(i_off + 6) + ", axis offset: " + offset + ", error: " + error + ", ncounts: " + ncounts);
      if(error < err_min) 
      {
        err_min = error;
        offset_opt=offset;
      }
    }
    //IJ.log("offset_opt: " + offset_opt + " err_min: " + err_min);
    // 18-11-2019
    // sinogram reconstruction with different value of  axis offset
    // if axisOffset is (sino_height-1)/2 the right image is the number 6
    //
    int nStack = 11;
    double offStep = 0.5,offset_dbl;
    ImageAccess focalSpot;
    ImageStack stack = new ImageStack(sino_height,sino_height);
    for(int i = 0; i < nStack; i++)
    {
      offset_dbl = (axisOffset + offset_opt) + (i - nStack/2)*offStep;
      focalSpot = FBP(sino,sino_height,nAngularStep,filterType,filterOrder,filterCutoff,offset_dbl);
      stack.addSlice("FocalSpot ",focalSpot.createFloatProcessor(),i);
      lstReport.add("Calculated " + i + "-mo stack element");
    }
    ImagePlus fs = new ImagePlus("FocalSpot Stack",stack);
    fs.show();
    fs.updateAndDraw();
  }
  /*******************************************************************************************************************/
  
  /**
   *  Function to calculate the 1st derivative of in array.
   *  The output 
   */ 
  public void applyFirstDerivativeFilter(double[] in, double[] out)
  {
    //
    // il filtro e' [-1 1] quindi
    //
    int n = in.length;
    out[0]   = 0.0;
    out[n-1] = 0.0;
    
    for(int i=1;i < n-1;i++)
      out[i] = in[i-1]-in[i];
  }
  
  public void scanLine(double arr[],double threshold,double[] aa)
  {
    int length = arr.length;
    //
    // scansione del vettore
    //    
    for(int i = 0; i < length;i++)
    {
      if (arr[i] < threshold)
      {
        aa[0] = i;
      }
      else 
      {
        break;
      }
    }
    
    for(int i = length - 1; i > 0;i--)
    {
      if(arr[i] < threshold)
      {
        aa[1] = i;
      }
      else 
      {
        break;
      }
    }
  }
  
  private double[] getLine(ImageProcessor ip, double x1, double y1, double x2, double y2, double[] data) 
  {
    double dx = x2-x1;
    double dy = y2-y1;
    int n = (int)Math.round(Math.sqrt(dx*dx + dy*dy));
    if (data==null)
      data = new double[n];
    double xinc = dx/n;
    double yinc = dy/n;
    double rx = x1;
    double ry = y1;
    for (int i=0; i<n; i++) {
      data[i] = (double)ip.getInterpolatedValue(rx, ry);
      rx += xinc;
      ry += yinc;
    }
    return data;
  } 
  
  void drawLine(double x1, double y1, double x2, double y2, ImagePlus imp) 
  {
    ImageCanvas ic = imp.getCanvas();
    if (ic==null) return;
    Graphics g = ic.getGraphics();
    g.setColor(Roi.getColor());
    g.setXORMode(Color.black);
    g.drawLine(ic.screenX((int)(x1+0.5)), ic.screenY((int)(y1+0.5)), ic.screenX((int)(x2+0.5)), ic.screenY((int)(y2+0.5)));
  }
  /*
   *  BackProjection Method for focal spot reconstruction 
   * 
   */
  public ImageAccess FBP(ImageAccess img,int ns,int nphi,int filtro,int ordine,double cutoff,double sOffset)
  {
    int nx   = ns;
    int ny   = ns;
    /*
     *  (length = 2^k) >= ns
     */
    int length = 1 << (int)(Math.log(ns)/Math.log(2.0) + 1.0);
    if(length != (ns * 2)) length = length * 2; // per applicare la convoluzione con l'uso della FFT
    //IJ.log("ns = " + ns + " length = " + length);
    int idummy,is;
    double deltaphi,deltas,deltax,deltay;
    double phiMin,sMin,xMin,yMin,offset;
    double phi,sPixel,xPixel,yPixel,cphi,sphi;
    double oldvalue,incvalue,newvalue;
    double[] profilo = new double[ns];
    double[] ar      = new double[length];
    double[] ai      = new double[length];
    double[] yfiltro = new double[length];
    double[] xfiltro = new double[length];
    double[] window  = new double[length];
    
    PlotWindow  plotWindow;
    ImageAccess sngFiltered = new ImageAccess(length,nphi);
    ImageAccess imgTmp      = new ImageAccess(nx,ny);
    //
    lstReport.add("length = " + length);
    lstReport.add("width  = " + img.getWidth());
    lstReport.add("height = " + img.getHeight());
    lstReport.add("nphi   = " + nphi);
    
    deltaphi = 0.0;
    phiMin = 0.0;
    deltaphi = 2.0*Math.PI/nphi;
    
    deltas = 1.0;
    deltax = (deltay = deltas);
    
    sMin = -(ns - 1)*deltas/2.0;
    xMin = -(nx - 1)*deltax/2.0;
    yMin = -(ny - 1)*deltay/2.0;
    idummy = (length - ns)/2;
    offset = sOffset + (length -ns)/2.0;
    
    //lstReport.add("idummy = " + idummy);
    lstReport.add("sOffset = " + sOffset);
    //lstReport.add("offset = " + offset);
    
    computeFilterRamLak(xfiltro,yfiltro,deltas,length);
    FS_FFT.computeFFT(1,length,xfiltro,yfiltro);
    computeFilterWindow(filtro, cutoff, ordine,window);
    
    //for(int i = 0; i < length;i++)
    //{
    //  double modulo = window[i]*Math.sqrt(xfiltro[i]*xfiltro[i]+yfiltro[i]*yfiltro[i]);
    //  IJ.log(i + " " + modulo);
    //}
    
    //lstReport.add("start filtering...");
    //ImageAccess imgNew = img.duplicate();
    //ImagePlus impResult = new ImagePlus("pippo", imgNew.createFloatProcessor());
    //impResult.show();
    
    for(int i = 0; i < nphi; i++)
    {
      img.getColumn(i,profilo);
      
      for(int j = 0; j < idummy; j++)
      {
        ar[j] = 0.0;
        ai[j] = 0.0;
      }
      for(int j = 0; j < ns; j++)
      {
        ar[j + idummy] = profilo[j];
        ai[j] = 0.0;
      }
      for(int j = 0; j < idummy; j++)
      {
        ar[j + idummy + ns] = 0.0;
        ai[j] = 0.0;
      }
      
      FS_FFT.computeFFT(1,length,ar,ai);
      double re,im; 
      for(int j = 0; j < length;j++)
      {
        re = window[j]*(ar[j]*xfiltro[j]-ai[j]*yfiltro[j]);
        im = window[j]*(ar[j]*yfiltro[j]+ai[j]*xfiltro[j]);
        ar[j]=re;ai[j]=im;
      }
      
      FS_FFT.computeFFT(-1,length,ar,ai);
      //lstReport.add("FFT-inverse ");
      
      sngFiltered.putRow(i,ar);
      
    }
    
    //sngFiltered.show("sngFiltered");
    
    // retroproiezione
    for(int i = 0; i < nphi;i++) {
      //lstReport.add("iphi = " + i);
      phi = phiMin + i * deltaphi;
      cphi = Math.cos(phi);
      sphi = Math.sin(phi);
      sngFiltered.getRow(i,ar);
      for(int ix = 0; ix < nx; ix++) {
        xPixel = xMin + ix * deltax;
        for(int iy = 0; iy < ny; iy++) {
          yPixel = yMin + iy  * deltay;
          sPixel = -xPixel *sphi + yPixel*cphi;
          is = (int)Math.round((sPixel + offset)/deltas);
          if( (is >= 0) && (is < ar.length-1)) {
            oldvalue = imgTmp.getPixel(ix,iy);
            incvalue = ar[is] + (ar[is + 1] - ar[is])/deltas*(sPixel + offset - is *deltas);
            newvalue = oldvalue + incvalue;
            imgTmp.putPixel(ix,iy,newvalue);
          }
        }
      }
    }
    
    return imgTmp;
  }
  
  /*
   *  method to compute the filter window W(v) -> H(v)=F(v)*W(v)
   */ 
  void computeFilterWindow(int filterType,double cutoff,int ordine,double[] window)
  {
    double deltanu,freq;
    int    n = window.length;
    
    deltanu = 1.0/n;
    
    switch(filterType)
    {
      case 0: // pure RAM-LAK
        window[0]   = 1.0;
        freq = (n/2)*deltanu;
        if( freq <= cutoff ) 
        {
          window[n/2] = 1.0;
        }
        else {
          window[n/2] = 0.0;
        }
        
        for(int i = 1;i < n/2;i++)
        {
          freq = i * deltanu;
          if( freq <= cutoff ) 
          {    
            window[i]     =  1.0;
            window[n - i] =  1.0;
          }
          else {
            window[i]     =  0.0;
            window[n - i] =  0.0;
          }
        }
        break;
      case 1: // Hann
        window[0] = 1.0;
        freq = (n/2)*deltanu;
        window[n/2] = (0.5 + 0.5 * Math.cos(Math.PI * freq/cutoff));
        
        for(int i = 1;i < n/2;i++)
        {
          freq = i * deltanu;
          window[i]     =  (0.5 + 0.5 * Math.cos(Math.PI * freq/cutoff));
          window[n - i] =  (0.5 + 0.5 * Math.cos(Math.PI * freq/cutoff));
        }
        break;
      case 2: //  Hamming
        window[0] = 1.0;
        freq = (n/2)*deltanu;
        window[n/2] = (0.5 + 0.46 * Math.cos(Math.PI * freq/cutoff));
        
        for(int i = 1;i < n/2;i++)
        {
          freq = i * deltanu;
          window[i]     =  (0.5 + 0.46 * Math.cos(Math.PI * freq/cutoff));
          window[n - i] =  (0.5 + 0.46 * Math.cos(Math.PI * freq/cutoff));
        }
        break;
      case 3: //  Butterworth
        window[0] = 1.0;
        freq = (n/2)*deltanu;
        window[n/2] = 1.0/Math.sqrt(1.0 + Math.pow((freq/cutoff),2.0*ordine));
        
        for(int i = 1;i < n/2;i++)
        {
          freq = i * deltanu;
          window[i]     =  1.0/Math.sqrt(1.0 + Math.pow((freq/cutoff),2.0*ordine));
          window[n - i] =  1.0/Math.sqrt(1.0 + Math.pow((freq/cutoff),2.0*ordine));
        }
        break;
    }                  
  }
  
  /*
   * Calcolo il filtro RamLak nello spazio immagine
   */
  void computeFilterRamLak(double[] xfiltro,double[] yfiltro,double deltas,int n)
  {
    double pi2 = Math.PI*Math.PI;
    double d,t2=1.0;
    
    if(deltas > 0.0) t2=1.0/(deltas*deltas);
    
    for(int i = 0; i < n;i++)
    {
      xfiltro[i]=0.0;
      yfiltro[i]=0.0;
    }
    
    xfiltro[0]=0.25*t2;
    
    for(int i = 1;i <= n/2;i=i+2)
    {
      d = -t2/(pi2*i*i);
      xfiltro[i]    =d;
      xfiltro[n - i]=d;
    }     
  }
}

/**
 This class contains the algorithm to perform a discrete FFT (Fast Fourier Transform)
 for real periodic functions. To use it, first choose a sample size N (which <b>must</b> 
 be a power of 2). Divide the period into N equally spaced points and evaluate a function 
 at those N points. Save these function values in a <tt>double array ar</tt>. Create a
 second array <tt>b</tt> of N zeros. Then use 
 
 <blockquote>
 <pre>
 FFT.computeFFT(ar, ai)
 </pre>
 </blockquote>
 
 This will return the FFT coefficients in the arrays <tt>ar</tt> and <tt>ai</tt>.
 <p>
 
 <b>Please note:</b> This class will not compile. Why not? How do you fix it? 
 Also, it is not immediately clear how the output values contained in <tt>ar</tt> 
 and <tt>ai</tt> relate to the Fourier coefficients of the Fourier Series for the
 given function <tt>f</tt>. To find that relationship, analyse the coefficients for
 the functions <tt>sin(t)</tt>, <tt>cos(t)</tt>, and 
 <tt>1 + 2*sin(t) + 3*cos(t) + 4*sin(2t) + 5*cos(2t)</tt>. That should clarify it.
 */

class FS_FFT
{ 
  private final static double EPSI = 10E-5;
  
  public static void computeFFT(int dir, int n, double ar[], double ai[]) 
  { 
    double scale = Math.sqrt(1.0/n);
    double tempr,tempi;
    double delta,w,wr,wi,tr,ti;
    int i, j, m;
    int mmax,istep;
    
    for (i = j = 0; i < n; ++i) 
    { 
      if (j >= i) 
      { 
        tempr = ar[j]*scale;
        tempi = ai[j]*scale;
        ar[j] = ar[i]*scale;
        ai[j] = ai[i]*scale;
        ar[i] = tempr;
        ai[i] = tempi;
      }
      m = n/2;
      while ((m >= 1) && (j >= m)) 
      { 
        j -= m;
        m /= 2;
      }
      j += m;
    }
    
    for (mmax = 1, istep = 2*mmax; mmax < n; mmax = istep, istep = 2*mmax) 
    { 
      delta = dir*Math.PI/(double)mmax;
      for (m = 0; m < mmax; ++m) 
      { 
        w = m*delta;
        wr = Math.cos(w);
        wi = Math.sin(w);
        for (i = m; i < n; i += istep) 
        { 
          j = i + mmax;
          tr = wr*ar[j] - wi*ai[j];
          ti = wr*ai[j] + wi*ar[j];
          ar[j] = ar[i] - tr;
          ai[j] = ai[i] - ti;
          ar[i] += tr;
          ai[i] += ti;
        }
      }
      mmax = istep;
    }
  }
}