package features;

import Utils.Utils;
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.measure.Calibration;
import ij.process.ByteProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import ij.process.ShortProcessor;
import math3d.Eigensystem2x2Double;
import math3d.Eigensystem2x2Float;
import math3d.Eigensystem3x3Double;
import math3d.Eigensystem3x3Float;
import math3d.JacobiDouble;
import math3d.JacobiFloat;

public class ComputeCurvatures implements Runnable {
   private boolean _3D;
   private ComputeCurvatures.FloatArray data;
   private double[][] hessianMatrix;
   private float[] eigenValues = new float[3];
   private ComputeCurvatures.FloatArray3D[] result3D;
   private ComputeCurvatures.FloatArray3D result2D;
   private double min = Double.MAX_VALUE;
   private double max = Double.MIN_VALUE;
   protected ImagePlus imp;
   protected double sigma;
   protected boolean useCalibration;
   protected GaussianGenerationCallback callback;
   private boolean cancelGeneration = false;

   public ComputeCurvatures() {
   }

   public ComputeCurvatures(ImagePlus imp, double sigma, GaussianGenerationCallback callback, boolean useCalibration) {
      this.imp = imp;
      this.sigma = sigma;
      this.callback = callback;
      this.useCalibration = useCalibration;
   }

   public void cancelGaussianGeneration() {
      this.cancelGeneration = true;
   }

   @Override
   public void run() {
      if (this.imp == null) {
         IJ.error("BUG: imp should not be null - are you using the right constructor?");
      } else {
         this.setup();
      }
   }

   public void setup() {
      try {
         if (this.imp == null) {
            IJ.error("BUG: imp should not be null - are you using the right constructor?");
         } else {
            if (this.callback != null) {
               this.callback.proportionDone(0.0);
            }

            if (this.imp.getStackSize() > 1) {
               ImageStack stack = this.imp.getStack();
               this._3D = true;
               this.data = this.StackToFloatArray(stack);
               if (this.data == null) {
                  return;
               }
            } else {
               this._3D = false;
               this.data = this.ImageToFloatArray(this.imp.getProcessor());
               if (this.data == null) {
                  return;
               }
            }

            boolean computeGauss = true;
            boolean showGauss = true;
            Calibration calibration = this.imp.getCalibration();
            if (!Utils.isReleaseVersion) {
               System.out.println("Computing Gauss image");
            }

            if (this._3D) {
               this.data = this.computeGaussianFastMirror(
                  (ComputeCurvatures.FloatArray3D)this.data, (float)this.sigma, this.callback, this.useCalibration ? calibration : null
               );
               if (this.data == null) {
                  if (this.callback != null) {
                     this.callback.proportionDone(-1.0);
                  }

                  return;
               }

               if (showGauss) {
                  this.FloatArrayToStack((ComputeCurvatures.FloatArray3D)this.data, "Gauss image", 0.0F, 255.0F).show();
               }
            } else {
               this.data = this.computeGaussianFastMirror(
                  (ComputeCurvatures.FloatArray2D)this.data, (float)this.sigma, this.callback, this.useCalibration ? calibration : null
               );
               if (this.data == null) {
                  if (this.callback != null) {
                     this.callback.proportionDone(-1.0);
                  }

                  return;
               }

               if (showGauss) {
                  FloatArrayToImagePlus((ComputeCurvatures.FloatArray2D)this.data, "Gauss image", 0.0F, 255.0F).show();
               }
            }

            if (this.callback != null) {
               this.callback.proportionDone(1.0);
            }
         }
      } catch (OutOfMemoryError var4) {
         long requiredMiB = (long)(this.imp.getWidth() * this.imp.getHeight() * this.imp.getStackSize() * 4 / 1048576);
         IJ.error("Out of memory when calculating the Gaussian convolution of the image (requires " + requiredMiB + "MiB");
         if (this.callback != null) {
            this.callback.proportionDone(-1.0);
         }
      }
   }

   public void setData(ImagePlus imp) {
      this.data = this.ImageToFloatArray(imp.getProcessor());
      if (this.data == null && !Utils.isReleaseVersion) {
         System.err.println("Data = null!!!!");
      }
   }

   public void setSigma(double sigma) {
      this.sigma = sigma;
   }

   public ComputeCurvatures.FloatArray getData() {
      return this.data;
   }

   public boolean hessianEigenvaluesAtPoint2D(
      int x, int y, boolean orderOnAbsoluteSize, float[] result, boolean normalize, boolean fixUp, float sepX, float sepY
   ) {
      long start = System.currentTimeMillis();
      if (this._3D) {
         if (!Utils.isReleaseVersion) {
            System.err.println("hessianEigenvaluesAtPoint2D( x, y, z, ... ) is only for 2D data.");
         }

         return false;
      } else {
         ComputeCurvatures.FloatArray2D data2D = (ComputeCurvatures.FloatArray2D)this.data;
         if (fixUp) {
            if (x == 0) {
               x = 1;
            }

            if (x == data2D.width - 1) {
               x = data2D.width - 2;
            }

            if (y == 0) {
               y = 1;
            }

            if (y == data2D.height - 1) {
               y = data2D.height - 2;
            }
         }

         float[][] hessianMatrix = this.computeHessianMatrix2DFloat(data2D, x, y, this.sigma, sepX, sepY);
         float[] eigenValues = this.computeEigenValues(hessianMatrix);
         if (eigenValues == null) {
            return false;
         } else {
            float e0 = eigenValues[0];
            float e1 = eigenValues[1];
            float e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
            float e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
            if (e0c <= e1c) {
               result[0] = e0;
               result[1] = e1;
            } else {
               result[0] = e1;
               result[1] = e0;
            }

            if (normalize) {
               float divideBy = Math.abs(result[1]);
               result[0] /= divideBy;
               result[1] /= divideBy;
            }

            return true;
         }
      }
   }

   public long[] hessianEigenvaluesAtPoint2D2(
      int x, int y, boolean orderOnAbsoluteSize, float[] result, boolean normalize, boolean fixUp, float sepX, float sepY
   ) {
      long start = System.currentTimeMillis();
      if (this._3D) {
         if (!Utils.isReleaseVersion) {
            System.err.println("hessianEigenvaluesAtPoint2D( x, y, z, ... ) is only for 2D data.");
         }

         return null;
      } else {
         ComputeCurvatures.FloatArray2D data2D = (ComputeCurvatures.FloatArray2D)this.data;
         if (fixUp) {
            if (x == 0) {
               x = 1;
            }

            if (x == data2D.width - 1) {
               x = data2D.width - 2;
            }

            if (y == 0) {
               y = 1;
            }

            if (y == data2D.height - 1) {
               y = data2D.height - 2;
            }
         }

         long hessianMatrixStart = System.currentTimeMillis();
         float[][] hessianMatrix = this.computeHessianMatrix2DFloat(data2D, x, y, this.sigma, sepX, sepY);
         long hessianMatrixEnd = System.currentTimeMillis();
         float[] eigenValues = this.computeEigenValues(hessianMatrix);
         long eigenValuesEnd = System.currentTimeMillis();
         if (eigenValues == null) {
            return null;
         } else {
            float e0 = eigenValues[0];
            float e1 = eigenValues[1];
            float e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
            float e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
            if (e0c <= e1c) {
               result[0] = e0;
               result[1] = e1;
            } else {
               result[0] = e1;
               result[1] = e0;
            }

            if (normalize) {
               float divideBy = Math.abs(result[1]);
               result[0] /= divideBy;
               result[1] /= divideBy;
            }

            long end1 = System.currentTimeMillis();
            long z1 = hessianMatrixEnd - hessianMatrixStart;
            long z2 = eigenValuesEnd - hessianMatrixEnd;
            long z3 = end1 - eigenValuesEnd;
            long end = System.currentTimeMillis();
            long z4 = end - start;
            return new long[]{z1, z2, z3, z4};
         }
      }
   }

   public float[] doubleToFloat(double[] a) {
      float[] a2 = new float[a.length];

      for(int r = 0; r < a.length; ++r) {
         a2[r] = (float)a[r];
      }

      return a2;
   }

   public double[][] floatToDouble(float[][] a) {
      double[][] a2 = new double[a.length][a[0].length];

      for(int r = 0; r < a.length; ++r) {
         for(int c = 0; c < a[0].length; ++c) {
            a2[r][c] = (double)a[r][c];
         }
      }

      return a2;
   }

   public boolean hessianEigenvaluesAtPoint2D(
      int x, int y, boolean orderOnAbsoluteSize, double[] result, boolean normalize, boolean fixUp, float sepX, float sepY
   ) {
      if (this._3D) {
         IJ.error("hessianEigenvaluesAtPoint2D( x, y, z, ... ) is only for 2D data.");
         return false;
      } else {
         ComputeCurvatures.FloatArray2D data2D = (ComputeCurvatures.FloatArray2D)this.data;
         if (fixUp) {
            if (x == 0) {
               x = 1;
            }

            if (x == data2D.width - 1) {
               x = data2D.width - 2;
            }

            if (y == 0) {
               y = 1;
            }

            if (y == data2D.height - 1) {
               y = data2D.height - 2;
            }
         }

         double[][] hessianMatrix = this.computeHessianMatrix2DDouble(data2D, x, y, this.sigma, sepX, sepY);
         double[] eigenValues = this.computeEigenValues(hessianMatrix);
         if (eigenValues == null) {
            return false;
         } else {
            double e0 = eigenValues[0];
            double e1 = eigenValues[1];
            double e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
            double e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
            if (e0c <= e1c) {
               result[0] = e0;
               result[1] = e1;
            } else {
               result[0] = e1;
               result[1] = e0;
            }

            if (normalize) {
               double divideBy = Math.abs(result[1]);
               result[0] /= divideBy;
               result[1] /= divideBy;
            }

            return true;
         }
      }
   }

   public boolean hessianEigenvaluesAtPoint3D(
      int x, int y, int z, boolean orderOnAbsoluteSize, float[] result, boolean normalize, boolean fixUp, float sepX, float sepY, float sepZ
   ) {
      if (!this._3D) {
         IJ.error("hessianEigenvaluesAtPoint3D( x, y, z, ... ) is only for 3D data.");
         return false;
      } else {
         ComputeCurvatures.FloatArray3D data3D = (ComputeCurvatures.FloatArray3D)this.data;
         if (fixUp) {
            if (x == 0) {
               x = 1;
            }

            if (x == data3D.width - 1) {
               x = data3D.width - 2;
            }

            if (y == 0) {
               y = 1;
            }

            if (y == data3D.height - 1) {
               y = data3D.height - 2;
            }

            if (z == 0) {
               z = 1;
            }

            if (z == data3D.depth - 1) {
               z = data3D.depth - 2;
            }
         }

         float[][] hessianMatrix = this.computeHessianMatrix3DFloat(data3D, x, y, z, this.sigma, sepX, sepY, sepZ);
         float[] eigenValues = this.computeEigenValues(hessianMatrix);
         if (eigenValues == null) {
            return false;
         } else {
            float e0 = eigenValues[0];
            float e1 = eigenValues[1];
            float e2 = eigenValues[2];
            float e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
            float e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
            float e2c = orderOnAbsoluteSize ? Math.abs(e2) : e2;
            if (e0c <= e1c) {
               if (e1c <= e2c) {
                  result[0] = e0;
                  result[1] = e1;
                  result[2] = e2;
               } else if (e0c <= e2c) {
                  result[0] = e0;
                  result[1] = e2;
                  result[2] = e1;
               } else {
                  result[0] = e2;
                  result[1] = e0;
                  result[2] = e1;
               }
            } else if (e0c <= e2c) {
               result[0] = e1;
               result[1] = e0;
               result[2] = e2;
            } else if (e1c <= e2c) {
               result[0] = e1;
               result[1] = e2;
               result[2] = e0;
            } else {
               result[0] = e2;
               result[1] = e1;
               result[2] = e0;
            }

            if (normalize) {
               float divideBy = Math.abs(result[2]);
               result[0] /= divideBy;
               result[1] /= divideBy;
               result[2] /= divideBy;
            }

            return true;
         }
      }
   }

   public boolean hessianEigenvaluesAtPoint3D(
      int x, int y, int z, boolean orderOnAbsoluteSize, double[] result, boolean normalize, boolean fixUp, float sepX, float sepY, float sepZ
   ) {
      if (!this._3D) {
         IJ.error("hessianEigenvaluesAtPoint3D( x, y, z, ... ) is only for 3D data.");
         return false;
      } else {
         ComputeCurvatures.FloatArray3D data3D = (ComputeCurvatures.FloatArray3D)this.data;
         if (fixUp) {
            if (x == 0) {
               x = 1;
            }

            if (x == data3D.width - 1) {
               x = data3D.width - 2;
            }

            if (y == 0) {
               y = 1;
            }

            if (y == data3D.height - 1) {
               y = data3D.height - 2;
            }

            if (z == 0) {
               z = 1;
            }

            if (z == data3D.depth - 1) {
               z = data3D.depth - 2;
            }
         }

         double[][] hessianMatrix = this.computeHessianMatrix3DDouble(data3D, x, y, z, this.sigma, sepX, sepY, sepZ);
         double[] eigenValues = this.computeEigenValues(hessianMatrix);
         if (eigenValues == null) {
            return false;
         } else {
            double e0 = eigenValues[0];
            double e1 = eigenValues[1];
            double e2 = eigenValues[2];
            double e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
            double e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
            double e2c = orderOnAbsoluteSize ? Math.abs(e2) : e2;
            if (e0c <= e1c) {
               if (e1c <= e2c) {
                  result[0] = e0;
                  result[1] = e1;
                  result[2] = e2;
               } else if (e0c <= e2c) {
                  result[0] = e0;
                  result[1] = e2;
                  result[2] = e1;
               } else {
                  result[0] = e2;
                  result[1] = e0;
                  result[2] = e1;
               }
            } else if (e0c <= e2c) {
               result[0] = e1;
               result[1] = e0;
               result[2] = e2;
            } else if (e1c <= e2c) {
               result[0] = e1;
               result[1] = e2;
               result[2] = e0;
            } else {
               result[0] = e2;
               result[1] = e1;
               result[2] = e0;
            }

            if (normalize) {
               double divideBy = Math.abs(result[2]);
               result[0] /= divideBy;
               result[1] /= divideBy;
               result[2] /= divideBy;
            }

            return true;
         }
      }
   }

   public static ImagePlus FloatArrayToImagePlus(ComputeCurvatures.FloatArray2D image, String name, float min, float max) {
      ImagePlus imp = IJ.createImage(name, "32-Bit Black", image.width, image.height, 1);
      FloatProcessor ip = (FloatProcessor)imp.getProcessor();
      FloatArrayToFloatProcessor(ip, image);
      if (min == max) {
         ip.resetMinAndMax();
      } else {
         ip.setMinAndMax((double)min, (double)max);
      }

      imp.updateAndDraw();
      return imp;
   }

   public static void FloatArrayToFloatProcessor(ImageProcessor ip, ComputeCurvatures.FloatArray2D pixels) {
      float[] data = new float[pixels.width * pixels.height];
      int count = 0;

      for(int y = 0; y < pixels.height; ++y) {
         for(int x = 0; x < pixels.width; ++x) {
            data[count] = pixels.data[count++];
         }
      }

      ip.setPixels(data);
      ip.resetMinAndMax();
   }

   public ImagePlus FloatArrayToStack(ComputeCurvatures.FloatArray3D image, String name, float min, float max) {
      int width = image.width;
      int height = image.height;
      int nstacks = image.depth;
      ImageStack stack = new ImageStack(width, height);

      for(int slice = 0; slice < nstacks; ++slice) {
         ImagePlus impResult = IJ.createImage(name, "32-Bit Black", width, height, 1);
         ImageProcessor ipResult = impResult.getProcessor();
         float[] sliceImg = new float[width * height];

         for(int x = 0; x < width; ++x) {
            for(int y = 0; y < height; ++y) {
               sliceImg[y * width + x] = image.get(x, y, slice);
            }
         }

         ipResult.setPixels(sliceImg);
         if (min == max) {
            ipResult.resetMinAndMax();
         } else {
            ipResult.setMinAndMax((double)min, (double)max);
         }

         stack.addSlice("Slice " + slice, ipResult);
      }

      return new ImagePlus(name, stack);
   }

   public double[] computeEigenValues(double[][] matrix) {
      if (matrix.length == 3 && matrix[0].length == 3) {
         Eigensystem3x3Double e = new Eigensystem3x3Double(matrix);
         boolean result = e.findEvalues();
         return result ? e.getEvaluesCopy() : null;
      } else if (matrix.length == 2 && matrix[0].length == 2) {
         Eigensystem2x2Double e = new Eigensystem2x2Double(matrix);
         boolean result = e.findEvalues();
         return result ? e.getEvaluesCopy() : null;
      } else {
         JacobiDouble jc = new JacobiDouble(matrix, 50);
         return jc.getEigenValues();
      }
   }

   public float[] computeEigenValues(float[][] matrix) {
      if (matrix.length == 3 && matrix[0].length == 3) {
         Eigensystem3x3Float e = new Eigensystem3x3Float(matrix);
         boolean result = e.findEvalues();
         return result ? e.getEvaluesCopy() : null;
      } else if (matrix.length == 2 && matrix[0].length == 2) {
         Eigensystem2x2Float e = new Eigensystem2x2Float(matrix);
         boolean result = e.findEvalues();
         return result ? e.getEvaluesCopy() : null;
      } else {
         JacobiFloat jc = new JacobiFloat(matrix, 50);
         return jc.getEigenValues();
      }
   }

   public double[][] computeHessianMatrix2DDouble(ComputeCurvatures.FloatArray2D laPlace, int x, int y, double sigma, float sepX, float sepY) {
      if (laPlace == null) {
         laPlace = (ComputeCurvatures.FloatArray2D)this.data;
      }

      double[][] hessianMatrix = new double[2][2];
      double temp = (double)(2.0F * laPlace.get(x, y));
      hessianMatrix[0][0] = (double)laPlace.get(x + 1, y) - temp + (double)laPlace.get(x - 1, y);
      hessianMatrix[1][1] = (double)laPlace.get(x, y + 1) - temp + (double)laPlace.get(x, y - 1);
      hessianMatrix[0][1] = hessianMatrix[1][0] = (double)(
         ((laPlace.get(x + 1, y + 1) - laPlace.get(x - 1, y + 1)) / 2.0F - (laPlace.get(x + 1, y - 1) - laPlace.get(x - 1, y - 1)) / 2.0F) / 2.0F
      );

      for(int i = 0; i < 2; ++i) {
         for(int j = 0; j < 2; ++j) {
            hessianMatrix[i][j] *= sigma * sigma;
         }
      }

      return hessianMatrix;
   }

   public float[][] computeHessianMatrix2DFloat(ComputeCurvatures.FloatArray2D laPlace, int x, int y, double sigma, float sepX, float sepY) {
      if (laPlace == null) {
         laPlace = (ComputeCurvatures.FloatArray2D)this.data;
      }

      float[][] hessianMatrix = new float[2][2];
      float temp = 2.0F * laPlace.get(x, y);
      hessianMatrix[0][0] = laPlace.get(x + 1, y) - temp + laPlace.get(x - 1, y);
      hessianMatrix[1][1] = laPlace.get(x, y + 1) - temp + laPlace.get(x, y - 1);
      hessianMatrix[0][1] = hessianMatrix[1][0] = (
            (laPlace.get(x + 1, y + 1) - laPlace.get(x - 1, y + 1)) / 2.0F - (laPlace.get(x + 1, y - 1) - laPlace.get(x - 1, y - 1)) / 2.0F
         )
         / 2.0F;

      for(int i = 0; i < 2; ++i) {
         for(int j = 0; j < 2; ++j) {
            hessianMatrix[i][j] = (float)((double)hessianMatrix[i][j] * sigma * sigma);
         }
      }

      return hessianMatrix;
   }

   public double[][] computeHessianMatrix3DDouble(ComputeCurvatures.FloatArray3D img, int x, int y, int z, double sigma, float sepX, float sepY, float sepZ) {
      if (img == null) {
         img = (ComputeCurvatures.FloatArray3D)this.data;
      }

      double[][] hessianMatrix = new double[3][3];
      double temp = (double)(2.0F * img.get(x, y, z));
      hessianMatrix[0][0] = (double)img.get(x + 1, y, z) - temp + (double)img.get(x - 1, y, z);
      hessianMatrix[1][1] = (double)img.get(x, y + 1, z) - temp + (double)img.get(x, y - 1, z);
      hessianMatrix[2][2] = (double)img.get(x, y, z + 1) - temp + (double)img.get(x, y, z - 1);
      hessianMatrix[0][1] = hessianMatrix[1][0] = (double)(
         ((img.get(x + 1, y + 1, z) - img.get(x - 1, y + 1, z)) / 2.0F - (img.get(x + 1, y - 1, z) - img.get(x - 1, y - 1, z)) / 2.0F) / 2.0F
      );
      hessianMatrix[0][2] = hessianMatrix[2][0] = (double)(
         ((img.get(x + 1, y, z + 1) - img.get(x - 1, y, z + 1)) / 2.0F - (img.get(x + 1, y, z - 1) - img.get(x - 1, y, z - 1)) / 2.0F) / 2.0F
      );
      hessianMatrix[1][2] = hessianMatrix[2][1] = (double)(
         ((img.get(x, y + 1, z + 1) - img.get(x, y - 1, z + 1)) / 2.0F - (img.get(x, y + 1, z - 1) - img.get(x, y - 1, z - 1)) / 2.0F) / 2.0F
      );

      for(int i = 0; i < 3; ++i) {
         for(int j = 0; j < 3; ++j) {
            hessianMatrix[i][j] *= sigma * sigma;
         }
      }

      return hessianMatrix;
   }

   public float[][] computeHessianMatrix3DFloat(ComputeCurvatures.FloatArray3D img, int x, int y, int z, double sigma, float sepX, float sepY, float sepZ) {
      if (img == null) {
         img = (ComputeCurvatures.FloatArray3D)this.data;
      }

      float[][] hessianMatrix = new float[3][3];
      float temp = 2.0F * img.get(x, y, z);
      hessianMatrix[0][0] = img.get(x + 1, y, z) - temp + img.get(x - 1, y, z);
      hessianMatrix[1][1] = img.get(x, y + 1, z) - temp + img.get(x, y - 1, z);
      hessianMatrix[2][2] = img.get(x, y, z + 1) - temp + img.get(x, y, z - 1);
      hessianMatrix[0][1] = hessianMatrix[1][0] = (
            (img.get(x + 1, y + 1, z) - img.get(x - 1, y + 1, z)) / 2.0F - (img.get(x + 1, y - 1, z) - img.get(x - 1, y - 1, z)) / 2.0F
         )
         / 2.0F;
      hessianMatrix[0][2] = hessianMatrix[2][0] = (
            (img.get(x + 1, y, z + 1) - img.get(x - 1, y, z + 1)) / 2.0F - (img.get(x + 1, y, z - 1) - img.get(x - 1, y, z - 1)) / 2.0F
         )
         / 2.0F;
      hessianMatrix[1][2] = hessianMatrix[2][1] = (
            (img.get(x, y + 1, z + 1) - img.get(x, y - 1, z + 1)) / 2.0F - (img.get(x, y + 1, z - 1) - img.get(x, y - 1, z - 1)) / 2.0F
         )
         / 2.0F;

      for(int i = 0; i < 3; ++i) {
         for(int j = 0; j < 3; ++j) {
            hessianMatrix[i][j] = (float)((double)hessianMatrix[i][j] * sigma * sigma);
         }
      }

      return hessianMatrix;
   }

   public static float[] createGaussianKernel1D(float sigma, boolean normalize) {
      float[] gaussianKernel;
      if (sigma <= 0.0F) {
         gaussianKernel = new float[]{0.0F, 1.0F, 0.0F};
      } else {
         int size = Math.max(3, 2 * (int)((double)(3.0F * sigma) + 0.5) + 1);
         float two_sq_sigma = 2.0F * sigma * sigma;
         gaussianKernel = new float[size];

         for(int x = size / 2; x >= 0; --x) {
            float val = (float)Math.exp((double)(-((float)(x * x)) / two_sq_sigma));
            gaussianKernel[size / 2 - x] = val;
            gaussianKernel[size / 2 + x] = val;
         }
      }

      if (normalize) {
         float sum = 0.0F;

         for(int i = 0; i < gaussianKernel.length; ++i) {
            sum += gaussianKernel[i];
         }

         for(int i = 0; i < gaussianKernel.length; ++i) {
            gaussianKernel[i] /= sum;
         }
      }

      return gaussianKernel;
   }

   public ComputeCurvatures.FloatArray2D computeGaussianFastMirror(
      ComputeCurvatures.FloatArray2D input, float sigma, GaussianGenerationCallback callback, Calibration calibration
   ) {
      long s0 = System.currentTimeMillis();
      ComputeCurvatures.FloatArray2D output = new ComputeCurvatures.FloatArray2D(input.width, input.height);
      float kernelsumX = 0.0F;
      float kernelsumY = 0.0F;
      float kernelsumZ = 0.0F;
      float pixelWidth = 1.0F;
      float pixelHeight = 1.0F;
      float pixelDepth = 1.0F;
      if (calibration != null) {
         pixelWidth = (float)calibration.pixelWidth;
         pixelHeight = (float)calibration.pixelHeight;
         pixelDepth = (float)calibration.pixelDepth;
      }

      long s1 = System.currentTimeMillis();
      float[] kernelX = createGaussianKernel1D(sigma / pixelWidth, true);
      float[] kernelY = createGaussianKernel1D(sigma / pixelHeight, true);
      int filterSizeX = kernelX.length;
      int filterSizeY = kernelY.length;
      long s2 = System.currentTimeMillis();

      for(int i = 0; i < kernelX.length; ++i) {
         kernelsumX += kernelX[i];
      }

      for(int i = 0; i < kernelY.length; ++i) {
         kernelsumY += kernelY[i];
      }

      long s3 = System.currentTimeMillis();
      double totalPoints = (double)(input.width * input.height * 2);
      long pointsDone = 0L;

      for(int x = 0; x < input.width; ++x) {
         if (this.cancelGeneration) {
            return null;
         }

         for(int y = 0; y < input.height; ++y) {
            float avg = 0.0F;
            if (x - filterSizeX / 2 >= 0 && x + filterSizeX / 2 < input.width) {
               for(int f = -filterSizeX / 2; f <= filterSizeX / 2; ++f) {
                  avg += input.get(x + f, y) * kernelX[f + filterSizeX / 2];
               }
            } else {
               for(int f = -filterSizeX / 2; f <= filterSizeX / 2; ++f) {
                  avg += input.getMirror(x + f, y) * kernelX[f + filterSizeX / 2];
               }
            }

            output.set(avg / kernelsumX, x, y);
         }

         pointsDone += (long)input.height;
         if (callback != null) {
            callback.proportionDone((double)pointsDone / totalPoints);
         }
      }

      long s4 = System.currentTimeMillis();

      for(int x = 0; x < input.width; ++x) {
         if (this.cancelGeneration) {
            return null;
         }

         float[] temp = new float[input.height];

         for(int y = 0; y < input.height; ++y) {
            float avg = 0.0F;
            if (y - filterSizeY / 2 >= 0 && y + filterSizeY / 2 < input.height) {
               for(int f = -filterSizeY / 2; f <= filterSizeY / 2; ++f) {
                  avg += output.get(x, y + f) * kernelY[f + filterSizeY / 2];
               }
            } else {
               for(int f = -filterSizeY / 2; f <= filterSizeY / 2; ++f) {
                  avg += output.getMirror(x, y + f) * kernelY[f + filterSizeY / 2];
               }
            }

            temp[y] = avg / kernelsumY;
         }

         for(int y = 0; y < input.height; ++y) {
            output.set(temp[y], x, y);
         }

         pointsDone += (long)input.height;
         if (callback != null) {
            callback.proportionDone((double)pointsDone / totalPoints);
         }
      }

      long s5 = System.currentTimeMillis();
      if (callback != null) {
         callback.proportionDone(1.0);
      }

      if (!Utils.isReleaseVersion) {
         System.err.println("s1=" + (s1 - s0) + "ms\ts2=" + (s2 - s3) + "ms\ts3=" + (s4 - s3) + "ms\ts4=" + (s5 - s4));
      }

      return output;
   }

   public ComputeCurvatures.FloatArray3D computeGaussianFastMirror(
      ComputeCurvatures.FloatArray3D input, float sigma, GaussianGenerationCallback callback, Calibration calibration
   ) {
      ComputeCurvatures.FloatArray3D output = new ComputeCurvatures.FloatArray3D(input.width, input.height, input.depth);
      float kernelsumX = 0.0F;
      float kernelsumY = 0.0F;
      float kernelsumZ = 0.0F;
      float pixelWidth = 1.0F;
      float pixelHeight = 1.0F;
      float pixelDepth = 1.0F;
      if (calibration != null) {
         pixelWidth = (float)calibration.pixelWidth;
         pixelHeight = (float)calibration.pixelHeight;
         pixelDepth = (float)calibration.pixelDepth;
      }

      float[] kernelX = createGaussianKernel1D(sigma / pixelWidth, true);
      float[] kernelY = createGaussianKernel1D(sigma / pixelHeight, true);
      float[] kernelZ = createGaussianKernel1D(sigma / pixelDepth, true);
      int filterSizeX = kernelX.length;
      int filterSizeY = kernelY.length;
      int filterSizeZ = kernelZ.length;

      for(int i = 0; i < kernelX.length; ++i) {
         kernelsumX += kernelX[i];
      }

      for(int i = 0; i < kernelY.length; ++i) {
         kernelsumY += kernelY[i];
      }

      for(int i = 0; i < kernelZ.length; ++i) {
         kernelsumZ += kernelZ[i];
      }

      double totalPoints = (double)(input.width * input.height * input.depth * 3);
      long pointsDone = 0L;

      for(int x = 0; x < input.width; ++x) {
         if (this.cancelGeneration) {
            return null;
         }

         for(int y = 0; y < input.height; ++y) {
            for(int z = 0; z < input.depth; ++z) {
               float avg = 0.0F;
               if (x - filterSizeX / 2 >= 0 && x + filterSizeX / 2 < input.width) {
                  for(int f = -filterSizeX / 2; f <= filterSizeX / 2; ++f) {
                     avg += input.get(x + f, y, z) * kernelX[f + filterSizeX / 2];
                  }
               } else {
                  for(int f = -filterSizeX / 2; f <= filterSizeX / 2; ++f) {
                     avg += input.getMirror(x + f, y, z) * kernelX[f + filterSizeX / 2];
                  }
               }

               output.set(avg / kernelsumX, x, y, z);
            }
         }

         pointsDone += (long)(input.height * input.depth);
         if (callback != null) {
            callback.proportionDone((double)pointsDone / totalPoints);
         }
      }

      for(int x = 0; x < input.width; ++x) {
         if (this.cancelGeneration) {
            return null;
         }

         for(int z = 0; z < input.depth; ++z) {
            float[] temp = new float[input.height];

            for(int y = 0; y < input.height; ++y) {
               float avg = 0.0F;
               if (y - filterSizeY / 2 >= 0 && y + filterSizeY / 2 < input.height) {
                  for(int f = -filterSizeY / 2; f <= filterSizeY / 2; ++f) {
                     avg += output.get(x, y + f, z) * kernelY[f + filterSizeY / 2];
                  }
               } else {
                  for(int f = -filterSizeY / 2; f <= filterSizeY / 2; ++f) {
                     avg += output.getMirror(x, y + f, z) * kernelY[f + filterSizeY / 2];
                  }
               }

               temp[y] = avg / kernelsumY;
            }

            for(int y = 0; y < input.height; ++y) {
               output.set(temp[y], x, y, z);
            }
         }

         pointsDone += (long)(input.depth * input.height);
         if (callback != null) {
            callback.proportionDone((double)pointsDone / totalPoints);
         }
      }

      for(int x = 0; x < input.width; ++x) {
         if (this.cancelGeneration) {
            return null;
         }

         for(int y = 0; y < input.height; ++y) {
            float[] temp = new float[input.depth];

            for(int z = 0; z < input.depth; ++z) {
               float avg = 0.0F;
               if (z - filterSizeZ / 2 >= 0 && z + filterSizeZ / 2 < input.depth) {
                  for(int f = -filterSizeZ / 2; f <= filterSizeZ / 2; ++f) {
                     avg += output.get(x, y, z + f) * kernelZ[f + filterSizeZ / 2];
                  }
               } else {
                  for(int f = -filterSizeZ / 2; f <= filterSizeZ / 2; ++f) {
                     avg += output.getMirror(x, y, z + f) * kernelZ[f + filterSizeZ / 2];
                  }
               }

               temp[z] = avg / kernelsumZ;
            }

            for(int z = 0; z < input.depth; ++z) {
               output.set(temp[z], x, y, z);
            }
         }

         pointsDone += (long)(input.height * input.depth);
         if (callback != null) {
            callback.proportionDone((double)pointsDone / totalPoints);
         }
      }

      if (callback != null) {
         callback.proportionDone(1.0);
      }

      return output;
   }

   public ComputeCurvatures.FloatArray3D StackToFloatArray(ImageStack stack) {
      Object[] imageStack = stack.getImageArray();
      int width = stack.getWidth();
      int height = stack.getHeight();
      int nstacks = stack.getSize();
      if (imageStack == null || imageStack.length == 0) {
         IJ.error("Image Stack is empty.");
         return null;
      } else if (imageStack[0] instanceof int[]) {
         IJ.error("RGB images not supported at the moment.");
         return null;
      } else {
         ComputeCurvatures.FloatArray3D pixels = new ComputeCurvatures.FloatArray3D(width, height, nstacks);
         if (imageStack[0] instanceof byte[]) {
            for(int countSlice = 0; countSlice < nstacks; ++countSlice) {
               byte[] pixelTmp = (byte[])imageStack[countSlice];
               int count = 0;

               for(int y = 0; y < height; ++y) {
                  for(int x = 0; x < width; ++x) {
                     pixels.data[pixels.getPos(x, y, countSlice)] = (float)(pixelTmp[count++] & 255);
                  }
               }
            }
         } else if (imageStack[0] instanceof short[]) {
            for(int countSlice = 0; countSlice < nstacks; ++countSlice) {
               short[] pixelTmp = (short[])imageStack[countSlice];
               int count = 0;

               for(int y = 0; y < height; ++y) {
                  for(int x = 0; x < width; ++x) {
                     pixels.data[pixels.getPos(x, y, countSlice)] = (float)(pixelTmp[count++] & '\uffff');
                  }
               }
            }
         } else {
            for(int countSlice = 0; countSlice < nstacks; ++countSlice) {
               float[] pixelTmp = (float[])imageStack[countSlice];
               int count = 0;

               for(int y = 0; y < height; ++y) {
                  for(int x = 0; x < width; ++x) {
                     pixels.data[pixels.getPos(x, y, countSlice)] = pixelTmp[count++];
                  }
               }
            }
         }

         return pixels;
      }
   }

   public ComputeCurvatures.FloatArray2D ImageToFloatArray(ImageProcessor ip) {
      Object pixelArray = ip.getPixels();
      int count = 0;
      ComputeCurvatures.FloatArray2D image;
      if (ip instanceof ByteProcessor) {
         image = new ComputeCurvatures.FloatArray2D(ip.getWidth(), ip.getHeight());
         byte[] pixels = (byte[])pixelArray;

         for(int y = 0; y < ip.getHeight(); ++y) {
            for(int x = 0; x < ip.getWidth(); ++x) {
               image.data[count] = (float)(pixels[count++] & 255);
            }
         }
      } else if (ip instanceof ShortProcessor) {
         image = new ComputeCurvatures.FloatArray2D(ip.getWidth(), ip.getHeight());
         short[] pixels = (short[])pixelArray;

         for(int y = 0; y < ip.getHeight(); ++y) {
            for(int x = 0; x < ip.getWidth(); ++x) {
               image.data[count] = (float)(pixels[count++] & '\uffff');
            }
         }
      } else if (ip instanceof FloatProcessor) {
         image = new ComputeCurvatures.FloatArray2D(ip.getWidth(), ip.getHeight());
         float[] pixels = (float[])pixelArray;

         for(int y = 0; y < ip.getHeight(); ++y) {
            for(int x = 0; x < ip.getWidth(); ++x) {
               image.data[count] = pixels[count++];
            }
         }
      } else {
         System.err.println("RGB images not supported");
         image = null;
      }

      return image;
   }

   public abstract class FloatArray {
      public float[] data = null;

      public abstract ComputeCurvatures.FloatArray clone();
   }

   public class FloatArray2D extends ComputeCurvatures.FloatArray {
      public float[] data = null;
      public int width = 0;
      public int height = 0;

      public FloatArray2D(int width, int height) {
         this.data = new float[width * height];
         this.width = width;
         this.height = height;
      }

      public FloatArray2D(float[] data, int width, int height) {
         this.data = data;
         this.width = width;
         this.height = height;
      }

      public ComputeCurvatures.FloatArray2D clone() {
         ComputeCurvatures.FloatArray2D clone = ComputeCurvatures.this.new FloatArray2D(this.width, this.height);
         System.arraycopy(this.data, 0, clone.data, 0, this.data.length);
         return clone;
      }

      public int getPos(int x, int y) {
         return x + this.width * y;
      }

      public float get(int x, int y) {
         return this.data[this.getPos(x, y)];
      }

      public float getMirror(int x, int y) {
         if (x >= this.width) {
            x = this.width - (x - this.width + 2);
         }

         if (y >= this.height) {
            y = this.height - (y - this.height + 2);
         }

         if (x < 0) {
            int tmp = 0;

            for(int dir = 1; x < 0; ++x) {
               tmp += dir;
               if (tmp == this.width - 1 || tmp == 0) {
                  dir *= -1;
               }
            }

            x = tmp;
         }

         if (y < 0) {
            int tmp = 0;

            for(int dir = 1; y < 0; ++y) {
               tmp += dir;
               if (tmp == this.height - 1 || tmp == 0) {
                  dir *= -1;
               }
            }

            y = tmp;
         }

         return this.data[this.getPos(x, y)];
      }

      public float getZero(int x, int y) {
         if (x >= this.width) {
            return 0.0F;
         } else if (y >= this.height) {
            return 0.0F;
         } else if (x < 0) {
            return 0.0F;
         } else {
            return y < 0 ? 0.0F : this.data[this.getPos(x, y)];
         }
      }

      public void set(float value, int x, int y) {
         this.data[this.getPos(x, y)] = value;
      }
   }

   public class FloatArray3D extends ComputeCurvatures.FloatArray {
      public float[] data = null;
      public int width = 0;
      public int height = 0;
      public int depth = 0;

      public FloatArray3D(float[] data, int width, int height, int depth) {
         this.data = data;
         this.width = width;
         this.height = height;
         this.depth = depth;
      }

      public FloatArray3D(int width, int height, int depth) {
         this.data = new float[width * height * depth];
         this.width = width;
         this.height = height;
         this.depth = depth;
      }

      public ComputeCurvatures.FloatArray3D clone() {
         ComputeCurvatures.FloatArray3D clone = ComputeCurvatures.this.new FloatArray3D(this.width, this.height, this.depth);
         System.arraycopy(this.data, 0, clone.data, 0, this.data.length);
         return clone;
      }

      public int getPos(int x, int y, int z) {
         return x + this.width * (y + z * this.height);
      }

      public float get(int x, int y, int z) {
         return this.data[this.getPos(x, y, z)];
      }

      public float getMirror(int x, int y, int z) {
         if (x >= this.width) {
            x = this.width - (x - this.width + 2);
         }

         if (y >= this.height) {
            y = this.height - (y - this.height + 2);
         }

         if (z >= this.depth) {
            z = this.depth - (z - this.depth + 2);
         }

         if (x < 0) {
            int tmp = 0;

            for(int dir = 1; x < 0; ++x) {
               tmp += dir;
               if (tmp == this.width - 1 || tmp == 0) {
                  dir *= -1;
               }
            }

            x = tmp;
         }

         if (y < 0) {
            int tmp = 0;

            for(int dir = 1; y < 0; ++y) {
               tmp += dir;
               if (tmp == this.height - 1 || tmp == 0) {
                  dir *= -1;
               }
            }

            y = tmp;
         }

         if (z < 0) {
            int tmp = 0;

            for(int dir = 1; z < 0; ++z) {
               tmp += dir;
               if (tmp == this.height - 1 || tmp == 0) {
                  dir *= -1;
               }
            }

            z = tmp;
         }

         return this.data[this.getPos(x, y, z)];
      }

      public void set(float value, int x, int y, int z) {
         this.data[this.getPos(x, y, z)] = value;
      }

      public ComputeCurvatures.FloatArray2D getXPlane(int x) {
         ComputeCurvatures.FloatArray2D plane = ComputeCurvatures.this.new FloatArray2D(this.height, this.depth);

         for(int y = 0; y < this.height; ++y) {
            for(int z = 0; z < this.depth; ++z) {
               plane.set(this.get(x, y, z), y, z);
            }
         }

         return plane;
      }

      public float[][] getXPlane_float(int x) {
         float[][] plane = new float[this.height][this.depth];

         for(int y = 0; y < this.height; ++y) {
            for(int z = 0; z < this.depth; ++z) {
               plane[y][z] = this.get(x, y, z);
            }
         }

         return plane;
      }

      public ComputeCurvatures.FloatArray2D getYPlane(int y) {
         ComputeCurvatures.FloatArray2D plane = ComputeCurvatures.this.new FloatArray2D(this.width, this.depth);

         for(int x = 0; x < this.width; ++x) {
            for(int z = 0; z < this.depth; ++z) {
               plane.set(this.get(x, y, z), x, z);
            }
         }

         return plane;
      }

      public float[][] getYPlane_float(int y) {
         float[][] plane = new float[this.width][this.depth];

         for(int x = 0; x < this.width; ++x) {
            for(int z = 0; z < this.depth; ++z) {
               plane[x][z] = this.get(x, y, z);
            }
         }

         return plane;
      }

      public ComputeCurvatures.FloatArray2D getZPlane(int z) {
         ComputeCurvatures.FloatArray2D plane = ComputeCurvatures.this.new FloatArray2D(this.width, this.height);

         for(int x = 0; x < this.width; ++x) {
            for(int y = 0; y < this.height; ++y) {
               plane.set(this.get(x, y, z), x, y);
            }
         }

         return plane;
      }

      public float[][] getZPlane_float(int z) {
         float[][] plane = new float[this.width][this.height];

         for(int x = 0; x < this.width; ++x) {
            for(int y = 0; y < this.height; ++y) {
               plane[x][y] = this.get(x, y, z);
            }
         }

         return plane;
      }

      public void setXPlane(ComputeCurvatures.FloatArray2D plane, int x) {
         for(int y = 0; y < this.height; ++y) {
            for(int z = 0; z < this.depth; ++z) {
               this.set(plane.get(y, z), x, y, z);
            }
         }
      }

      public void setXPlane(float[][] plane, int x) {
         for(int y = 0; y < this.height; ++y) {
            for(int z = 0; z < this.depth; ++z) {
               this.set(plane[y][z], x, y, z);
            }
         }
      }

      public void setYPlane(ComputeCurvatures.FloatArray2D plane, int y) {
         for(int x = 0; x < this.width; ++x) {
            for(int z = 0; z < this.depth; ++z) {
               this.set(plane.get(x, z), x, y, z);
            }
         }
      }

      public void setYPlane(float[][] plane, int y) {
         for(int x = 0; x < this.width; ++x) {
            for(int z = 0; z < this.depth; ++z) {
               this.set(plane[x][z], x, y, z);
            }
         }
      }

      public void setZPlane(ComputeCurvatures.FloatArray2D plane, int z) {
         for(int x = 0; x < this.width; ++x) {
            for(int y = 0; y < this.height; ++y) {
               this.set(plane.get(x, y), x, y, z);
            }
         }
      }

      public void setZPlane(float[][] plane, int z) {
         for(int x = 0; x < this.width; ++x) {
            for(int y = 0; y < this.height; ++y) {
               this.set(plane[x][y], x, y, z);
            }
         }
      }
   }

   static class TrivialProgressDisplayer implements GaussianGenerationCallback {
      @Override
      public void proportionDone(double proportion) {
         if (proportion < 0.0) {
            IJ.showProgress(1.0);
         } else {
            IJ.showProgress(proportion);
         }
      }
   }
}