/*
* Copyright (C) 2014 - Martin Berger
* Copyright (C) 2017 - Christopher Syben
* 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.ImageJ;
import edu.stanford.rsl.conrad.data.generic.datatypes.Complex;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
public class ComplexGrid3D extends ComplexGrid {
/** Linear float array stores the complex numbers in an alternating manner [Re1, Imag1, Re2, Imag2, ..., ReN,ImagN]. */
private float[] buffer;
private ComplexGrid2D[] subGrids = null;
public ComplexGrid3D(int width, int height, int depth) {
assert width > 0 && height > 0 && depth > 0;
this.size = new int[] {width, height, depth};
this.spacing = new double[]{1,1,1};
this.origin = new double[3];
buffer = new float[width*height*depth*2];
initialize(width, height, depth);
notifyAfterWrite();
}
public ComplexGrid3D(ComplexGrid3D input){
this.size = new int[] {input.size[0], input.size[1], input.size[2]};
this.spacing = new double[] {input.spacing[0],input.spacing[1], input.spacing[2]};
this.origin = new double[] {input.origin[0],input.origin[1], input.origin[2]};
buffer = new float[input.getNumberOfElements()*2];
System.arraycopy(input.getBuffer(), input.getOffset()*2, buffer, 0, buffer.length);
initialize(input.size[0], input.size[1], input.size[2]);
if(input.openCLactive)
activateCL();
}
public ComplexGrid3D(Grid3D input){
this.spacing = new double[] {input.getSpacing()[0],input.getSpacing()[1], input.getSpacing()[2]};
this.origin = new double[] {input.getOrigin()[0],input.getOrigin()[1], input.getOrigin()[2]};
this.size = new int[] {input.getSize()[0], input.getSize()[1], input.getSize()[2]};
buffer = new float[input.getNumberOfElements()*2];
initialize(input.getSize()[0], input.getSize()[1], input.getSize()[2]);
for (int k = 0; k < input.getSize()[2]; k++) {
for (int j = 0; j < input.getSize()[1]; j++) {
for (int i = 0; i < input.getSize()[0]; i++) {
this.setAtIndex(i, j, k, input.getAtIndex(i, j, k));
}
}
}
}
public ComplexGrid3D(Grid3D realInput, Grid3D imagInput){
this.spacing = new double[] {realInput.getSpacing()[0],realInput.getSpacing()[1], realInput.getSpacing()[2]};
this.origin = new double[] {realInput.getOrigin()[0],realInput.getOrigin()[1], realInput.getOrigin()[2]};
this.size = new int[] {realInput.getSize()[0], realInput.getSize()[1], realInput.getSize()[2]};
buffer = new float[realInput.getNumberOfElements()*2];
initialize(realInput.getSize()[0], realInput.getSize()[1], realInput.getSize()[2]);
for(int k = 0; k < realInput.getSize()[2]; k++){
for(int j = 0; j < realInput.getSize()[1]; j++){
for(int i = 0; i < realInput.getSize()[0]; i++ ){
this.setAtIndex(i, j, k, realInput.getAtIndex(i, j, k), imagInput.getAtIndex(i, j, k));
}
}
}
}
protected void initialize(int width, int height, int depth)
{
subGrids = new ComplexGrid2D[depth];
for (int i = 0; i < subGrids.length; i++) {
subGrids[i] = new ComplexGrid2D(buffer, i*width*height, width, height);
}
}
public void fftshift()
{
ComplexGrid2D[] unshifted = subGrids.clone();
for(int i = 0; i < subGrids.length;i++ ){
unshifted[i].fftshift();
}
int depth = this.getSize()[2];
int pD = (int)Math.ceil(depth/2.0);
for(int k = pD; k < this.getSize()[2];k++)
{
this.setSubGrid(k-pD,unshifted[k]);
}
for(int k = 0; k < pD;k++)
{
this.setSubGrid(k+depth-pD,unshifted[k]);
}
}
public float[] getBuffer(){
notifyBeforeRead();
return buffer;
}
public void multiplyAtIndex(int i, int j, int k, float val) {
multiplyAtIndex(i, j, k, val, 0);
}
public void multiplyAtIndex(int i, int j, int k, float real, float imag) {
Complex in = new Complex(real, imag);
multiplyAtIndex(i, j, k, in);
}
public void multiplyAtIndex(int i, int j, int k, Complex in) {
setAtIndex(i, j, k, getAtIndex(i, j, k).mul(in));
}
public void divideAtIndex(int i, int j, int k, float val) {
divideAtIndex(i, j, k, val, 0);
}
public void divideAtIndex(int i, int j, int k, float real, float imag) {
Complex in = new Complex(real, imag);
divideAtIndex(i, j, k, in);
}
public void divideAtIndex(int i, int j, int k, Complex in) {
setAtIndex(i, j, k, getAtIndex(i, j, k).div(in));
}
public void addAtIndex(int i, int j, int k, float val) {
addAtIndex(i, j, k, val, 0);
}
public void addAtIndex(int i, int j, int k, float real, float imag) {
Complex in = new Complex(real, imag);
addAtIndex(i, j, k, in);
}
public void addAtIndex(int i, int j, int k, Complex in) {
setAtIndex(i, j, k, getAtIndex(i, j, k).add(in));
}
public void subtractAtIndex(int i, int j, int k, float val) {
subtractAtIndex(i, j, k, val, 0);
}
public void subtractAtIndex(int i, int j, int k, float real, float imag) {
Complex in = new Complex(real, imag);
subtractAtIndex(i, j, k, in);
}
public void subtractAtIndex(int i, int j, int k, Complex in) {
setAtIndex(i, j, k, getAtIndex(i, j, k).sub(in));
}
public void setAtIndex(int i, int j, int k, float val){
setAtIndex(i, j, k, val, 0);
}
public void setAtIndex(int i, int j, int k, float real, float imag){
Complex in = new Complex(real, imag);
setAtIndex(i, j, k, in);
}
public void setAtIndex(int i, int j, int k, Complex val){
this.getSubGrid(k).getSubGrid(j).setAtIndex(i,(float)val.getReal(),(float)val.getImag());
notifyAfterWrite();
}
public Complex getAtIndex(int i, int j, int k){
notifyBeforeRead();
return getSubGrid(k).getSubGrid(j).getAtIndex(i);
}
public ComplexGrid2D getSubGrid(int k){
notifyBeforeRead();
return subGrids[k];
}
public void setSubGrid(int k, ComplexGrid2D slice){
subGrids[k]=slice;
notifyAfterWrite();
}
@Override
public Complex getValue(int... idx) {
return getAtIndex(idx[0],idx[1],idx[2]);
}
@Override
public void setValue(Complex val, int... idx) {
setAtIndex(idx[0], idx[1], idx[2], val);
}
@Override
public ComplexGrid clone() {
notifyBeforeRead();
return new ComplexGrid3D(this);
}
@Override
public String toString() {
String result = new String();
result += "[";
for (int j = 0; j < size[1]; ++j) {
if (j != 0) result += "; ";
result += getSubGrid(j).toString();
}
result += "]";
return result;
}
@Override
public Grid3D getRealGrid() {
return getRealSubGrid(0, 0, 0, size[0], size[1], size[2]);
}
@Override
public Grid3D getImagGrid() {
return getImagSubGrid(0, 0, 0, size[0], size[1], size[2]);
}
@Override
public Grid3D getMagGrid() {
return getMagSubGrid(0, 0, 0, size[0], size[1], size[2]);
}
@Override
public Grid3D getPhaseGrid() {
return getPhaseSubGrid(0, 0, 0, size[0], size[1], size[2]);
}
public Grid3D getRealSubGrid(final int startX, final int startY, final int startZ, final int lengthX, final int lengthY, final int lengthZ){
Grid3D subgrid = new Grid3D(lengthX,lengthY,lengthZ);
for (int k=0; k < lengthZ; ++k){
for (int j=0; j < lengthY; ++j){
for (int i=0; i < lengthX; ++i){
subgrid.setAtIndex(i, j, k, this.getRealAtIndex(i+startX,j+startY,k+startZ));
}
}
}
return subgrid;
}
public Grid3D getImagSubGrid(final int startX, final int startY, final int startZ, final int lengthX, final int lengthY, final int lengthZ){
Grid3D subgrid = new Grid3D(lengthX,lengthY,lengthZ);
for (int k=0; k < lengthZ; ++k){
for (int j=0; j < lengthY; ++j){
for (int i=0; i < lengthX; ++i){
subgrid.setAtIndex(i, j, k, this.getImagAtIndex(i+startX,j+startY,k+startZ));
}
}
}
return subgrid;
}
public Grid3D getMagSubGrid(final int startX, final int startY, final int startZ, final int lengthX, final int lengthY, final int lengthZ){
Grid3D subgrid = new Grid3D(lengthX,lengthY,lengthZ);
for (int k=0; k < lengthZ; ++k){
for (int j=0; j < lengthY; ++j){
for (int i=0; i < lengthX; ++i){
float x = this.getRealAtIndex(i+startX,j+startY,k+startZ);
float y = this.getImagAtIndex(i+startX,j+startY,k+startZ);
subgrid.setAtIndex(i,j,k, (float)Math.sqrt(x*x+y*y));
}
}
}
return subgrid;
}
public Grid3D getPhaseSubGrid(final int startX, final int startY, final int startZ, final int lengthX, final int lengthY, final int lengthZ){
Grid3D subgrid = new Grid3D(lengthX,lengthY,lengthZ);
for (int k=0; k < lengthZ; ++k){
for (int j=0; j < lengthY; ++j){
for (int i=0; i < lengthX; ++i){
float x = this.getRealAtIndex(i+startX,j+startY,k+startZ);
float y = this.getImagAtIndex(i+startX,j+startY,k+startZ);
subgrid.setAtIndex(i,j,k, (float)Math.atan2(y,x));
}
}
}
return subgrid;
}
@Override
public float[] getAslinearMemory() {
notifyBeforeRead();
return buffer;
}
@Override
public void setAslinearMemory(float[] buffer) {
this.buffer = buffer;
notifyAfterWrite();
}
@Override
public float getRealAtIndex(int... idx) {
return (float)getAtIndex(idx[0],idx[1],idx[2]).getReal();
}
@Override
public void setRealAtIndex(float val, int... idx) {
setAtIndex(idx[0],idx[1],idx[2], val);
}
@Override
public float getImagAtIndex(int... idx) {
return (float)getAtIndex(idx[0],idx[1],idx[2]).getImag();
}
@Override
public void setImagAtIndex(float val, int... idx) {
Complex cval = getAtIndex(idx[0],idx[1],idx[2]);
cval.setImag(val);
setAtIndex(idx[0],idx[1],idx[2], cval);
}
@Override
public int getOffset() {
return 0;
}
public static void main(String[] args) {
new ImageJ();
int pow2 = 226;
ComplexGrid3D bla = new ComplexGrid3D(pow2,pow2,pow2);
int max = bla.getNumberOfElements();
for (int k = 0; k < max; k++) {
float h = (float)((double)k*Math.sqrt(2)/bla.getNumberOfElements());
bla.getBuffer()[2*k]=h;
bla.getBuffer()[2*k+1]=-1*h;
}
ComplexGrid3D gc = null;
ComplexGrid3D gc1 = null;
ComplexGrid3D gc2 = null;
long[] ellapsed = new long[3];
gc = new ComplexGrid3D(bla);
gc1 = new ComplexGrid3D(gc);
gc2 = new ComplexGrid3D(gc);
bla.show("GT");
Fourier ft = new Fourier();
int reps = 1;
for (int i = 0; i < reps; i++) {
long start = System.nanoTime();
ft.fft(gc);
ft.ifft(gc);
ellapsed[0] += System.nanoTime() - start;
start = System.nanoTime();
ft.fft2(gc1);
ft.ifft2(gc1);
ellapsed[1] += System.nanoTime() - start;
start = System.nanoTime();
ft.fft3(gc2);
ft.ifft3(gc2);
ellapsed[2] += System.nanoTime() - start;
}
for (int i = 0; i < ellapsed.length; i++) {
System.out.println(((double)ellapsed[i])/1e6/(double)reps);
}
gc.show("FFT1");
gc1.show("FFT2");
gc2.show("FFT3");
}
}