package com.github.celldynamics.quimp.plugin.dic;
   
   import org.slf4j.Logger;
   import org.slf4j.LoggerFactory;
   
   import com.github.celldynamics.quimp.plugin.utils.ImageProcessorPlus;
   
   import ij.ImagePlus;
   import ij.process.FloatProcessor;
  import ij.process.ImageProcessor;
  import ij.process.ImageStatistics;
  
  /**
   * Implementation of Line Integration and Deconvolution algorithm proposed by Kam.
   * 
   * <p><h1>Principles</h1>
   * This algorithm uses corrected line integration that runs in both directions from current point of
   * image (<i>x_n,y_n</i>) to <i>r_1</i> and end <i>r_2</i>. Vector <i>dr</i> is set to be parallel
   * to 0X axis. Final formula for reconstruction of pixel (<i>x_n,y_n</i>) is given by:<br>
   * <img src="doc-files/eq.png"/><br>
   * 
   * <p><h1>Algorithm</h1>
   * The algorithm perform the following steps to reconstruct the whole image:
   * <ol>
   * <li>Input is converted to 16bit and stored in private field ExtraImageProcessor
   * srcImageCopyProcessor
   * <li>The value of shift is added to all pixels of input image. All operations are performed on
   * copy of input ImageProcessor or ImagePlus. Rotate image by angle to bring bas-reliefs
   * perpendicular to OX axis. Rotated image has different dimensions and it is padded by 0. After
   * this operation 0 pixels are those that belong to background. For every line in rotated image
   * position of true pixels is calculated. True pixels are those that belong to original image
   * excluding background added during rotation. For every line position of first and last true pixel
   * is noted in table ranges
   * <li>Decay factors are pre-calculated and stored in decays table.
   * <li>Final reconstruction is performed.
   * </ol>
   * Image for reconstruction is passed during construction of DICReconstruction object. For this
   * object ranges and decays are evaluated and then user can call reconstructionDicLid() method to
   * get reconstruction. getRanges() method also rotates private object thus rotation in
   * reconstructionDicLid() is not necessary (flagged by DICReconstruction::isRotated). Different
   * situation happens when the whole stack is reconstructed. To prevent creating new instance of
   * DICReconstruction for every slice the setIp(ImageProcessor) method is used for connecting new
   * slice. In this case it is assumed that ImageProcessor objects are similar and they have the same
   * geometry. ranges are filled only once on DICReconstruction constructing thus images connected by
   * setIp(ImageProcessor) are not rotated. This situation is detected in
   * {@link #reconstructionDicLid()} by <tt>isRotated</tt> flag.
   * 
   * <p>Privates:
   * <ul>
   * <li><i>ranges</i> - true pixels begin and end on x axis. Set by getRanges(). [r][0] - x of first
   * pixel of line r of image, [r][1] - x of last pixel of image of line r
   * <li><i>maxWidth</i> - Set by getRanges()
   * <li><i>decays</i> - set by generateDeacy(double, int)
   * <li><i>srcImageCopyProcessor</i> - local <b>copy</b> of input ImageProcessor passed to object.
   * Set by constructors and setIp(ImageProcessor)
   * <li><i>isRotated</i> - It is set by getRanges() that rotates object to get true pixels position
   * and cancelled by setIP
   * </ul>
   * 
   * <p><b>Warning</b>
   * Currently this class supports only 8bit images. It can support also 16bit input but in this case
   * algorithm used for detection of true pixels may not work correctly for certain cases - when
   * maximum intensity will be max-shift
   * 
   * <p>The input image can be prefitlered before processing. This is running mean filter of given
   * mask
   * size applied at angle perpendicular to the shear (this angle is given by caller).
   * 
   * <p><h1>References</h1>
   * <ul>
   * <li>Z. Kam, "Microscopic differential interference contrast image processing by line integration
   * (LID) and deconvolution," Bioimaging, vol. 6, no. 4, pp. 166–176, 1998.
   * <li>B. Heise, A. Sonnleitner, and E. P. Klement, "DIC image reconstruction on large cell scans.,"
   * Microsc. Res. Tech., vol. 66, no. 6, pp. 312–320, 2005.
   * </ul>
   * 
   * @author p.baniukiewicz
   *
   */
  public class LidReconstructor {
  
    /**
     * The Constant LOGGER.
     */
    static final Logger LOGGER = LoggerFactory.getLogger(LidReconstructor.class.getName());
    private final int shift = 1; // shift added to original image eliminate 0s
  
    private ImageProcessor srcIp; // local reference of ImageProcessor (const)
    private double decay;
    private double angle;
    private double[] decays; // reference to preallocated decay data created
    private int maxWidth; // Width of image after rotation. Set by getRanges()
    private int[][] ranges; // true pixels begin and end on x axis.
    private ImageProcessor srcImageCopyProcessor; // local copy of input
    private boolean isRotated; // true if srcImageCopyProcessor is rotated
    private ImageStatistics is;
    private ImageProcessorPlus ipp; // helper class for rotating images
    private String prefilterangle;
    private int masksize;
 
   /**
    * Default constructor that accepts ImagePlus. It does not support stacks.
    * 
    * <p>Input srcImage is not modified
    * 
    * @param srcImage Image to process
    * @param decay algorithm constant
    * @param angle DIC angle
    * 
    * @throws DicException Throws exception after generateRanges()
    */
   public LidReconstructor(final ImagePlus srcImage, double decay, double angle)
           throws DicException {
     this(srcImage.getProcessor(), decay, angle);
   }
 
   /**
    * Default constructor that accepts ImageProcessor.
    * 
    * <p>Input ip is not modified
    * 
    * @param ip ImageProcessor to process
    * @param decay algorithm constant
    * @param angle DIC angle
    * 
    * @throws DicException Throws exception after generateRanges()
    */
   public LidReconstructor(final ImageProcessor ip, double decay, double angle) throws DicException {
     this(ip, decay, angle, "0", 0);
   }
 
   /**
    * Default constructor.
    * 
    * @param ip image to process
    * @param decay algorithm constant
    * @param angle DIC angle
    * @param prefilterangle Supported angle of prefiltering
    * @param masksize uneven mask size, 0 switches off filtering
    * @throws DicException wrong image
    */
   public LidReconstructor(final ImageProcessor ip, double decay, double angle,
           String prefilterangle, int masksize) throws DicException {
     LOGGER.trace("Use params: ip:" + ip + " decay:" + decay + " angle:" + angle + " filterangle:"
             + prefilterangle + " masksize:" + masksize);
     this.angle = angle;
     this.decay = decay;
     this.isRotated = false;
     this.prefilterangle = prefilterangle;
     this.masksize = masksize;
     ipp = new ImageProcessorPlus();
     setIp(ip);
     recalculate();
   }
 
   /**
    * Sets new reconstruction parameters for current object.
    * 
    * @param decay algorithm constant
    * @param angle DIC angle
    * @throws DicException Throws exception after generateRanges()
    */
   public void setParams(double decay, double angle) throws DicException {
     this.angle = angle;
     this.decay = decay;
     recalculate();
   }
 
   /**
    * Assigns ImageProcessor for reconstruction to current object. Releases previous one.
    * 
    * <p>This method can be used for changing image connected to DICReconstruction object. New image
    * should have the same architecture as image passed in constructor. Typically this method is
    * used for passing next slice from stack.
    * 
    * <p>Input ip is not modified
    * 
    * @param ip New ImageProcessor containing image for reconstruction.
    */
   public void setIp(final ImageProcessor ip) {
     this.srcIp = ip;
     // make copy of original image to not modify it - converting to 16bit
     this.srcImageCopyProcessor = srcIp.convertToShort(true);
     new ImageProcessorPlus().runningMean(srcImageCopyProcessor, prefilterangle, masksize);
     // ensure that minmax will be recalculated (usually they are stored in class field) set
     // interpolation
     srcImageCopyProcessor.resetMinAndMax();
     srcImageCopyProcessor.setInterpolationMethod(ImageProcessor.BICUBIC);
     // Rotating image - set 0 background
     srcImageCopyProcessor.setBackgroundValue(0.0);
     // getting mean value
     is = srcImageCopyProcessor.getStatistics();
     this.isRotated = false; // new Processor not rotated yet
   }
 
   /**
    * Recalculates true pixels range and new size of image after rotation. Setup private class
    * fields.
    * 
    * <p>Modifies private class fields: <tt>maxWidth</tt>, <tt>ranges</tt>,
    * <tt>srcImageCopyProcessor</tt>
    * 
    * <p><tt>maxWidth</tt> holds width of image after rotation
    * 
    * <p><tt>ranges</tt> table holds first and last x position of image line (first and last pixel of
    * image on background after rotation)
    * 
    * <p><tt>srcImageCopyProcessor</tt> is rotated and shifted.
    * 
    * @throws DicException when input image is close to saturation e.g. has values of 65536-shift.
    *         This is due to applied algorithm of detection image pixels after rotation.
    * @see #reconstructionDicLid()
    */
   private void getRanges() throws DicException {
     double maxpixel; // minimal pixel value
     int lastpixel; // first and last pixel of image in line
     int firstpixel;
     // check condition for removing 0 value from image
     maxpixel = srcImageCopyProcessor.getMax();
     if (maxpixel > 65535 - shift) {
       LOGGER.error("Possible image clipping - check if image is saturated");
       throw new DicException(String.format(
               "Possible image clipping - input image has at leas one" + " pixel with value %d",
               65535 - shift));
     }
     // scale pixels by adding 1 - we remove any 0 value from source image
     srcImageCopyProcessor.add(shift);
     srcImageCopyProcessor.resetMinAndMax();
     // rotate image with extending it. borders have the same value as
     // background
     srcImageCopyProcessor = ipp.rotate(srcImageCopyProcessor, angle, true);
     isRotated = true; // current object was rotated
     // get references of covered object for optimisation purposes
     int newWidth = srcImageCopyProcessor.getWidth();
     int newHeight = srcImageCopyProcessor.getHeight();
     // set private fields - size of image after rotation
     maxWidth = newWidth;
     ranges = new int[newHeight][2];
     // get true pixels positions for every row
     for (int r = 0; r < newHeight; r++) {
       // to not process whole line, detect where starts and ends pixels of
       // image (reject background added during rotation)
       for (firstpixel = 0; firstpixel < newWidth
               && srcImageCopyProcessor.get(firstpixel, r) == 0; firstpixel++) {
         ;
       }
       for (lastpixel = newWidth - 1; lastpixel >= 0
               && srcImageCopyProcessor.get(lastpixel, r) == 0; lastpixel--) {
         ;
       }
       // if empty row, reject it all. reconstructionDicLid process pixels in ranges. such
       // borders reject processing at all
       if (firstpixel >= newWidth) {
         firstpixel = 1;
         lastpixel = 0;
       }
       ranges[r][0] = firstpixel;
       ranges[r][1] = lastpixel;
     }
   }
 
   /**
    * Recalculates tables on demand. Calculates new ranges for true pixels and new decay table.
    * 
    * @throws DicException Throws exception after generateRanges()
    */
   private void recalculate() throws DicException {
     // calculate preallocated decay data
     // generateRanges() must be called first as it initializes fields used
     // by generateDecay()
     getRanges();
     generateDeacy(decay, maxWidth);
   }
 
   /**
    * Reconstruct DIC image by LID method. This is main method used to reconstruct passed ip
    * object.
    * 
    * <p>The reconstruction algorithm assumes that input image bas-reliefs are oriented horizontally,
    * thus correct angle should be provided.
    * 
    * <p><b>Warning</b>
    * 
    * <p>Used optimization with detecting of image pixels based on their value may not be accurate
    * when input image will be 16-bit and it will contain saturated pixels
    * 
    * @return Return reconstruction of srcImage as 16-bit image
    */
   public ImageProcessor reconstructionDicLid() {
     double cumsumup;
     double cumsumdown;
     int c;
     int u;
     int d;
     int r; // loop indexes
     int linindex = 0; // output table linear index
     if (!isRotated) { // rotate if not rotated in getRanges
       srcImageCopyProcessor.add(shift); // we use different IP so shift must be added
       srcImageCopyProcessor = ipp.rotate(srcImageCopyProcessor, angle, true);
     }
     // dereferencing for optimization purposes
     int newWidth = srcImageCopyProcessor.getWidth();
     int newHeight = srcImageCopyProcessor.getHeight();
     // create array for storing results - 32bit float as imageprocessor
     ImageProcessor outputArrayProcessor = new FloatProcessor(newWidth, newHeight);
     float[] outputPixelArray = (float[]) outputArrayProcessor.getPixels();
 
     // do for every row - bas-relief is oriented horizontally
     for (r = 0; r < newHeight; r++) {
       // ranges[r][0] - first image pixel in line r
       // ranges[r][1] - last image pixel in line r
       linindex = linindex + ranges[r][0];
       // for every point apply KAM formula
       for (c = ranges[r][0]; c <= ranges[r][1]; c++) {
         // up
         cumsumup = 0;
         for (u = c; u >= ranges[r][0]; u--) {
           cumsumup += (srcImageCopyProcessor.get(u, r) - shift - is.mean) * decays[Math.abs(u - c)];
         }
         // down
         cumsumdown = 0; // cumulative sum from point r to the end of column
         for (d = c; d <= ranges[r][1]; d++) {
           cumsumdown +=
                   (srcImageCopyProcessor.get(d, r) - shift - is.mean) * decays[Math.abs(d - c)];
         }
         // integral
         outputPixelArray[linindex] = (float) (cumsumup - cumsumdown);
         linindex++;
       }
       linindex = linindex + newWidth - ranges[r][1] - 1;
     }
     // rotate back output processor
     outputArrayProcessor.setBackgroundValue(0.0);
     outputArrayProcessor.rotate(-angle);
     // crop it back to original size
     outputArrayProcessor =
             ipp.cropImageAfterRotation(outputArrayProcessor, srcIp.getWidth(), srcIp.getHeight());
 
     return outputArrayProcessor.convertToShort(true); // return reconstruction
   }
 
   /**
    * Generates decay table with exponential distances between pixels multiplied by decay
    * coefficient.
    * 
    * @param decay The value of decay coefficient
    * @param length Length of table, usually equals to longest processed line on image
    */
   private void generateDeacy(double decay, int length) {
     decays = new double[length];
     for (int i = 0; i < length; i++) {
       decays[i] = Math.exp(-decay * i);
     }
   }
 
 }