/*
* Copyright (C) 2014 Marcel Pohlmann
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
package edu.stanford.rsl.tutorial.interpolation;
import java.util.ArrayList;
import ij.ImageJ;
import edu.stanford.rsl.conrad.data.generic.datatypes.Complex;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.Grid2DComplex;
import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
import edu.stanford.rsl.tutorial.fan.FanBeamBackprojector2D;
import edu.stanford.rsl.tutorial.fan.FanBeamProjector2D;
import edu.stanford.rsl.tutorial.fan.redundancy.BinaryWeights;
import edu.stanford.rsl.tutorial.fan.redundancy.CompensationWeights;
import edu.stanford.rsl.tutorial.fan.redundancy.ParkerWeights;
import edu.stanford.rsl.tutorial.fan.redundancy.SilverWeights;
import edu.stanford.rsl.tutorial.filters.RamLakKernel;
import edu.stanford.rsl.tutorial.filters.SheppLoganKernel;
import edu.stanford.rsl.tutorial.fourierConsistency.wedgefilter.DoubleWedgeFilterFanES;
/**
* This class implements different methods that can be used to estimate
* missing projections of a CT scan. The input parameter g represents the
* under-sampled sinogram, that already holds missing projections, w
* is the missing projection image, where 0:=missing projection
* and 1:=measured projection.
*
* @author Marcel Pohlmann
*
*/
public abstract class Interpolation {
/**
* Method that in-paints missing projections into the sparse sinogram (parallel-beam geometry), by
* iteratively reprojecting done in the projection space as described in the Paper:
* "Projection Space Iteration Reconstruction-Reprojection",
* J. H. Kim, K. Y. KWAK, IEEE Transactions on Medical Imaging, Vol. MI-4, No.3, September 1985
* This method assumes a parallel beam sinogram as input parameter!
*
* @param g Sparse sinogram
* @param w Missing projection mask
* @param params Fan-Beam parameters
* @param maxIter Number of maximum iterations
* @return In-painted sinogram
*/
public static Grid2D PSIRR(Grid2D g, Grid2D w, double rp, int maxIter, FanParameters params){
// TODO: CODE STILL BUGGY!
boolean debug = false;
if(debug){
new ImageJ();
}
double maxT = params.getMaxT(),
deltaT = params.getDeltaT();
double delta_l = g.getSpacing()[0],
delta_theta = g.getSpacing()[1];
Grid2D map = new Grid2D(g.getWidth(), g.getHeight());
Grid2D g_updated = (Grid2D) g.clone();
for(int i = 0; i < maxIter; ++i){
Grid2D g_filt = (Grid2D) g_updated.clone();
RamLakKernel ramLak = new RamLakKernel((int) (maxT / deltaT), deltaT);
for (int theta = 0; theta < g_filt.getSize()[1]; ++theta) {
ramLak.applyToGrid(g_filt.getSubGrid(theta));
}
for(int m = 0; m < w.getHeight(); ++m){
if(w.getAtIndex(0, m) == 0.0){
// found missing projection
for(int n = 0; n < g.getWidth(); ++n){
double theta_b = Math.acos(((n*delta_l - (maxT/2))/rp));
double newVal = 2*rp*Math.sin(theta_b)*g_filt.getAtIndex(n, m)*delta_l*delta_theta;
double sum = 0.0;
for(int j = 0; j < g.getHeight(); ++j){
int n1 = (int) Math.round(((rp*Math.cos(m*delta_theta - j*delta_theta + theta_b) + (maxT/2))/delta_l)),
n2 = (int) Math.round(((rp*Math.cos(m*delta_theta - j*delta_theta - theta_b) + (maxT/2))/delta_l));
if(debug){
map.setAtIndex(n1, j, 1.0f);
map.setAtIndex(n2, j, 2.0f);
}
for(int k = n1; k < n2; ++k){
if(j !=k){
sum += (1/(Math.sin(m*delta_theta - j*delta_theta)))*g_filt.getAtIndex(k, j);
}
}
}
newVal += sum*delta_theta*delta_l;
g_updated.setAtIndex(n, m, (float) newVal);
if(debug && (n%100 == 0)){
map.clone().show();
}
map = new Grid2D(g.getWidth(), g.getHeight());
}
//ramLak.applyToGrid(g_updated.getSubGrid(m));
}
}
for(int p = 0; p < g_updated.getHeight(); p++) {
for(int q = 0; q < g_updated.getWidth(); q++) {
if(Float.isNaN(g_updated.getAtIndex(q, p))){
g_updated.setAtIndex(q, p, 0.0f);
}
}
}
}
return g_updated;
}
/**
* Method that in-paints missing projections into the sparse sinogram, by iteratively reprojecting as
* described in the Paper:
* "Iterative Reconstruction-Reprojection: An Algorithm for Limited Data Cardiac-Computed Tomography",
* M. Nassi, IEEE Transactions on Biomedical Engineering, Vol. BME-29, No.5, May 1982
* @param g Sparse sinogram
* @param w Missing projection mask
* @param params Fan-Beam parameters
* @param maxIter Number of maximum iterations
* @param gpuDriven Enable OpenCL projector and backprojector
* @param imSize Image size of the reconstruction
* @param imSpacing Image spacing of the reconstruction
* @return In-painted sinogram
*/
public static Grid2D IRR(Grid2D g, Grid2D w, FanParameters params, int maxIter, boolean gpuDriven, int[] imSize, float[] imSpacing) {
boolean debug = false;
// extract FanParameters container
double maxT = params.getMaxT(),
deltaT = params.getDeltaT(),
focalLength = params.getFocalLength(),
maxBeta = params.getMaxBeta(),
deltaBeta = params.getDeltaBeta();
// extract image size and spacing for the reconstruction
int imgSzXMM = imSize[0], // [mm]
imgSzYMM = imSize[1]; // [mm]
float pxSzXMM = imSpacing[0],
pxSzYMM = imSpacing[1];
// make instance of backprojector
FanBeamBackprojector2D fbp = new FanBeamBackprojector2D(focalLength,
deltaT, deltaBeta, imgSzXMM, imgSzYMM);
// make instance of projector
FanBeamProjector2D fanBeamProjector = new FanBeamProjector2D(
focalLength, maxBeta, deltaBeta, maxT, deltaT);
double[] spacing = {deltaT, deltaBeta};
g.setSpacing(spacing);
int weightType = 0;
Grid2D RedundancyWeights;
switch (weightType) {
case 0:
RedundancyWeights = new ParkerWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
break;
case 1:
RedundancyWeights = new SilverWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
break;
case 2:
RedundancyWeights = new CompensationWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
break;
case 3:
RedundancyWeights = new BinaryWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
break;
default:
RedundancyWeights = new CompensationWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
}
Grid2D reference = (Grid2D) g.clone();
for(int iter = 0; iter < maxIter; ++iter){
// apply redundancy weights
Grid2D gFilt = (Grid2D) NumericPointwiseOperators.multipliedBy(g, RedundancyWeights);
gFilt.setSpacing(spacing);
// apply SheppLogan kernel
SheppLoganKernel sheppLogan = new SheppLoganKernel((int) (maxT / deltaT), deltaT);
for (int j=0; j < gFilt.getHeight(); j++){
sheppLogan.applyToGrid(gFilt.getSubGrid(j));
}
// do the backprojection
Grid2D reco = fbp.backprojectPixelDriven(gFilt);
// set negative values caused by the backprojection to zero
NumericPointwiseOperators.removeNegative(reco);
// set origin and spacing
reco.setOrigin(-imgSzXMM/2, -imgSzYMM/2);
reco.setSpacing(pxSzXMM, pxSzYMM);
Grid2D gHat;
// do the forward projection
if(gpuDriven == true){
gHat = fanBeamProjector.projectRayDrivenCL(reco);
}else{
gHat = fanBeamProjector.projectRayDriven(reco);
}
// insert current estimations of missing projections
insertionHelper(g, gHat, w);
if(debug && iter == (maxIter - 1)){
reference.clone().show("reference sinogram - should look like sparse sino");
g.clone().show("sparse sinogram inpainted with proj estimates");
reco.clone().show("latest reconstruction");
gHat.clone().show("FP of latest reconstruction");
Grid2D g_test = (Grid2D) g.clone();
NumericPointwiseOperators.multiplyBy(g_test, RedundancyWeights);
for (int j=0; j < gFilt.getHeight(); j++){
sheppLogan.applyToGrid(g_test.getSubGrid(j));
}
Grid2D test = fbp.backprojectPixelDrivenCL(g_test);
NumericPointwiseOperators.removeNegative(test);
test.clone().show("DEBUGGING");
}
}
return g;
}
/**
* Method that estimates missing projections by iteratively applying a wedge-filter in the Fourier Space.
*
* @see class DoubleWedgeFilterFanES
* @param g Sparse sinogram
* @param w Missing projection mask (0 denotes a missing projection)
* @param maxIter Number of maximum iterations
* @param zeroPad Apply zero-padding before FFT
* @param fanParams Array containing maximum object radius rp, source-isocenter distance L, detector-isocenter distance D
* @return sinogram holding missing projections
*/
public static Grid2D iterativeWedgeFilter(Grid2D g, Grid2D w, int maxIter, boolean zeroPad, double[] fanParams) {
// extract fan beam parameters
double rp = fanParams[0];
double _L = fanParams[1];
double _D = fanParams[2];
// initialize Complex Grid2D that will hold the Fourier transform of the sinogram
Grid2DComplex c = new Grid2DComplex(g, zeroPad);
// compute and set frequency axis
int[] size = c.getSize();
double[] spatialSpacing = g.getSpacing();
double[] fMax = {0.5*(1/spatialSpacing[0]), 0.5*(1/spatialSpacing[1])};
double[] fSpacing = {2*fMax[0]/size[0], 2*fMax[1]/size[1]};
double[] fOrigin = {-fMax[0]+fSpacing[0], -fMax[1]+fSpacing[1]};
// compute the double wedge filter
DoubleWedgeFilterFanES wedgeFilter = new DoubleWedgeFilterFanES(c.getSize(), fSpacing, fOrigin, rp, _L, _D);
// erode wedge filter to avoid discretization problems
wedgeFilter.erode(7);
for(int iter = 0; iter < maxIter; ++iter) {
// compute 2D-Fourier transform of the current estimate of the sinogram
Grid2DComplex cCopy = c.clone();
cCopy.transformForward();
// FFT-shift the 2D-FT to fit the double wedge filter axis', where the zero frequency is centered
cCopy.fftshift();
// apply wedge filter on 2D-FT
for(int i = 0; i < wedgeFilter.getHeight(); ++i) {
for(int j = 0; j < wedgeFilter.getWidth(); ++j) {
if(wedgeFilter.getPixelValue(j, i) > 0.0) {
}else {
cCopy.putPixelValue(2*j, i, 0.0);
cCopy.putPixelValue(2*j+1, i, 0.0);
}
}
}
// invert the FFT-shift
cCopy.ifftshift();
// apply the inverse 2D-FT
cCopy.transformInverse();
Grid2D magnitude = cCopy.getMagnSubGrid(0, 0, g.getWidth(), g.getHeight());
// c2 now holds estimated projections;
// note that the measured projections are also affected by the double wedge, therefore ..
Grid2DComplex c2 = new Grid2DComplex(magnitude, zeroPad);
// .. insert ONLY the new projection estimates
insertionHelper(c, c2, w);
}
return c.getMagnSubGrid(0, 0, g.getWidth(), g.getHeight());
}
/**
* Method that estimates missing projections by iteratively decreasing the current rp and applying a wedge-filter
* in the Fourier Space
*
* @see class DoubleWedgeFilterFanES
* @param g Sparse sinogram
* @param w Missing projection mask (0 denotes a missing projection)
* @param maxIter Number of maximum iterations
* @param zeroPad Apply zero-padding before FFT
* @param fanParams Array containing maximum object radius rp, source-isocenter distance L, detector-isocenter distance D
* @return In-painted sinogram
*/
public static Grid2D iterativeWedgeFilter_opti(Grid2D g, Grid2D w, int maxIter, boolean zeroPad, double[] fanParams) {
// extract Fan Beam Parameters
double rp = fanParams[0];
double _L = fanParams[1];
double _D = fanParams[2];
// initialize output Grid2D as the same size as g
Grid2D output = new Grid2D(g.getWidth(), g.getHeight());
// define step size for increasing rp
double delta_rp = rp/2;
// iterate over decreasing rp
for(double cur_rp = rp; cur_rp > 0.0; cur_rp -= delta_rp){
Grid2D g_original = (Grid2D) g.clone();
// initialize Complex Grid2D that will hold the Fourier transform of the sinogram
Grid2DComplex c = new Grid2DComplex(g_original, zeroPad);
// compute and set frequency axis
int[] size = c.getSize();
double[] spatialSpacing = g.getSpacing();
double[] fMax = {0.5*(1/spatialSpacing[0]), 0.5*(1/spatialSpacing[1])};
double[] fSpacing = {2*fMax[0]/size[0], 2*fMax[1]/size[1]};
double[] fOrigin = {-fMax[0]+fSpacing[0], -fMax[1]+fSpacing[1]};
// compute the double wedge filter
DoubleWedgeFilterFanES wedgeFilter = new DoubleWedgeFilterFanES(c.getSize(), fSpacing, fOrigin, cur_rp, _L, _D);
// erode wedge filter to avoid errors caused by discretization
wedgeFilter.erode(7);
for(int iter = 0; iter < maxIter; ++iter) {
// compute 2D-Fourier transform of the current estimate of the sinogram
Grid2DComplex cCopy = c.clone();
cCopy.transformForward();
// FFT-shift the 2D-FT to fit the double wedge filter axis', where the zero frequency is centered
cCopy.fftshift();
// apply wedge filter on 2D-FT
for(int i = 0; i < wedgeFilter.getHeight(); ++i) {
for(int j = 0; j < wedgeFilter.getWidth(); ++j) {
if(wedgeFilter.getPixelValue(j, i) > 0.0) {
}else {
cCopy.putPixelValue(2*j, i, 0.0);
cCopy.putPixelValue(2*j+1, i, 0.0);
}
}
}
// invert the FFT-shift
cCopy.ifftshift();
// apply the inverse 2D-FT
cCopy.transformInverse();
Grid2D magnitude = cCopy.getMagnSubGrid(0, 0, g.getWidth(), g.getHeight());
// c2 now holds estimated projections;
// note that the measured projections are also affected by the double wedge, therefore ..
Grid2DComplex c2 = new Grid2DComplex(magnitude, zeroPad);
// .. insert ONLY the new projection estimates
insertionHelper(c, c2, w);
}
Grid2D c_magnitude = c.getMagnSubGrid(0, 0, g.getWidth(), g.getHeight());
// fill estimated projections into corresponding positions in the sinogram
for(int theta = 0; theta < output.getHeight(); ++theta){
for(int s = (int) ((output.getWidth()/2) - cur_rp); s <= ((output.getWidth()/2) + cur_rp); s++){
output.setAtIndex(s, theta, c_magnitude.getAtIndex(s, theta));
}
}
}
return output;
}
/**
* Method that fixes image defects by patchwise Spectral Deconvolution as decribed in the paper:
* "Defect Interpolation in Digital Radiography - How Object-Oriented Transform Coding Helps",
* T. Aach, V. Metzler, SPIE Vol. 4322: Medical Imaging, February 2001
*
* @see method Grid2D SpectralDeconvoltion(..)
* @param g_original Image with defect
* @param w Defect mask (0 denotes the defect)
* @param patchSize Holds the desired size of the patch
* @param maxIter Number of maximum iterations
* @param zeroPadSignal Apply zero-padding on image before FFT
* @return In-painted image
*/
public static Grid2D ApplyPatchwiseSpectralDeconvolution(Grid2D g_original, Grid2D w, int[] patchSize, int maxIter,
boolean zeroPadSignal){
Grid2D g = (Grid2D) g_original.clone();
// move patch over image with an overlay of half the patch size (in both directions)
for(int i = 0; i < g.getHeight() - (patchSize[1]/2); i += (patchSize[1]/2)){
for(int j = 0; j < g.getWidth() - (patchSize[0]/2); j += (patchSize[0]/2)){
// patch will hold the current image sector
Grid2D patch = new Grid2D(patchSize[0], patchSize[1]);
// patch mask will hold the corresponding sector of the defect mask
Grid2D patch_mask = new Grid2D(patchSize[0], patchSize[1]);
// fill patch and patch_mask with values
for(int l = 0; l < patchSize[0]; ++l){
for(int k = 0; k < patchSize[1]; ++k){
patch.setAtIndex(l, k, g_original.getAtIndex(l + j, k + i));
patch_mask.setAtIndex(l, k, w.getAtIndex(l + j, k + i));
}
}
// apply Spectral Deconvolution on the current patch
Grid2D result = SpectralDeconvoltion((Grid2D) patch, (Grid2D) patch_mask, maxIter, zeroPadSignal);
// write the result back into g
for(int l = 0; l < patchSize[0]; ++l){
for(int k = 0; k < patchSize[1]; ++k){
g.setAtIndex(l + j, k + i, result.getAtIndex(l, k));
}
}
}
}
return g;
}
/**
* Method that fixes image defects by Spectral Deconvolution as decribed in the paper:
* "Defect Interpolation in Digital Radiography - How Object-Oriented Transform Coding Helps",
* T. Aach, V. Metzler, SPIE Vol. 4322: Medical Imaging, February 2001
*
* @param g_original Image with defect
* @param w Defect mask (0 denotes the defect)
* @param maxIter Number of maximum iterations
* @param zeroPadSignal Apply zero-padding on image before FFT
* @return In-painted image
*/
public static Grid2D SpectralDeconvoltion(Grid2D g_original, Grid2D w, int maxIter, boolean zeroPadSignal) {
Grid2D g = new Grid2D(g_original);
double maxDeltaE_G_Ratio = Double.POSITIVE_INFINITY;
double maxDeltaE_G_Ratio_Tres = 1.0e-6;
Grid2DComplex G = new Grid2DComplex(g, zeroPadSignal);
G.transformForward();
Grid2DComplex W = new Grid2DComplex(w, zeroPadSignal);
W.transformForward();
int[] dim = G.getSize();
int[] halfDim = {dim[0]/2, dim[1]/2};
// Initialization
Grid2DComplex FHat = new Grid2DComplex(dim[0], dim[1], false);
Grid2DComplex FHatNext = new Grid2DComplex(dim[0], dim[1], false);
double lastEredVal = 0;
for(int i = 0; i < maxIter; i++) {
// Check convergence criterion
if(maxDeltaE_G_Ratio <= maxDeltaE_G_Ratio_Tres) {
System.out.println("Break in " + i + "th iteration");
break;
}
// In the i-th iteration select the line pair s1,t1
// which maximizes the energy reduction [Paragraph after Eq. (16) in the paper]
double maxDeltaE_G = Double.NEGATIVE_INFINITY;
ArrayList<Integer[]> sj1 = null;
// Find the entries with the maximum values and write them in sj1
for(int j = 0; j < dim[0]; j++) {
for(int k = 0; k < dim[1]; k++) {
double val = G.getAtIndex(j, k);
if(val > maxDeltaE_G) {
sj1 = new ArrayList<Integer[]>();
sj1.add(new Integer[]{j,k});
maxDeltaE_G = val;
}else if(val == maxDeltaE_G) {
sj1.add(new Integer[]{j, k});
}
}
}
int idx = (int) Math.floor(Math.random() * sj1.size());
int s1 = sj1.get(idx)[0];
int t1 = sj1.get(idx)[1];
// Calculate the ratio of energy reduction in comparison to the last iteration
if(i > 0) {
maxDeltaE_G_Ratio = Math.abs((maxDeltaE_G - lastEredVal) / maxDeltaE_G);
}
// Save the last energy reduction value for next iteration
lastEredVal = maxDeltaE_G;
// Compute the corresponding linepair s2, t2:
// mirror the positions at halfDim
int s2 = (s1 > 0) ? dim[0] - (s1 % dim[0]) : s1;
int t2 = (t1 > 0) ? dim[1] - (t1 % dim[1]) : t1;
// This we require in the next step:
// [Paragraph after Eq. (17) in the paper]
int twice_s1 = (2 * s1) % dim[0];
int twice_t1 = (2 * t1) % dim[1];
boolean specialCase = false;
if( (s1 == 0 && t1 == 0) ||
(s1 == 0 && t1 == halfDim[1])||
(s1 == halfDim[0] && t1 == 0) ||
(s1 == halfDim[0] && t1 == halfDim[1])) {
//SPECIAL CASES
specialCase = true;
Complex FHatNextVal = new Complex(FHatNext.getRealAtIndex(s1, t1), FHatNext.getImagAtIndex(s1, t1));
Complex GVal = new Complex(G.getRealAtIndex(s1, t1), G.getImagAtIndex(s1, t1));
Complex WVal = new Complex(W.getRealAtIndex(0, 0), W.getImagAtIndex(0, 0));
Complex res = FHatNextVal.add(GVal.mul(dim[0] * dim[1]).div(WVal));
FHatNext.setRealAtIndex(s1, t1, (float) res.getReal());
FHatNext.setImagAtIndex(s1, t1, (float) res.getImag());
}else {
// GENERAL CASE
Complex FHatNextVal_s1t1 = new Complex(FHatNext.getRealAtIndex(s1, t1), FHatNext.getImagAtIndex(s1, t1));
Complex FHatNextVal_s2t2 = new Complex(FHatNext.getRealAtIndex(s2, t2), FHatNext.getImagAtIndex(s2, t2));
Complex GVal = new Complex(G.getRealAtIndex(s1, t1), G.getImagAtIndex(s1, t1));
Complex WVal00 = new Complex(W.getRealAtIndex(0, 0), W.getImagAtIndex(0, 0));
Complex WValTwice = new Complex(W.getRealAtIndex(twice_s1, twice_t1), W.getImagAtIndex(twice_s1, twice_t1));
Complex tVal = ((GVal.mul(WVal00)).sub((GVal.getConjugate().mul(WValTwice)))).mul((dim[0] * dim[1]));
tVal = tVal.div(WVal00.getMagn() * WVal00.getMagn() - WValTwice.getMagn() * WValTwice.getMagn());
Complex res1 = FHatNextVal_s1t1.add(tVal);
Complex res2 = FHatNextVal_s2t2.add(tVal.getConjugate());
FHatNext.setRealAtIndex(s1, t1, (float) res1.getReal());
FHatNext.setImagAtIndex(s1, t1, (float) res1.getImag());
FHatNext.setRealAtIndex(s2, t2, (float) res2.getReal());
FHatNext.setImagAtIndex(s2, t2, (float) res2.getImag());
}
// Form the new error spectrum
updateSpectrum(G, FHatNext, FHat, W, s1, t1, specialCase);
G.setAtIndex(s1, t1, 0);
if(!specialCase) {
G.setAtIndex(s2, t2, 0);
}
FHat = new Grid2DComplex(FHatNext);
}
// Compute the inverse Fourier Transform of the estimated image
FHat.transformInverse();
// Insert missing pixels
for(int j = 0; j < g.getSize()[1]; j++) {
for(int i = 0; i < g.getSize()[0]; i++) {
if(w.getAtIndex(i, j) == 0) {
g.setAtIndex(i, j, FHat.getRealAtIndex(i, j));
}
}
}
return g;
}
private static void insertionHelper(Grid2D c, Grid2D cHat , Grid2D w) {
// replace missing lines in c with samples of cHat
for(int i = 0; i < w.getHeight(); ++i) {
for(int j = 0; j < w.getWidth(); ++j){
if(w.getAtIndex(j, i) == 0.0) {
c.setAtIndex(j, i, cHat.getAtIndex(j, i));
}
}
}
}
private static void updateSpectrum(Grid2DComplex G, Grid2DComplex FHatNext, Grid2DComplex FHat, Grid2DComplex W, int s1, int t1, boolean specialCase) {
int[] sz = FHatNext.getSize();
// Accumulation: Update pair (s1,t1),(s2,t2)
Complex Fst = new Complex(FHatNext.getRealAtIndex(s1, t1) - FHat.getRealAtIndex(s1, t1), FHatNext.getImagAtIndex(s1, t1) - FHat.getImagAtIndex(s1, t1));
Complex Fstc = Fst.getConjugate();
int divNr = sz[0] * sz[1];
// Compute the new error spectrum
for(int j = 0; j < sz[1]; j++) {
for(int i = 0; i < sz[0]; i++) {
if(specialCase) {
int xneg = (i - s1) % sz[0];
int yneg = (j - t1) % sz[1];
if(xneg < 0) {
xneg = sz[0] + xneg;
}
if(yneg < 0) {
yneg = sz[1] + yneg;
}
Complex GVal = new Complex(G.getRealAtIndex(i, j), G.getImagAtIndex(i, j));
G.setRealAtIndex(i, j, (float) GVal.getReal());
G.setImagAtIndex(i, j, (float) GVal.getImag());
}else {
int xpos = (i + s1) % sz[0];
int ypos = (j + t1) % sz[1];
int xneg = (i - s1) % sz[0];
int yneg = (j - t1) % sz[1];
if(xneg < 0) {
xneg = sz[0] + xneg;
}
if(yneg < 0) {
yneg = sz[1] + yneg;
}
Complex WPos = new Complex(W.getRealAtIndex(xpos, ypos), W.getImagAtIndex(xpos, ypos));
Complex WNeg = new Complex(W.getRealAtIndex(xneg, yneg), W.getImagAtIndex(xneg, yneg));
Complex GVal = new Complex(G.getRealAtIndex(i, j), G.getImagAtIndex(i, j));
GVal = GVal.sub(((Fst.mul(WNeg)).add(Fstc.mul(WPos))).div(divNr));
G.setRealAtIndex(i, j, (float) GVal.getReal());
G.setImagAtIndex(i, j, (float) GVal.getImag());
}
}
}
}
}