// Direct Fourier Reconstruction
// Jean-Christophe Taveau
// http://crazybiocomputing.blogspot.com

// 1- Get sinogram and create the final stack
var imp = IJ.getImage();
var ip=imp.getProcessor();
var w=imp.getWidth();
var nbProj=imp.getHeight();
var step=180/nbProj;
var orgx = w/2;
var orgy = w/2;

// Output complex image of sinogram
var FTsino = IJ.createImage("Complex of sinogram", "32-bit Black", w, nbProj, 2);
FTsino.getStack().setSliceLabel("Real",1);
FTsino.getStack().setSliceLabel("Imaginary",2);
var FTsino_p = FTsino.getStack();

//
// M A I N 
//

IJ.log("Part I: 1D-FFT(sinogram) i.e. stack of 1D-FFT(row)");
for (var y=0;y<nbProj;y++)
{
  IJ.log("Slice #"+y);
  // 1- Extract a 1D-projection (one row of sinogram)
  var row = getRow(ip,y);

  // 2- Init complex numbers of this row
  var tmp = new Array(row.length);
  for (var j in row)
    tmp[j] = {"re":row[j], "im": 0.0};

  // 3- 1D-DFT
  row = DFT(tmp,false);

  // 4- Swap data
  for (var j=0;j<w;j++)
  {
    var r = (j + w - w/2) % w;
    if (r % 2 == 0) {
      FTsino_p.setVoxel(r,y,0,row[j].re);
      FTsino_p.setVoxel(r,y,1,row[j].im);
    }
    else {
      FTsino_p.setVoxel(r,y,0,-row[j].re);
      FTsino_p.setVoxel(r,y,1,-row[j].im);
    }
  }    

}

FTsino.show();

// F U N C T I O N S

// Slow Fourier Transform
// From http://arachnoid.com/signal_processing/dft.html
function DFT(input,inverse)
{
  var N = input.length;
  var dft = new Array(N);
  var pi2 = (inverse)? 2.0 * Math.PI:-2.0 * Math.PI;
  var a=0.0;
  var invs = 1.0 / N;
  for(var y = 0;y < N;y++) {
    dft[y] = {"re":0.0, "im":0.0};
    for(var x = 0;x < N;x++) {
      a = pi2 * y * x * invs;
      dft[y].re += input[x].re * Math.cos(a) - input[x].im * Math.sin(a);
      dft[y].im += input[x].re * Math.sin(a) + input[x].im * Math.cos(a);
    }

    if(!inverse) {
      dft[y].re *= invs;
      dft[y].im *= invs;
    }

  }
  return dft;
}