/*
* Copyright (C) 2014 - Martin Berger
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
package edu.stanford.rsl.conrad.data.generic.complex;
import ij.IJ;
import ij.ImageJ;
import edu.emory.mathcs.jtransforms.fft.FloatFFT_1D;
import edu.emory.mathcs.jtransforms.fft.FloatFFT_2D;
import edu.emory.mathcs.jtransforms.fft.FloatFFT_3D;
import edu.stanford.rsl.conrad.data.generic.datatypes.Complex;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.utils.ImageUtil;
public class Fourier {
public void fft(ComplexGrid input){
fft(input,0);
}
public void fft(ComplexGrid input, int dim){
if (input.getSize().length <= dim){
throw new IllegalArgumentException("FFT along provided dimension not possible: Exceeds data dimension");
}
FloatFFT_1D fftInstance = new FloatFFT_1D(input.getSize()[dim]);
// If dim == 0, the memory is in right shape and we can transform right away
if (dim == 0){
for (int i = 0; i < input.getNumberOfElements()/input.getSize()[dim]; i++) {
fftInstance.complexForward(input.getAslinearMemory(), i*input.getSize()[dim]*2);
}
input.notifyAfterWrite();
}
else if(input.getSize().length==2){
// Only FFT over columns are possible - Otherwise dim would be 0 and row-wise would be selected.
for (int i = 0; i < input.getSize()[0]; i++){
float[] cplxData = dataToStridedLine(input, new int[]{i,0}, dim);
fftInstance.complexForward(cplxData);
stridedLineToData(input, new int[]{i,0}, dim, cplxData);
}
}
else if(input.getSize().length==3){
// FFTs over columns and depth are possible - Rows would yield dim == 0 --> first case
if (dim == 1){ // columns case
for (int i = 0; i < input.getSize()[0]; i++){
for (int k = 0; k < input.getSize()[2]; k++){
float[] cplxData = dataToStridedLine(input, new int[]{i,0,k}, dim);
fftInstance.complexForward(cplxData);
stridedLineToData(input, new int[]{i,0,k}, dim, cplxData);
}
}
}
else{ // z-case
for (int i = 0; i < input.getSize()[0]; i++){
for (int j = 0; j < input.getSize()[1]; j++){
float[] cplxData = dataToStridedLine(input, new int[]{i,j,0}, dim);
fftInstance.complexForward(cplxData);
stridedLineToData(input, new int[]{i,j,0}, dim, cplxData);
}
}
}
}
}
public void ifft(ComplexGrid input){
ifft(input,0);
}
public void ifft(ComplexGrid input, int dim){
if (input.getSize().length <= dim){
throw new IllegalArgumentException("FFT along provided dimension not possible: Exceeds data dimension");
}
FloatFFT_1D fftInstance = new FloatFFT_1D(input.getSize()[dim]);
// If dim == 0, the memory is in right shape and we can transform right away
if (dim == 0){
for (int i = 0; i < input.getNumberOfElements()/input.getSize()[dim]; i++) {
fftInstance.complexInverse(input.getAslinearMemory(), i*input.getSize()[dim]*2,true);
}
input.notifyAfterWrite();
}
else if(input.getSize().length==2){
// Only FFT over columns are possible - Otherwise dim would be 0 and row-wise would be selected.
for (int i = 0; i < input.getSize()[0]; i++){
float[] cplxData = dataToStridedLine(input, new int[]{i,0}, dim);
fftInstance.complexInverse(cplxData,true);
stridedLineToData(input, new int[]{i,0}, dim, cplxData);
}
}
else if(input.getSize().length==3){
// FFTs over columns and depth are possible - Rows would yield dim == 0 --> first case
if (dim == 1){ // columns case
for (int i = 0; i < input.getSize()[0]; i++){
for (int k = 0; k < input.getSize()[2]; k++){
float[] cplxData = dataToStridedLine(input, new int[]{i,0,k}, dim);
fftInstance.complexInverse(cplxData,true);
stridedLineToData(input, new int[]{i,0,k}, dim, cplxData);
}
}
}
else{ // z-case
for (int i = 0; i < input.getSize()[0]; i++){
for (int j = 0; j < input.getSize()[1]; j++){
float[] cplxData = dataToStridedLine(input, new int[]{i,j,0}, dim);
fftInstance.complexInverse(cplxData,true);
stridedLineToData(input, new int[]{i,j,0}, dim, cplxData);
}
}
}
}
}
public void fft2(ComplexGrid input){
fft2(input,0);
}
public void fft2(ComplexGrid input, int dim){
if (input.getSize().length < 2){
throw new IllegalArgumentException("FFT along provided dimension not possible: Exceeds data dimension");
}
if (input.getSize().length == 2){
FloatFFT_2D fftInstance = new FloatFFT_2D(input.getSize()[1], input.getSize()[0]);
fftInstance.complexForward(input.getAslinearMemory());
input.notifyAfterWrite();
}else{ // 3D case
FloatFFT_2D fftInstance=null;
switch(dim){
case 0: // FFT along x/y planes
fftInstance = new FloatFFT_2D(input.getSize()[1], input.getSize()[0]);
for (int k = 0; k < input.getSize()[2]; k++) {
float[] cplx = new float[input.getSize()[0]*input.getSize()[1]*2];
System.arraycopy(input.getAslinearMemory(), k*cplx.length, cplx, 0, cplx.length);
fftInstance.complexForward(cplx);
System.arraycopy(cplx, 0, input.getAslinearMemory(), k*cplx.length, cplx.length);
}
input.notifyAfterWrite();
break;
case 1: // FFT along x/z planes
fftInstance = new FloatFFT_2D(input.getSize()[2], input.getSize()[0]);
for (int j = 0; j < input.getSize()[1]; j++) {
float[] cplx = dataToPlane(input,j*input.getSize()[0],dim);
fftInstance.complexForward(cplx);
planeToData(input,j*input.getSize()[0],dim,cplx);
}
break;
case 2: // FFT along y/z planes
fftInstance = new FloatFFT_2D(input.getSize()[2], input.getSize()[1]);
for (int i = 0; i < input.getSize()[0]; i++) {
float[] cplx = dataToPlane(input,i,dim);
//new ComplexGrid2D(cplx, 0, input.getSize()[1], input.getSize()[2]).show();
fftInstance.complexForward(cplx);
planeToData(input,i,dim,cplx);
}
break;
default:
break;
}
}
}
public void ifft2(ComplexGrid input){
ifft2(input,0);
}
public void ifft2(ComplexGrid input, int dim){
if (input.getSize().length < 2){
throw new IllegalArgumentException("FFT along provided dimension not possible: Exceeds data dimension");
}
if (input.getSize().length == 2){
FloatFFT_2D fftInstance = new FloatFFT_2D(input.getSize()[1], input.getSize()[0]);
fftInstance.complexInverse(input.getAslinearMemory(),true);
input.notifyAfterWrite();
}else{ // 3D case
FloatFFT_2D fftInstance=null;
switch(dim){
case 0: // FFT along x/y planes
fftInstance = new FloatFFT_2D(input.getSize()[1], input.getSize()[0]);
for (int k = 0; k < input.getSize()[2]; k++) {
float[] cplx = new float[input.getSize()[0]*input.getSize()[1]*2];
System.arraycopy(input.getAslinearMemory(), k*cplx.length, cplx, 0, cplx.length);
fftInstance.complexInverse(cplx,true);
System.arraycopy(cplx, 0, input.getAslinearMemory(), k*cplx.length, cplx.length);
}
input.notifyAfterWrite();
break;
case 1: // FFT along x/z planes
fftInstance = new FloatFFT_2D(input.getSize()[2], input.getSize()[0]);
for (int j = 0; j < input.getSize()[1]; j++) {
float[] cplx = dataToPlane(input,j*input.getSize()[0],dim);
fftInstance.complexInverse(cplx,true);
planeToData(input,j*input.getSize()[0],dim,cplx);
}
break;
case 2: // FFT along y/z planes
fftInstance = new FloatFFT_2D(input.getSize()[2], input.getSize()[1]);
for (int i = 0; i < input.getSize()[0]; i++) {
float[] cplx = dataToPlane(input,i,dim);
fftInstance.complexInverse(cplx,true);
planeToData(input,i,dim,cplx);
}
break;
default:
break;
}
}
}
public void fft3(ComplexGrid input){
FloatFFT_3D fftInstance = new FloatFFT_3D(input.getSize()[2], input.getSize()[1], input.getSize()[0]);
fftInstance.complexForward(input.getAslinearMemory());
input.notifyAfterWrite();
}
public void ifft3(ComplexGrid input){
FloatFFT_3D fftInstance = new FloatFFT_3D(input.getSize()[2], input.getSize()[1], input.getSize()[0]);
fftInstance.complexInverse(input.getAslinearMemory(),true);
input.notifyAfterWrite();
}
private float[] dataToStridedLine(ComplexGrid grid, int[] idx, int dim){
float[] out = new float[grid.getSize()[dim]*2];
for (int i = 0; i < grid.getSize()[dim]; i++) {
out[i*2]= (float)grid.getValue(idx).getReal();
out[i*2+1]=(float)grid.getValue(idx).getImag();
idx[dim]++;
}
return out;
}
private void stridedLineToData(ComplexGrid grid, int[] idx, int dim, float[] fftBuffer){
for (int i = 0; i < grid.getSize()[dim]; i++) {
grid.setValue(new Complex(fftBuffer[i*2],fftBuffer[i*2+1]),idx);
idx[dim]++;
}
}
private float[] dataToPlane(ComplexGrid grid, int startPos, int dim){
float[] out = null;
switch(dim){
case 1: //over x/z planes
out = new float[grid.getSize()[0]*grid.getSize()[2]*2];
for (int k = 0; k < grid.getSize()[2]; k++) {
System.arraycopy(grid.getAslinearMemory(), (k*grid.getSize()[0]*grid.getSize()[1]+startPos)*2, out, k*grid.getSize()[0]*2, grid.getSize()[0]*2);
}
break;
case 2: //over y/z planes
out = new float[grid.getSize()[1]*grid.getSize()[2]*2];
for (int k = 0; k < grid.getSize()[2]; k++) {
for (int j = 0; j < grid.getSize()[1]; j++) {
out[(k*grid.getSize()[1]+j)*2] = grid.getAslinearMemory()[2*(startPos+j*grid.getSize()[0]+k*grid.getSize()[0]*grid.getSize()[1])];
out[(k*grid.getSize()[1]+j)*2+1] = grid.getAslinearMemory()[2*(startPos+j*grid.getSize()[0]+k*grid.getSize()[0]*grid.getSize()[1])+1];
}
}
break;
default:
break;
}
return out;
}
private void planeToData(ComplexGrid grid, int startPos, int dim, float[] cplx){
switch(dim){
case 1: //over x/z planes
for (int k = 0; k < grid.getSize()[2]; k++) {
//System.arraycopy(grid.getAslinearMemory(), (k*grid.getSize()[0]*grid.getSize()[1]+startPos)*2, out, k*grid.getSize()[0], grid.getSize()[0]);
System.arraycopy(cplx, k*grid.getSize()[0]*2, grid.getAslinearMemory(), (k*grid.getSize()[0]*grid.getSize()[1]+startPos)*2, grid.getSize()[0]*2);
}
grid.notifyAfterWrite();
break;
case 2: //over y/z planes
for (int k = 0; k < grid.getSize()[2]; k++) {
for (int j = 0; j < grid.getSize()[1]; j++) {
grid.getAslinearMemory()[2*(startPos+j*grid.getSize()[0]+k*grid.getSize()[0]*grid.getSize()[1])]=cplx[(k*grid.getSize()[1]+j)*2];
grid.getAslinearMemory()[2*(startPos+j*grid.getSize()[0]+k*grid.getSize()[0]*grid.getSize()[1])+1]=cplx[(k*grid.getSize()[1]+j)*2+1];
}
}
grid.notifyAfterWrite();
break;
default:
break;
}
}
public static void main(String[] args) {
new ImageJ();
Grid3D blubb = ImageUtil.wrapImagePlus(IJ.openImage("D:\\Data\\ClackdoylePhantom\\Clackdoyle3D.tif"));
/*Grid3D blubb = new Grid3D(4,4,4);
DoubleFunction fct = (x -> Math.abs((double)x-1.5));
for (int i = 0; i < blubb.getSize()[0]; i++) {
for (int j = 0; j < blubb.getSize()[1]; j++) {
for (int k = 0; k < blubb.getSize()[2]; k++) {
if (fct.f(i) < 1 && fct.f(j) < 1 && fct.f(k) < 1)
blubb.setAtIndex(i, j, k, 1);
}
}
}*/
Fourier ft = new Fourier();
ComplexGrid cg3d = new ComplexGrid3D(blubb);
ComplexGrid cg3dRef = (ComplexGrid3D)cg3d.clone();
ComplexGrid cg2d = ((ComplexGrid3D)cg3d).getSubGrid(300).clone();
ComplexGrid cg2dRef = (ComplexGrid2D)cg2d.clone();
ComplexGrid cg1d = ((ComplexGrid2D)cg2d).getSubGrid(300).clone();
ComplexGrid cg1dRef = (ComplexGrid1D)cg1d.clone();
// 1D FFT testing
ft.fft(cg1d);ft.ifft(cg1d);
ft.fft(cg2d,0);ft.ifft(cg2d,0);ft.fft(cg2d,1);ft.ifft(cg2d,1);
ft.fft(cg3d,0);ft.ifft(cg3d,0);ft.fft(cg3d,1);ft.ifft(cg3d,1);ft.fft(cg3d,2);ft.ifft(cg3d,2);
ft.fft2(cg2d);ft.ifft2(cg2d);
ft.fft2(cg3d,0);ft.ifft2(cg3d,0);ft.fft2(cg3d,1);ft.ifft2(cg3d,1);ft.fft2(cg3d,2);ft.ifft2(cg3d,2);
ft.fft3(cg3d);
ft.ifft3(cg3d);
ComplexGridOperator cgo = new ComplexGridOperator();
cgo.subtractBy(cg1dRef, cg1d);
cg1dRef.show("1D-Diff");
cgo.subtractBy(cg2dRef, cg2d);
cg2dRef.show("2D-Diff");
cgo.subtractBy(cg3dRef, cg3d);
cg3dRef.show("3D-Diff");
}
}