package com.github.celldynamics.quimp.plugin.randomwalk;
   
   import java.awt.Color;
   import java.io.File;
   import java.nio.file.Path;
   import java.nio.file.Paths;
   import java.security.InvalidParameterException;
   import java.util.ArrayList;
   import java.util.Arrays;
  import java.util.List;
  import java.util.concurrent.atomic.AtomicInteger;
  import java.util.stream.IntStream;
  
  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
  import org.apache.commons.math3.linear.ArrayRealVector;
  import org.apache.commons.math3.linear.MatrixDimensionMismatchException;
  import org.apache.commons.math3.linear.MatrixUtils;
  import org.apache.commons.math3.linear.RealMatrix;
  import org.apache.commons.math3.linear.RealMatrixChangingVisitor;
  import org.apache.commons.math3.stat.StatUtils;
  import org.slf4j.Logger;
  import org.slf4j.LoggerFactory;
  
  import com.github.celldynamics.quimp.QuimP;
  import com.github.celldynamics.quimp.utils.QuimPArrayUtils;
  
  import ij.IJ;
  import ij.ImagePlus;
  import ij.ImageStack;
  import ij.plugin.ImageCalculator;
  import ij.process.BinaryProcessor;
  import ij.process.ByteProcessor;
  import ij.process.ImageProcessor;
  
  /*
   * //!>
   * @startuml doc-files/RandomWalkSegmentation_1_UML.png
   * usecase IN as "Initialisation
   * --
   * Apply when main object is initialised.
   * Object is always connected to the image
   * that will be segmented."
   * 
   * usecase RUN as "Run
   * --
   * Related to running segmentation on
   * the object. The result is segmented
   * image"
   * 
   * usecase CON as "Convert
   * --
   * Related to various conversions
   * that are available as static
   * methods. They allow to change
   * representation of objects among
   * dataformats used by this object.
   * "
   * 
   * usecase PAR as "Parameters
   * --
   * Related to creating object
   * representing parameters
   * of the segmentation process.
   * "
   * Dev -> IN
   * Dev -> RUN
   * Dev --> CON
   * Dev --> PAR
   * 
   * RUN ..> IN : <<include>>
   * PAR -- IN
   * CON -- RUN
   * @enduml
   * //!<
   */
  /*
   * //!>
   * @startuml doc-files/RandomWalkSegmentation_2_UML.png
   * start
   * :Assign inputs
   * to internal fields;
   * note right
   * Inputs are **referenced**
   * end note
   * :Validate input
   * image;
   * :Convert image to
   * ""Matrix"";
   * note right
   * Or Matrix to image
   * depending on constructor
   * end note
   * stop
   * @enduml
   * //!<
   */
  /*
   * //!>
   * @startuml doc-files/RandomWalkSegmentation_3_UML.png
  * start
  * :Compute gradients;
  * note left
  * Gradients depend
  * on image structure
  * and they are tabelarised
  * for speed
  * end note
  * :solve problem;
  * -> //Result of segmentation, FG and BG probability maps//;
  * note right
  * Main procedure for
  * RW segmentation.
  * Return segmented image.
  * end note
  * if (next sweep?) then (yes)
  * :Run
  * intermediate
  * filter;
  * note left
  * Optional, run on
  * previous results
  * end note
  * :Prepare to
  * next sweep;
  * note left
  * Convert previous result
  * into seed, increment
  * ""currentSweep""
  * end note
  * :solve problem;
  * endif
  * :Compare FG and BG;
  * note left
  * Foreground and 
  * background probabilities
  * are compared to get 
  * unique result
  * end note
  * -> //Binary image//;
  * :Run
  * post-filtering;
  * note right 
  * Optional
  * end note
  * stop
  * @enduml
  * //!<
  */
 /*
  * //!>
  * @startuml doc-files/RandomWalkSegmentation_4_UML.png
  * actor User as user
  * participant RandomWalkSegmentation as rw
  * user -> rw : run(..)
  * activate rw
  * rw -> rw : precomputeGradients()
  *  activate rw
  *  rw --> rw : RealMatrix[]
  *  deactivate rw
  * rw -> rw : solver(..) 
  *  note right: Solver always return\nprobability map for\nforeground unless input\nFG was empty
  *  activate rw
  *  rw --> rw : Map<Seeds, RealMatrix> solved
  *  deactivate rw
  * alt next sweep?
  *  rw -> rw :  rollNextSweep(solved)
  *   activate rw
  *   rw --> rw : Map<Seeds, ImageProcessor> seedsNext
  *   deactivate rw
  *  rw -> rw : solver(seedsNext) 
  *   activate rw
  *   rw --> rw : Map<Seeds, RealMatrix> solved
  *   deactivate rw
  * end 
  * rw -> rw : compare(solved)
  *  note right: Compare FG and BG maps\nto get segmentation result
  *  activate rw
  *  rw --> rw : RealMatrix result
  *  deactivate rw
  * alt final filtering == yes
  *  rw -> Converter : RealMatrix
  *   activate Converter
  *   Converter --> rw : ImageProcessor
  *   deactivate Converter
  *  rw -> BinaryFilters : filter(ImageProcessor) 
  *   activate BinaryFilters
  *   BinaryFilters --> rw : filtered ImageProcessor
  *   deactivate BinaryFilters
  * else No final filtering
  *  rw -> Converter : RealMatrix
  *   activate Converter
  *   Converter --> rw : ImageProcessor
  *   deactivate Converter  
  * end 
  * note right: Except filtering, output can be\ncut by raw mask here, See run method.
  * rw --> user : ImageProcessor
  * deactivate rw
  * @enduml
  * //!<
  */
 /**
  * This is implementation of Matlab version of Random Walk segmentation algorithm.
  * 
  * <p>Available use cases are: <br>
  * <img src="doc-files/RandomWalkSegmentation_1_UML.png"/><br>
  * 
  * <h1>Parameters</h1>
  * In this Use Case user creates a set of segmentation parameters. They are hold in
  * {@link RandomWalkOptions}
  * class. This class also contains some pre- and post-processing settings used by
  * {@link RandomWalkSegmentation} object. Default constructor sets recommended values to all numeric
  * options skipping those related to pre- post-processing objects.
  * 
  * <h1>Initialisation</h1>
  * In this step the image that is supposed to be segmented is assigned to the object as well as
  * parameters of the process. There are two constructors available accepting <tt>ImageProcessor</tt>
  * and <tt>RelaMatrix</tt> image. These formats can be converted to each other using static
  * {@link QuimPArrayUtils#imageProcessor2RealMatrix(ImageProcessor)} and
  * {@link QuimPArrayUtils#realMatrix2ImageProcessor(RealMatrix)}.
  * <img src="doc-files/RandomWalkSegmentation_2_UML.png"/><br>
  * 
  * <h1>Convert</h1>
  * This Use Case contains preparatory static methods used for converting data into required format.
  * {@link RandomWalkSegmentation} uses {@link Seeds} structure to store information about foreground
  * (FG) and background (BG) seeds. Initially these seeds can be provided in different formats such
  * as:
  * <ol>
  * <li>RGB image where FG seeds are {@link Color#RED} and BG is {@link Color#GREEN}. In this mode
  * only these two seeds can be stored.
  * <li>Grayscale image, each object is labelled by different gray level. In this mode each object is
  * segmented separately, setting other objects to background.
  * <li>Binary image, works as grayscales.
  * <li>ROIs - rois selected on image.
  * </ol>
  * All above formats are converted to {@link Seeds} (the only accepted input by
  * {@link #run(Seeds)} by numerous methods from {@link SeedProcessor}. The
  * {@link Seeds} contains two keys important for segmentation, {@link SeedTypes#FOREGROUNDS} and
  * {@link SeedTypes#BACKGROUND}, {@link SeedTypes#FOREGROUNDS} is mandatory whereas
  * {@link SeedTypes#BACKGROUND} can be empty. Along one key seeds are stored as separate
  * <b>binary</b>
  * images that contain only pixels that belong to specified seed (e.g. red color for RGB or
  * specified grayscale for multi-cell segmentation). The optionl key {@link SeedTypes#ROUGHMASK}
  * contains one binary image that stands for initial estimation of object larger than the object. It
  * is used by local mean feature only ({@link RandomWalkOptions#useLocalMean} set <tt>true</tt>).
  * 
  * <h1>Run</h1>
  * In this Use Case the image assigned in <b>Initialisation</b> step is segmented and processed
  * according to parameters set in <b>Parameters</b>. The process is started by calling
  * {@link #run(Seeds)} method. For RGB seeds FG and BG are always defined but if the input is
  * converted from grayscale image, {@link SeedTypes#BACKGROUND} in {@link Seeds} is empty. The
  * method returns probability maps {@link ProbabilityMaps} that are organised similarly to
  * {@link Seeds}. Again for multi-cell segmentation returned {@link ProbabilityMaps} may not contain
  * {@link SeedTypes#BACKGROUND} key. The algorithm assumes that background is on probability 0 and
  * it is {@link Color#BLACK}. Probabilities are compared in {@link #compare(ProbabilityMaps)}
  * method that returns 2d matrix being an output image.
  * 
  * <p>Here is brief look for {@link #run(Seeds)} method activity:
  * <img src="doc-files/RandomWalkSegmentation_3_UML.png"/><br>
  * 
  * <p>Sequence diagram of {@link #run(Seeds)} giving the order of called methods:
  * <img src="doc-files/RandomWalkSegmentation_4_UML.png"/><br>
  * 
  * <p>See: src/test/Resources-static/Matlab/rw_laplace4.m
  * 
  * @author p.baniukiewicz
  */
 public class RandomWalkSegmentation {
 
   /**
    * How often we will compute relative error.
    */
   final int relErrStep = 20;
 
   /**
    * Keep information about current sweep that {@link #solver(Seeds, RealMatrix[])} is doing.
    * 
    * <p>Affect indexing parameters arrays that keep different params for different sweeps. Used also
    * in {@link #solver(Seeds, RealMatrix[]) for computing actual number of iterations (second sweep
    * uses half of defined)}
    * 
    * @see RandomWalkOptions
    */
   private int currentSweep = 0;
 
   /**
    * Squared maximal theoretical intensity value.
    * 
    * <p>For 8-bit images it is 255^2, for 16-bit 65535^2. Set in constructor.
    */
   private int maxTheoreticalIntSqr;
 
   /**
    * Reasons of stopping diffusion process.
    * 
    * @author p.baniukiewicz
    *
    */
   private enum StoppedBy {
     /**
      * Maximum number of iterations reached.
      */
     ITERATIONS(0),
     /**
      * Found NaN in solution.
      */
     NANS(1),
     /**
      * Found Inf in solution.
      */
     INFS(2),
     /**
      * Relative error smaller than limit.
      */
     RELERR(3);
 
     private final int value;
 
     private StoppedBy(int value) {
       this.value = value;
     }
 
     public int getValue() {
       return value;
     }
   }
 
   /**
    * Define foreground and background indexes enums with numerical indexes.
    * 
    * @author p.baniukiewicz
    *
    */
   public enum SeedTypes {
     /**
      * Denote foreground related data. Usually on index 0.
      * 
      * <p>FG data can be grayscale images up to 8-bits.
      */
     FOREGROUNDS(0),
     /**
      * Denote background related data. Usually on index 1.
      * 
      * <p>BG data can be grayscale image up to 8-bits.
      */
     BACKGROUND(1),
     /**
      * Rough mask used for computing local mean. Used only if {@link RandomWalkOptions#useLocalMean}
      * is true. This mask should be always binary.
      * 
      * @see RandomWalkSegmentation#solver(Seeds, RealMatrix[])
      * @see RandomWalkSegmentation#getMeanSeedLocal(ImageProcessor, int)
      */
     ROUGHMASK(2);
 
     /** The index. */
     private final int index;
 
     /**
      * Instantiates a new seed types.
      *
      * @param index the index
      */
     private SeedTypes(int index) {
       this.index = index;
     }
 
     /**
      * Convert enum to int. Used when addressing arrays.
      * 
      * @return index assigned to enum value.
      */
     public int getIndex() {
       return index;
     }
   }
 
   /**
    * The Constant LOGGER.
    */
   static final Logger LOGGER = LoggerFactory.getLogger(RandomWalkSegmentation.class.getName());
 
   /**
    * Direction of circshift coded as in Matlab.
    */
   public static final int RIGHT = -10;
   /**
    * Direction of circshift coded as in Matlab.
    */
   public static final int LEFT = 10;
   /**
    * Direction of circshift coded as in Matlab.
    */
   public static final int TOP = -01;
   /**
    * Direction of circshift coded as in Matlab.
    */
   public static final int BOTTOM = 01;
 
   /**
    * Image to process in 8bit greyscale. Converted to RealMatrix
    * 
    * @see #ip
    */
   private RealMatrix image;
   /**
    * Original image to process.
    * 
    * @see #image
    */
   private ImageProcessor ip;
   /**
    * User provided parameters.
    */
   private RandomWalkOptions params;
   /**
    * Probability map obtained in {@link #run(Seeds)}.
    */
   private ProbabilityMaps solved = null;
 
   /**
    * Construct segmentation object from ImageProcessor.
    * 
    * @param ip image to segment
    * @param params parameters
    * @throws RandomWalkException on wrong image format
    */
   public RandomWalkSegmentation(ImageProcessor ip, RandomWalkOptions params)
           throws RandomWalkException {
     if (ip.getBitDepth() != 8 && ip.getBitDepth() != 16) {
       throw new RandomWalkException("Only 8-bit or 16-bit images are supported");
     }
     this.ip = ip;
     this.image = QuimPArrayUtils.imageProcessor2RealMatrix(ip);
     this.params = params;
     setMaxTheoreticalIntSqr(ip);
   }
 
   /**
    * Construct segmentation object from 2D RealMatrix representing image.
    * 
    * <p>It is assumed that this input image is 8-bit. Use
    * {@link #RandomWalkSegmentation(ImageProcessor, RandomWalkOptions)} for support 8 and 16-bit
    * images.
    * Passing wrong image can have effect to results as {@link #solver(Seeds, RealMatrix[])}
    * normalises
    * image intensities to maximal theoretical intensity. See
    * {@link #setMaxTheoreticalIntSqr(ImageProcessor)} and {@link #solver(Seeds, RealMatrix[])}
    * 
    * @param image image to segment
    * @param params parameters
    */
   public RandomWalkSegmentation(RealMatrix image, RandomWalkOptions params) {
     this.image = image;
     this.ip = QuimPArrayUtils.realMatrix2ImageProcessor(image);
     this.params = params;
     setMaxTheoreticalIntSqr(ip);
   }
 
   /**
    * Main runner, does segmentation.
    * 
    * <p>Requires defined FOREGROUNDS seeds and one BACKGROUND. If background is not given it assumes
    * during weighting ({@link #compare(ProbabilityMaps)}) background probability map with
    * probability 0. IT should wor even if there is only one seed FG and no BG because most arrays
    * are 0-filled by default.
    * 
    * @param seeds Seed arrays from {@link SeedProcessor}
    * @return Segmented image as ByteProcessor or null if segmentation failed due to e.g. empty seeds
    * @throws RandomWalkException On wrong seeds
    */
   public ImageProcessor run(Seeds seeds) throws RandomWalkException {
     LOGGER.debug("Running with options: " + params.toString());
     if (seeds.get(SeedTypes.FOREGROUNDS) == null) {
       return null; // no FG maps - no segmentation
     }
     RealMatrix[] precomputed = precomputeGradients(); // precompute gradients
     solved = solver(seeds, precomputed);
     if (params.intermediateFilter != null && params.gamma[1] != 0) { // do second sweep
       LOGGER.debug("Running next sweep: " + params.intermediateFilter.getClass().getName());
       Seeds seedsNext = rollNextSweep(solved);
       if (seeds.get(SeedTypes.ROUGHMASK) != null) {
         seedsNext.put(SeedTypes.ROUGHMASK, seeds.get(SeedTypes.ROUGHMASK, 1));
       }
       solved = solver(seedsNext, precomputed);
     }
     RealMatrix result = compare(solved); // result as matrix
     ImageProcessor resultim =
             QuimPArrayUtils.realMatrix2ImageProcessor(result).convertToByteProcessor(true);
     // cut mask - cut segmentation result by initial ROUGHMASK if present
     if (params.maskLimit == true && seeds.get(SeedTypes.ROUGHMASK) != null) {
       ImageCalculator ic = new ImageCalculator();
       ImagePlus retc = ic.run("and create",
               new ImagePlus("", new BinaryProcessor(
                       (ByteProcessor) seeds.get(SeedTypes.ROUGHMASK, 1).convertToByte(true))),
               new ImagePlus("", new BinaryProcessor((ByteProcessor) resultim)));
       resultim = retc.getProcessor();
     }
     // do final filtering
     if (params.finalFilter != null) {
       return params.finalFilter.filter(resultim.convertToByteProcessor(true));
     } else {
       return resultim.convertToByteProcessor(true);
     }
 
   }
 
   /**
    * Prepare seeds from results of previous solver.
    * 
    * @param solved results from {@link #solver(Seeds, RealMatrix[])}
    * @return new seed taken from previous solution.
    * @throws RandomWalkException on unsupported image or empty seed list after decoding
    */
   private Seeds rollNextSweep(ProbabilityMaps solved) throws RandomWalkException {
     final double weight = 1e20;
     Seedsamics/quimp/plugin/randomwalk/Seeds.html#Seeds">Seeds ret = new Seeds(2);
     RealMatrix solvedWeighted;
     // make copy of input results to weight them
     for (int m = 0; m < solved.get(SeedTypes.FOREGROUNDS).size(); m++) {
       solvedWeighted = solved.get(SeedTypes.FOREGROUNDS).get(m).copy();
 
       // seed_fg = FGl>1e20*BGl
       solvedWeighted.walkInOptimizedOrder(
               new MatrixCompareWeighted(solved.get(SeedTypes.BACKGROUND).get(0), weight));
 
       // convert weighted results to images
       ImageProcessor fg1 =
               QuimPArrayUtils.realMatrix2ImageProcessor(solvedWeighted).convertToByte(true);
 
       if (QuimP.SUPER_DEBUG) { // save intermediate results
         LOGGER.debug("Saving intermediate results");
         String tmpdir = System.getProperty("java.io.tmpdir") + File.separator;
         IJ.saveAsTiff(new ImagePlus("", fg1), tmpdir + "fg1_" + m + "_QuimP.tif");
       }
       // filter them (if we are here intermediate filter can't be null)
       fg1 = params.intermediateFilter.filter(fg1);
       ret.put(SeedTypes.FOREGROUNDS, fg1);
     }
     // increment sweep pointer to point correct parameters (if they are different for next sweep)
     currentSweep++;
     solvedWeighted = solved.get(SeedTypes.BACKGROUND).get(0).copy();
     // flatten foregrounds to weight background
     double[][] fl = flatten(solved, SeedTypes.FOREGROUNDS);
     // seed_bg = BGl>1e20*FGl;
     solvedWeighted.walkInOptimizedOrder(
             new MatrixCompareWeighted(MatrixUtils.createRealMatrix(fl), weight));
     ImageProcessor bg1 =
             QuimPArrayUtils.realMatrix2ImageProcessor(solvedWeighted).convertToByte(true);
     if (QuimP.SUPER_DEBUG) { // save intermediate results
       String tmpdir = System.getProperty("java.io.tmpdir") + File.separator;
       IJ.saveAsTiff(new ImagePlus("", bg1), tmpdir + "bg1_\"+m+\"_QuimP.tif");
     }
     bg1.invert();
     bg1 = params.intermediateFilter.filter(bg1);
     bg1.invert();
     ret.put(SeedTypes.BACKGROUND, bg1);
 
     return ret;
   }
 
   /**
    * Compare probabilities from and create segmentation matrix depending on winner.
    * 
    * @param solved foreground and background seeds. If there is no background seed this procedure
    *        assume probability map of size of foreground filled with 0s
    * @return matrix of size of input image where objects are labelled with constitutive numbers
    *         starting from 1
    */
   RealMatrix compare(ProbabilityMaps solved) {
     double backgroundColorValue = 0.0; // value for fill background
     double[][][] fgmaps3d = solved.convertTo3dMatrix(SeedTypes.FOREGROUNDS);
     double[][][] bgmaps3d = solved.convertTo3dMatrix(SeedTypes.BACKGROUND);
 
     if (fgmaps3d == null) {
       return null;
     }
     int rows = fgmaps3d[0].length;
     int cols = fgmaps3d[0][0].length;
     // assume probability of 0 for BG. If there are only FG objects and probability of being FG1 is
     // 0 as well as for FG2, this will impose background
     if (bgmaps3d == null) {
       bgmaps3d = new double[1][rows][cols];
     }
     if (rows == 0 || cols == 0 || rows != bgmaps3d[0].length || cols != bgmaps3d[0][0].length) {
       return null;
     }
 
     RealMatrix ret = MatrixUtils.createRealMatrix(rows, cols);
 
     int[][] fgMaxMap = flattenInd(fgmaps3d);
     int[][] bgMaxMap = flattenInd(bgmaps3d);
 
     for (int r = 0; r < rows; r++) {
       for (int c = 0; c < cols; c++) {
         int fi = fgMaxMap[r][c]; // z-index of max element in FG
         int bi = bgMaxMap[r][c]; // z-index of max element in bg
         if (fgmaps3d[fi][r][c] > bgmaps3d[bi][r][c]) { // FG>BG - object
           ret.setEntry(r, c, fi + 1); // set object coded at index z+1 (1-based)
         } else {
           ret.setEntry(r, c, backgroundColorValue); // background
         }
       }
     }
     return ret;
 
   }
 
   /**
    * Compare values along z dimension for each x,y and return index of max value.
    * 
    * @param in matrix to flatten
    * @return 2d matrix of indexes (0-based) of maximal values for each x,y along z
    */
   private int[][] flattenInd(double[][][] in) {
     int rows = in[0].length;
     int cols = in[0][0].length;
     int[][] ret = new int[rows][cols];
 
     if (in.length == 1) {
       return ret; // just return 0 indexes
     }
     for (int r = 0; r < rows; r++) {
       for (int c = 0; c < cols; c++) {
         double max = in[0][r][c];
         for (int z = 0; z < in.length; z++) {
           if (in[z][r][c] > max) {
             max = in[z][r][c];
             ret[r][c] = z;
           }
         }
       }
     }
     return ret;
   }
 
   /**
    * Produce 2d matrix of max values taken along z.
    * 
    * @param maps maps to flatten to flatten
    * @param key seed key
    * 
    * @return 2d matrix of max values (0-based) of maximal values for each x,y along z
    */
   double[][] flatten(ProbabilityMaps maps, SeedTypes key) {
     double[][][] in = maps.convertTo3dMatrix(key);
     if (in == null) {
       return null;
     }
     int rows = in[0].length;
     int cols = in[0][0].length;
     double[][] ret = new double[rows][cols];
 
     if (in.length == 1) {
       return in[0]; // just return input
     }
     for (int r = 0; r < rows; r++) {
       for (int c = 0; c < cols; c++) {
         double max = in[0][r][c];
         for (int z = 0; z < in.length; z++) {
           if (in[z][r][c] > max) {
             max = in[z][r][c];
           }
           ret[r][c] = max;
         }
       }
     }
     return ret;
   }
 
   /**
    * Shift image in given direction with step of one and circular boundary conditions.
    *
    * @param input Image to be shifted
    * @param direction Shift direction. This method is adjusted to be compatible with Matlab code
    *        and to keep Matlab naming (rw_laplace4.m) thus the shift direction names are not
    *        adequate to shift direction.
    * @return Copy of input shifted by one pixel in direction.
    */
   RealMatrix circshift(RealMatrix input, int direction) {
     // Routine cuts part of matrix that does not change, paste it into new matrix shifting in given
     // direction and then add missing looped row/column (last or first depending on shift
     // direction).
     double[][] sub; // part of matrix that does no change but is shifted
     int rows = input.getRowDimension(); // cache sizes
     int cols = input.getColumnDimension();
     Array2DRowRealMatrix out = new Array2DRowRealMatrix(rows, cols); // output matrix, shifted
     switch (direction) {
       case BOTTOM:
         // rotated right - last column become first
         // cut submatrix from first column to before last
         sub = new double[rows][cols - 1];
         input.copySubMatrix(0, rows - 1, 0, cols - 2, sub); // cols-2 - cols is size not last ind
         // create new matrix - paste submatrix but shifted right
         out.setSubMatrix(sub, 0, 1);
         // copy last column to first
         out.setColumnVector(0, input.getColumnVector(cols - 1));
         break;
       case TOP:
         // rotated left - first column become last
         // cut submatrix from second column to last
         sub = new double[rows][cols - 1];
         input.copySubMatrix(0, rows - 1, 1, cols - 1, sub);
         // create new matrix - paste submatrix but shifted right
         out.setSubMatrix(sub, 0, 0);
         // copy first column to last
         out.setColumnVector(cols - 1, input.getColumnVector(0));
         break;
       case RIGHT:
         // rotated top - first row become last
         // cut submatrix from second row to last
         sub = new double[rows - 1][cols];
         input.copySubMatrix(1, rows - 1, 0, cols - 1, sub);
         // create new matrix - paste submatrix but shifted up
         out.setSubMatrix(sub, 0, 0);
         // copy first row to last
         out.setRowVector(rows - 1, input.getRowVector(0));
         break;
       case LEFT:
         // rotated bottom - last row become first
         // cut submatrix from first row to before last
         sub = new double[rows - 1][cols];
         input.copySubMatrix(0, rows - 2, 0, cols - 1, sub);
         // create new matrix - paste submatrix but shifted up
         out.setSubMatrix(sub, 1, 0);
         // copy last row to first
         out.setRowVector(0, input.getRowVector(rows - 1));
         break;
       default:
         throw new IllegalArgumentException("circshift: Unknown direction");
     }
     return out;
   }
 
   /**
    * Compute (a-b)^2.
    * 
    * @param a left operand
    * @param b right operand
    * @return Image (a-b).^2
    */
   RealMatrix getSqrdDiffIntensity(RealMatrix a, RealMatrix b) {
     RealMatrix s = a.subtract(b);
     s.walkInOptimizedOrder(new MatrixElementPower());
     return s;
   }
 
   /**
    * Pre-compute normalised gradient matrices.
    * 
    * <p>This computes squared intensity differences.These correspond to horizontal and vertical
    * gradients assuming distance of one between pixels.
    * 
    * @return Array of precomputed data in the following order: [0] - gRight2 [1] - gTop2 [2] -
    *         gLeft2 [3] - gBottom2
    */
   private RealMatrix[] precomputeGradients() {
     // setup shifted images assuming shift of one pixel and periodic boundary conditions.
     RealMatrix right = circshift(image, RIGHT);
     RealMatrix top = circshift(image, TOP);
     // compute squared intensity differences (a-b).^2
     RealMatrix gradRight2 = getSqrdDiffIntensity(image, right);
     RealMatrix gradTop2 = getSqrdDiffIntensity(image, top);
     // compute maximum of horizontal and vertical intensity gradients (for normalization)
     double maxGright2 = QuimPArrayUtils.getMax(gradRight2);
     LOGGER.debug("maxGright2 " + maxGright2);
     double maxGtop2 = QuimPArrayUtils.getMax(gradTop2);
     LOGGER.debug("maxGtop2 " + maxGtop2);
     // find maximum of horizontal and vertical maxima
     double maxGrad2 = maxGright2 > maxGtop2 ? maxGright2 : maxGtop2;
     LOGGER.debug("maxGrad2max " + maxGrad2);
     if (maxGrad2 == 0) {
       maxGrad2 = 1.0;
     }
     // Normalize squared gradients to maxGrad
     gradRight2.walkInOptimizedOrder(new MatrixElementDivide(maxGrad2));
     gradTop2.walkInOptimizedOrder(new MatrixElementDivide(maxGrad2));
     // assign outputs
     RealMatrix[] out = new RealMatrix[4];
     // get remaining directions - use already normalized squared gradients
     RealMatrix gradLeft2 = circshift(gradRight2, LEFT);
     out[2] = gradLeft2;
     RealMatrix gradBottom2 = circshift(gradTop2, BOTTOM);
     out[3] = gradBottom2;
     out[0] = gradRight2;
     out[1] = gradTop2;
 
     return out;
   }
 
   /**
    * Compute mean value from image only seeded pixels.
    * 
    * @param seeds coordinates of points used to calculate their mean intensity
    * @return mean value for points
    * @see Seeds#convertToList(SeedTypes)
    */
   protected double getMeanSeedGlobal(List<Point> seeds) {
     return StatUtils.mean(getValues(image, seeds).getDataRef());
   }
 
   /**
    * Calculate local mean intensity for input image but only for areas masked by specified mask.
    * 
    * <p>The mean over segmented image intensity is evaluated within square window of configurable
    * size and only for those pixels that are masked by binary mask given to this method. Mean is
    * calculated respectfully to the number of masked pixels within window. Window is moved over the
    * whole image.
    * 
    * <p>This method works similarly to the convolution with the difference that the kernel is
    * normalised for each position of the window to the number of masked pixels (within it). If for
    * any position of the window there are no masked pixels inside, value 0.0 is set as result.
    * 
    * @param mask Binary mask of segmented image. Mask must contain only pixels with intensity 0 or
    *        255 (according to definition of binary image in IJ)
    * @param localMeanMaskSize Odd size of kernel
    * @return Averaged image. Average is computed for every location of the kernel for its center
    *         utlising only masked pixels. For not masked pixels the mean is set to 0.0
    */
   protected RealMatrix getMeanSeedLocal(ImageProcessor mask, int localMeanMaskSize) {
     // IJ convolution is utilised. Computations are done twice, first time only on mask to get the
     // number of masked pixels covered by it for each image pixel
     // second time for segmented image to get sum of intensities within the kernel.
     // in both cases kernels are normalised (default behaviour of ImageProcessor.convolve) but then
     // "denormalised" by multiplying by their size. Kernels are 1-filled.
     if (localMeanMaskSize % 2 == 0) {
       throw new IllegalArgumentException("Kernel sie must be odd");
     }
     // must be binary as we later assume that mx intensity is 255 (like in IJ binary images)
     if (!mask.isBinary()) {
       throw new IllegalArgumentException("Mask must be binary");
     }
     if (mask.getWidth() != ip.getWidth() || mask.getHeight() != ip.getHeight()) {
       throw new IllegalArgumentException("Mask must have size of processed image");
     }
     ImageProcessor maskc = mask.duplicate(); // local copy of mask to not modify it
     maskc.subtract(254); // if it is IJ binary it contains [0,255], scale to 0 - 1.0
     // make mask copy for getting number of pixels in mask for each position of kernel
     ImageProcessor numofpix = maskc.duplicate();
 
     // generate 1-filled kernel of given size
     float[] kernel = new float[localMeanMaskSize * localMeanMaskSize];
     Arrays.fill(kernel, 1.0f);
     // get number of mask pixesl for ach position of kernel
     numofpix = maskc.convertToFloat(); // must convert here
     numofpix.convolve(kernel, localMeanMaskSize, localMeanMaskSize);
     numofpix.multiply(kernel.length); // convolution normalises kernel - revert it 1)
     // cut to mask - mulitply by 1.0 mask - reject pixels that are not masked
     numofpix = new ImageCalculator()
             .run("mul create", new ImagePlus("", numofpix), new ImagePlus("", maskc))
             .getProcessor();
 
     // get sum of intensities within the kernel for segmented image
     // cut image to input mask - set all other pixels to 0. Must be done here to not count nonmasked
     // pixels
     ImageProcessor cutImage = new ImageCalculator()
             .run("mul create float", new ImagePlus("", ip), new ImagePlus("", maskc))
             .getProcessor();
     // convolve cut image
     cutImage.setCalibrationTable(null);
     cutImage.convolve(kernel, localMeanMaskSize, localMeanMaskSize);
     cutImage.multiply(kernel.length); // denormalise result
     // deal with edges of image after convolution - remove them
     cutImage = new ImageCalculator()
             .run("mul create float", new ImagePlus("", cutImage), new ImagePlus("", maskc))
             .getProcessor();
     // rounding floats to integer values - convolution works for float images only
     float[] cutImageRaw = (float[]) cutImage.getPixels();
     for (int i = 0; i < cutImage.getPixelCount(); i++) {
       cutImageRaw[i] = Math.round(cutImageRaw[i]);
     }
     float[] numofpixRaw = (float[]) numofpix.getPixels();
     for (int i = 0; i < numofpix.getPixelCount(); i++) {
       numofpixRaw[i] = Math.round(numofpixRaw[i]);
     }
     // proper mean - use only pixels inside mask. we use cutImageRaw that contains sums of
     // intensities for each location of kernel and numofpixRaw that contains number of masked pixels
     RealMatrix meanseedFg = new Array2DRowRealMatrix(cutImage.getHeight(), cutImage.getWidth());
     int count = 0;
     for (int r = 0; r < cutImage.getHeight(); r++) {
       for (int c = 0; c < cutImage.getWidth(); c++) {
         if (numofpix.getPixelValue(c, r) != 0) { // skip 0 pixels to not divide by 0
           meanseedFg.setEntry(r, c, ((double) cutImageRaw[count] / numofpixRaw[count]));
         } else {
           meanseedFg.setEntry(r, c, 0.0); // if no pixels in mask for current r,c - set mean to 0
         }
         count++;
       }
     }
 
     if (QuimP.SUPER_DEBUG) { // save intermediate results
       LOGGER.debug("Saving intermediate results");
       String tmpdir = System.getProperty("java.io.tmpdir") + File.separator;
       IJ.saveAsTiff(new ImagePlus("", numofpix), tmpdir + "maskc_QuimP.tif");
       IJ.saveAsTiff(new ImagePlus("", cutImage), tmpdir + "imagecc_QuimP.tif");
       IJ.saveAsTiff(new ImagePlus("", QuimPArrayUtils.realMatrix2ImageProcessor(meanseedFg)),
               tmpdir + "meanseedFg_QuimP.tif");
       LOGGER.trace("meanseedFg M[183;289] " + meanseedFg.getEntry(183 - 1, 289 - 1));
       LOGGER.trace("meanseedFg M[242;392] " + meanseedFg.getEntry(242 - 1, 392 - 1));
       LOGGER.trace("cutImage M[192;279] " + cutImage.getPixelValue(279 - 1, 192 - 1));
     }
     return meanseedFg;
   }
 
   /*
    * //!>
    * @startuml doc-files/RandomWalkSegmentation_5_UML.png
    * start
    * :Convert seeds to coordinates;
    * :Combine FG and BG maps together;
    * repeat :solve for object **i**
    * :Set object **i** as FG;
    * :Set all other objects as BG;
    * if (local mean?) then (yes)
    *  :compute **local** mean for **FG**;
    *  note left
    *  Result is an image
    *  end note
    *  :compute **global** mean for **BG**;
    *  note left
    *  Result is a number
    *  end note
    *  :subtract FG and BG means\nfrom image;
    * else (no)
    *  :compute **global** mean for **FG**;
    *  :compute **global** mean for **BG**;
    *  :subtract FG and BG means\nfrom image;
    * endif
    * -> Here we have Image-meanseed;
    * :compute normalised\n""(Image-meanseed).^2"";
    * note left
    * Square ""Image-meanseed"" and
    * normalise to maximal theoretical
    * intensity value
    * end note
    * :compute weights;
    * note right
    * Based on intensity of
    * pixels in stencil and 
    * intensity gradient
    * end note
    * :compute average of weights;
    * note left
    * Separately for horizontal
    * and vertical differences.
    * Used later for equalising
    * grid
    * end note
    * :compute diffusion constant;
    * note left
    * To obey stability
    * criterion
    * end note
    * :compute averaged weights;
    * note right
    * Equalising differentiation
    * grid.
    * end note
    * :compute FG;
    * note left
    * Solving Laplace equation
    * by Euler method (loop, number of iter
    * depends on sweep number)
    * end note
    * :Store map;
    * note right
    * Generally algorithm assumes that 
    * last map is BG and can detect it 
    * and store under proper key in Seeds
    * end note
    * repeat while (more objects?)
    * stop
    * @enduml
    * //!<
    */
   /**
    * Run Random Walk segmentation.
    * 
    * <p>The activity diagram for solver is as follows:
    * <img src="doc-files/RandomWalkSegmentation_5_UML.png"/><br>
    * 
    * <p>Note 1: that solver treats FG and BG seeds like equal objects. There is no difference
    * between foreground and background seeds.
    * 
    * <p>Note 2:number of iterations for BG object is limited to maximum number of iterations that
    * occurred for FG objects.
    * 
    * @param seeds seed array returned from {@link SeedProcessor}
    * @param gradients pre-computed gradients returned from {@link #precomputeGradients()}
    * @return Computed probabilities for background and foreground. Returned structure can be also
    *         empty if there are not FG seeds provided on input (no FG map in {@link Seeds}) or may
    *         not contain BG map.
    */
   protected ProbabilityMaps solver(Seeds seeds, RealMatrix[] gradients) {
     ImageStack debugPm = null;
     RealMatrix diffIfg = null; // normalised squared differences to mean seed intensities for FG
     // store number of iterations performed for each object. These numbers are used for stopping
     // iterations earlier for background object
     ArrayList<Integer> iterations = new ArrayList<>();
 
     ProbabilityMapsp/plugin/randomwalk/ProbabilityMaps.html#ProbabilityMaps">ProbabilityMaps ret = new ProbabilityMaps(); // keep output probab map for each object
 
     if (seeds.get(SeedTypes.FOREGROUNDS) == null) { // if no FG maps (e.g. all disappeared)
       return ret;
     }
     // seed images as list of points
     List<List<Point>> seedsPointsFg = seeds.convertToList(SeedTypes.FOREGROUNDS);
     // user selected background (it is solved like other objects but e.g. local mean does not apply
     // for it so we need to know what was selected by user as background). There is possible
     // that there will be no BCK key, if input is given as GraySclale image
     List<List<Point>> userBckPoints = seeds.convertToList(SeedTypes.BACKGROUND);
     // add background at the end - it will be solved as regular object
     seedsPointsFg.addAll(userBckPoints);
     // background points used in conjunction with current foreground. Background points are all
     // other points which are not current foreground (e.g. other cells + user background, or all
     // cells if we solve for user background)
     List<List<Point>> seedsPointsBg = new ArrayList<>();
 
     // solve problem for each label in FOREGROUNDS, other labels (if any) are merged with BACKGROUND
     for (int cell = 0; cell < seedsPointsFg.size(); cell++) { // over each seed label
       // some maps for FOREGROUNDS key can be empty, so lists will be too. Note that decodeSeeds
       // throw exception when all maps for specified key are empty. Other situations are allowed.
       if (seedsPointsFg.get(cell).isEmpty()) {
         continue;
       }
       // make copy of objects seeds - need of removing current one and integrate remaining with bck
       ArrayList<List<Point>> tmpSeedsPointsFg = new ArrayList<List<Point>>(seedsPointsFg);
       tmpSeedsPointsFg.remove(cell); // remove current object seed
       // clear Bg, it always contains all objects except current (seedsPointsFg[cell])
       seedsPointsBg.clear();
       seedsPointsBg.addAll(tmpSeedsPointsFg); // add all remaining foregrounds to current bck
 
       // decide whether to use local mean or global mean. Local mean is computed within square mask
       // of configurable size whereas the global mean is a mean intensity of all seeded pixels.
       // Local mean evaluated only for FG objects.
       if (params.useLocalMean && seeds.get(SeedTypes.ROUGHMASK) != null
               && !userBckPoints.contains(seedsPointsFg.get(cell))) { // skip BG object
         RealMatrix localMeanFg =
                 getMeanSeedLocal(seeds.get(SeedTypes.ROUGHMASK, 1), params.localMeanMaskSize);
         diffIfg = image.subtract(localMeanFg);
       } else { // global for whole seeds
         // compute intensity means for image points labelled by seeds
         double meanseed = getMeanSeedGlobal(seedsPointsFg.get(cell));
         LOGGER.debug("meanseed_fg=" + meanseed);
         // compute normalised squared differences to mean seed intensities (Image-meanseed).^2
         diffIfg = image.scalarAdd(-meanseed);
       }
       // normalize (Image-meanseed).^2 to maximal (theoretical) value which is 255^2 for 8-bit
       // images. Have it as private field as we support 16 images as well
       diffIfg.walkInOptimizedOrder(new MatrixElementPowerDiv(maxTheoreticalIntSqr));
       LOGGER.trace("fseeds size: " + seedsPointsFg.get(cell).size());
       LOGGER.trace("bseeds size: " + seedsPointsBg.stream().mapToInt(p -> p.size()).sum());
 
       // compute weights for diffusion in all four directions, dependent on local gradients and
       // differences to mean intensities of seeds, for FG maps
       Array2DRowRealMatrix wrfg = (Array2DRowRealMatrix) computeweights(diffIfg, gradients[0]);
       Array2DRowRealMatrix wlfg = (Array2DRowRealMatrix) computeweights(diffIfg, gradients[2]);
       Array2DRowRealMatrix wtfg = (Array2DRowRealMatrix) computeweights(diffIfg, gradients[1]);
       Array2DRowRealMatrix wbfg = (Array2DRowRealMatrix) computeweights(diffIfg, gradients[3]);
 
       // compute averaged weights, left/right and top/bottom used when computing second spatial
       // derivative from first one, avgwx_fg = 0.5*(wl_fg+wr_fg) - for FG
       RealMatrix avgwxfg = wlfg.add(wrfg); // wl_fg+wr_fg - (left+right)
       avgwxfg.walkInOptimizedOrder(new MatrixElementMultiply(0.5)); // 0.5*(wl_fg+wr_fg)
       RealMatrix avgwyfg = wtfg.add(wbfg); // wt_fg+wb_fg - (top+bottom)
       avgwyfg.walkInOptimizedOrder(new MatrixElementMultiply(0.5)); // 0.5*(wt_fg+wb_fg)
 
       // Compute diffusion coefficient that will obey stability criterion
       double diffusion = getDiffusionConst(wrfg, wlfg, wtfg, wbfg, avgwxfg, avgwyfg);
       LOGGER.debug("D=" + diffusion);
 
       // get average "distance" between weights multiplying w = w.*avgw, this is only for
       // optimisation purposes.
       // does not apply for FG if we use local mean, applied for BG always (better results)
       if (params.useLocalMean == false || userBckPoints.contains(seedsPointsFg.get(cell))) {
         wrfg.walkInOptimizedOrder(new MatrixDotProduct(avgwxfg));
         wlfg.walkInOptimizedOrder(new MatrixDotProduct(avgwxfg));
         wtfg.walkInOptimizedOrder(new MatrixDotProduct(avgwyfg));
         wbfg.walkInOptimizedOrder(new MatrixDotProduct(avgwyfg));
       }
 
       // initialize FG probability maps, they are outputs from this routine
       Array2DRowRealMatrix fg =
               new Array2DRowRealMatrix(image.getRowDimension(), image.getColumnDimension());
       // this temporary array will keep solution from n-1 iteration used for computing rel error
       double[][] tmpFglast2d = new double[image.getRowDimension()][image.getColumnDimension()];
 
       // dereferencing arrays, all computations are done on underlying [][] arrays but not on
       // RelaMatrix objects, optimisation again
       double[][] wrfg2d = wrfg.getDataRef(); // weight to right for FG
       double[][] wlfg2d = wlfg.getDataRef(); // weight to left for FG
       double[][] wtfg2d = wtfg.getDataRef(); // weight to top for FG
       double[][] wbfg2d = wbfg.getDataRef(); // weight to bottom for FG
 
       double[][] fg2d = fg.getDataRef(); // reference to FG probabilities
 
       StoppedBy stoppedReason = StoppedBy.ITERATIONS; // default assumption
       int i; // iteration counter
       // compute correct number of iterations. Second sweep uses 0.5*user
       int iter;
       // use less iterations when diffuse background. Background needs much more iterations to reach
       // specified relError and after weighting it dominates leaving only original object seed as
       // segmented object. Here we stop segmenting background after certain number of iterations but
       // not relErr. This is how we have it solved in MAtlab
       if (true && (userBckPoints.contains(seedsPointsFg.get(cell)) && iterations.size() > 0)) {
         // just use average of iters for BCK
         iter = iterations.stream().mapToInt(Integer::intValue).max().getAsInt();
         // FIXME This can be disabled, then BCK will need more iterations but sometimes results are
         // better
         // potential pitfall is if user mark BG far from cell, then small number of iters is not
         // enough to flood whole background (but it will work because during comparison bck is on 0
         // and FG segmentation rather does not leave object
         iter /= (currentSweep + 1);
       } else { // object - use specified number of iters
         iter = params.iter / (currentSweep + 1);
       }
       // main loop here we simulate diffusion process in time
       outerloop: for (i = 0; i < iter; i++) {
         if (i % relErrStep == 0) {
           LOGGER.info("Iter: " + i);
         } else {
           LOGGER.trace("Iter: " + i);
         }
         // fill seed pixels explicitly with probability 1 for FG and BG
         ArrayRealVector tmp = new ArrayRealVector(1); // filled with 0, setValues() needs that input
         double[] tmpref = tmp.getDataRef(); // just get reference to underlying array
         // set probability to 0 of being FG for BG seeds and vice versa
         for (List<Point> b : seedsPointsBg) {
           setValues(fg, b, tmp); // set 0 all seed pixel currently considered as BG
         }
         // set probability to 1 for of being FG for FG seeds and vice versa
         tmpref[0] = 1;
         setValues(fg, seedsPointsFg.get(cell), tmp); // set 1
         tmp = null;
 
         // ------------------- Computation of FG map ----------------------------------------------
         // groups for long term for FG. Get all four neighbours to current pixel
         Array2DRowRealMatrix fgcircright = (Array2DRowRealMatrix) circshift(fg, RIGHT);
         Array2DRowRealMatrix fgcircleft = (Array2DRowRealMatrix) circshift(fg, LEFT);
         Array2DRowRealMatrix fgcirctop = (Array2DRowRealMatrix) circshift(fg, TOP);
         Array2DRowRealMatrix fgcircbottom = (Array2DRowRealMatrix) circshift(fg, BOTTOM);
 
         // extract required references from RealMatrix objects.
         double[][] fgcircright2d = fgcircright.getDataRef(); // reference to FG prob shifted right
         double[][] fgcircleft2d = fgcircleft.getDataRef(); // reference to FG probab. shifted left
         double[][] fgcirctop2d = fgcirctop.getDataRef(); // reference to FG probab. shifted top
         double[][] fgcircbottom2d = fgcircbottom.getDataRef(); // reference to FG prob. shifted bot
 
         AtomicInteger at = new AtomicInteger(stoppedReason.getValue());
         // Traverse all pixels in FG map and update them according to diffusion from 4 neighbours of
         // each pixel
         int nrows = fg.getRowDimension();
         int ncols = fg.getColumnDimension();
         IntStream.range(0, nrows * ncols).parallel().forEach(ii -> {
           int c = ii % ncols;
           int r = ii / ncols;
           fg2d[r][c] += params.dt * (diffusion * (((fgcircright2d[r][c] - fg2d[r][c]) / wrfg2d[r][c]
                   - (fg2d[r][c] - fgcircleft2d[r][c]) / wlfg2d[r][c])
                   + ((fgcirctop2d[r][c] - fg2d[r][c]) / wtfg2d[r][c]
                           - (fg2d[r][c] - fgcircbottom2d[r][c]) / wbfg2d[r][c])));
           // - params.gamma[currentSweep] * fg2d[r][c] * bg2d[r][c] - disabled
           // validate numerical quality. This flags will stop iterations (break outerloop) but
           // after updating all pixels in FG maps
           if (Double.isNaN(fg2d[r][c])) { // if at least one NaN in solution
             at.set(StoppedBy.NANS.getValue());
           }
           if (Double.isInfinite(fg2d[r][c])) { // or Inf (will overwrite previous flag of course)
             at.set(StoppedBy.INFS.getValue());
           }
         });
 
         // Test state of the flag. Stop iteration if there is NaN or Inf. Iterations are stopped
         // after full looping over FG maps.
         if (stoppedReason == StoppedBy.NANS || stoppedReason == StoppedBy.INFS) {
           break outerloop;
         }
         // check error every relErrStep number of iterations
         if ((i + 1) % relErrStep == 0) {
           double rele = computeRelErr(tmpFglast2d, fg2d);
           LOGGER.info("Relative error for object " + cell + " = " + rele);
           if (rele < params.relim[currentSweep]) {
             stoppedReason = StoppedBy.RELERR;
             // store number of iters for object, required for limiting iterations for BCK
             if (!userBckPoints.contains(seedsPointsFg.get(cell))) {
               iterations.add(i);
             }
             break outerloop;
           }
         }
         // store probabilities over iterations
         if (QuimP.SUPER_DEBUG) {
           debugPm =
                   (debugPm == null) ? new ImageStack(fg.getColumnDimension(), fg.getRowDimension())
                           : debugPm;
           if (i > 1000) {
             if (i % 50 == 0) {
               debugPm.addSlice(QuimPArrayUtils.realMatrix2ImageProcessor(fg));
             }
           } else {
             debugPm.addSlice(QuimPArrayUtils.realMatrix2ImageProcessor(fg));
           }
         }
         // remember FG map for this iteration to use it to compute relative error in next iteration
         QuimPArrayUtils.copy2darray(fg2d, tmpFglast2d);
       } // iter
       LOGGER.info("Sweep " + currentSweep + " for object " + cell + " stopped by " + stoppedReason
               + " after " + i + " iteration from " + iter);
       if (userBckPoints.contains(seedsPointsFg.get(cell))) { // we processed background seeds
         ret.put(SeedTypes.BACKGROUND, fg); // store it in separate key - needed for proper compar.
       } else {
         ret.put(SeedTypes.FOREGROUNDS, fg);
       }
       // save stack of probability maps (over iterations) for each processed object separately
       if (QuimP.SUPER_DEBUG) {
         if (debugPm != null) {
           ImagePlus debugIm = new ImagePlus("debug", debugPm);
           String tmp = System.getProperty("java.io.tmpdir");
           Path p = Paths.get(tmp, "Rw_ProbMap-cell_" + cell);
           IJ.saveAsTiff(debugIm, p.toString());
           debugPm = null; // next object
         }
       }
     } // cell
     return ret;
   }
 
   /**
    * Compute relative error between current foreground and foreground from previous iteration.
    * 
    * <p>Assumes that both matrixes have the same size and they are regular arrays.
    * 
    * @param fglast foreground matrix from previous iteration
    * @param fg current foreground
    * @return relative mean error sum[2* |fg - fglast|/(fg + fglast)]/numofel
    */
   double computeRelErr(double[][] fglast, double[][] fg) {
     int rows = fglast.length;
     int cols = fglast[0].length;
     double rel = 0; // result to sum up
     double tmp = 0; // temporary - current element
     for (int r = 0; r < rows; r++) { // iterate over every element
       for (int c = 0; c < cols; c++) {
         double denominator = fg[r][c] + fglast[r][c]; // compute denominator
         if (denominator == 0.0) { // not divide by 0
           tmp = 0.0;
         } else { // get relative error
           tmp = 2 * Math.abs(fg[r][c] - fglast[r][c]) / denominator;
         }
         rel += tmp; // sum it up to get mean at the end
       }
     }
     double rele = rel / (rows * cols);
     return rele; // return mean error
   }
 
   /**
    * Compute diffusion weights using difference to mean intensity and gradient.
    * 
    * @param diffI2 normalized squared differences to mean seed intensities
    * @param grad2 normalized the squared gradients by the maximum gradient
    * @return wr_fg = exp(alpha*diffI2 + beta*grad2);
    */
   private RealMatrix computeweights(RealMatrix diffI2, RealMatrix grad2) {
     double alpha = params.alpha; // user provided segmentation parameter
     double beta = params.beta; // user provided segmentation parameter
     double[][] diffI2d; // temporary references to intensity to mean values
     double[][] grad22d; // and gradient
     // convert RealMatrix to 2D array, approach depends on the RealMatrix type (see RealMatrix doc)
     if (diffI2 instanceof Array2DRowRealMatrix) {
       diffI2d = ((Array2DRowRealMatrix) diffI2).getDataRef();
     } else {
       diffI2d = diffI2.getData();
     }
     if (grad2 instanceof Array2DRowRealMatrix) {
       grad22d = ((Array2DRowRealMatrix) grad2).getDataRef();
     } else {
       grad22d = grad2.getData();
     }
     Array2DRowRealMatrix w =
             new Array2DRowRealMatrix(diffI2.getRowDimension(), diffI2.getColumnDimension()); // out
     double[][] w2d = w.getDataRef(); // reference of w to skip get/set element from RealMatrix
     // calculate weights
     for (int r = 0; r < diffI2.getRowDimension(); r++) {
       for (int c = 0; c < diffI2.getColumnDimension(); c++) {
         w2d[r][c] = Math.exp(diffI2d[r][c] * alpha + grad22d[r][c] * beta);
       }
     }
 
     return w;
   }
 
   /**
    * Multiply two matrices (element by element) and return minimum value from result.
    * 
    * @param a left operand (matrix)
    * @param b right operand (matrix of size of left operand)
    * @return min(a.*b)
    * @see #getDiffusionConst(RealMatrix, RealMatrix, RealMatrix, RealMatrix, RealMatrix, RealMatrix)
    */
   private double getMinofDotProduct(RealMatrix a, RealMatrix b) {
     RealMatrix cp = a.copy();
     cp.walkInOptimizedOrder(new MatrixDotProduct(b));
     return QuimPArrayUtils.getMin(cp);
   }
 
   /**
    * Compute diffusion constant that obeys stability criterion.
    * 
    * <p>Diffusion constant that meets CFL stability criterion (when solving Euler scheme):
    * D*dt/(dx^2) should be <<1/4
    * It is computed in the following way:
    * 
    * <pre>
    * <code>
    * drl2 = min ( min(min (wr_fg.*avgwx_fg)) , min(min (wl_fg.*avgwx_fg)) );
    * dtb2 = min ( min(min (wt_fg.*avgwy_fg)) , min(min (wb_fg.*avgwy_fg)) );
    * 
    * D=0.25*min(drl2,dtb2);
    * </code>
    * </pre>
    * 
    * @param wrfg weights for diffusion in right direction
    * @param wlfg weights for diffusion in left direction
    * @param wtfg weights for diffusion in top direction
    * @param wbfg weights for diffusion in bottom direction
    * @param avgwxfg averaged left -- right weights
    * @param avgwyfg averaged bottom -- top weights
    * @return diffusion constant that obeys stability criterion
    */
   private double getDiffusionConst(RealMatrix wrfg, RealMatrix wlfg, RealMatrix wtfg,
           RealMatrix wbfg, RealMatrix avgwxfg, RealMatrix avgwyfg) {
     // diffusion constant along x
     double tmp1 = getMinofDotProduct(wrfg, avgwxfg); // min(wrfg .* avgwxfg)
     double tmp2 = getMinofDotProduct(wlfg, avgwxfg); // min(wlfg .* avgwxfg)
     double drl2 = tmp1 < tmp2 ? tmp1 : tmp2; // min(tmp1,tmp2)
     // diffusion constant along y
     tmp1 = getMinofDotProduct(wtfg, avgwyfg); // min(wtfg .* avgwyfg)
     tmp2 = getMinofDotProduct(wbfg, avgwyfg); // min(wbfg .* avgwyfg)
     double dtb2 = tmp1 < tmp2 ? tmp1 : tmp2; // min(tmp1,tmp2)
     double diffusion = drl2 < dtb2 ? drl2 : dtb2; // D = min(drl2,dtb2)
     diffusion *= 0.25; // D = 0.25*min(drl2,dtb2)
     LOGGER.debug("drl2=" + drl2 + " dtb2=" + dtb2);
     return diffusion;
   }
 
   /**
    * Set maximum theoretical squared intensity depending on image type.
    * 
    * @param ip segmented image.
    */
   private void setMaxTheoreticalIntSqr(ImageProcessor ip) {
     switch (ip.getBitDepth()) {
       case 16:
         maxTheoreticalIntSqr = 65535 * 65535;
         break;
       case 8:
       default:
         maxTheoreticalIntSqr = 255 * 255; // default in case of RealMatrix on input
     }
   }
 
   /**
    * Return values from in matrix that are on indexes ind.
    * 
    * @param in Input matrix 2D
    * @param ind List of indexes
    * @return Values that are on indexes ind
    */
   ArrayRealVector getValues(RealMatrix in, List<Point> ind) {
     ArrayRealVector out = new ArrayRealVector(ind.size());
     int l = 0;
     for (Point p : ind) {
       out.setEntry(l++, in.getEntry(p.row, p.col));
     }
     return out;
   }
 
   /**
    * Set values from val on indexes ind in array in.
    * 
    * @param in Input matrix. Will be modified
    * @param ind List of indexes
    * @param val List of values, length must be the same as ind or 1
    */
   void setValues(RealMatrix in, List<Point> ind, ArrayRealVector val) {
     if (ind.size() != val.getDimension() && val.getDimension() != 1) {
       throw new InvalidParameterException(
               "Vector with data must contain 1 element or the same as indexes");
     }
     int delta; // step we move in val across elements, can be 1 or 0
     int l = 0; // will point current element in val
     if (val.getDimension() == 1) {
       delta = 0; // still point to the same element (because it is only one)
     } else {
       delta = 1; // move to the next one
     }
     for (Point p : ind) { // iterate over points to fill with values val
       in.setEntry(p.row, p.col, val.getDataRef()[l]);
       l += delta; // skip to next value (points are iterated in for)
     }
   }
 
   /**
    * Return probability maps for each object (foreground is last). Valid after running the plugin.
    * 
    * @return probability maps
    */
   public ProbabilityMaps getProbabilityMaps() {
     return solved;
   }
 
   /**
    * Sub and then div in-place this matrix and another.
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixDotSubDiv implements RealMatrixChangingVisitor {
 
     /** The sub. */
     RealMatrix sub;
 
     /** The div. */
     RealMatrix div;
 
     /**
      * Instantiates a new matrix dot sub div.
      *
      * @param sub the sub
      * @param div the div
      */
     public MatrixDotSubDiv(RealMatrix sub, RealMatrix div) {
       this.sub = sub;
       this.div = div;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
 
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
 
       return (arg2 - sub.getEntry(arg0, arg1)) / div.getEntry(arg0, arg1);
     }
 
   }
 
   /**
    * Perform element-wise multiplication by value (.*val in Matlab). Done in-place.
    * 
    * @author p.baniukiewicz
    */
   static class MatrixElementMultiply implements RealMatrixChangingVisitor {
 
     /** The multiplier. */
     private double multiplier;
 
     /**
      * Assign multiplier.
      * 
      * @param d multiplier
      */
     MatrixElementMultiply(double d) {
       this.multiplier = d;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
     }
 
     /**
      * Multiply entry by value.
      *
      * @param arg0 the arg 0
      * @param arg1 the arg 1
      * @param arg2 the arg 2
      * @return the double
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return arg2 * multiplier;
     }
 
   }
 
   /**
    * Perform element-wise division by value (./val in Matlab). Done in-place.
    * 
    * @author p.baniukiewicz
    */
   static class MatrixElementDivide implements RealMatrixChangingVisitor {
 
     /** The div. */
     private double div;
 
     /**
      * Assign multiplier.
      * 
      * @param d multiplier
      */
     MatrixElementDivide(double d) {
       this.div = d;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
     }
 
     /**
      * Multiply entry by value.
      *
      * @param arg0 the arg 0
      * @param arg1 the arg 1
      * @param arg2 the arg 2
      * @return the double
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return arg2 / div;
     }
 
   }
 
   /**
    * Perform element-wise exp. Done in-place.
    * 
    * @author p.baniukiewicz
    */
   static class MatrixElementExp implements RealMatrixChangingVisitor {
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
 
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return Math.exp(arg2);
     }
 
   }
 
   /**
    * Perform element-wise power (.^2 in Matlab). Done in-place.
    * 
    * @author p.baniukiewicz
    */
   static class MatrixElementPower implements RealMatrixChangingVisitor {
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
     }
 
     /**
      * Multiply entry by itself.
      *
      * @param arg0 the arg 0
      * @param arg1 the arg 1
      * @param arg2 the arg 2
      * @return the double
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return arg2 * arg2;
     }
   }
 
   /**
    * Multiply in-place this matrix by another.
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixDotProduct implements RealMatrixChangingVisitor {
 
     /** The matrix. */
     RealMatrix matrix;
 
     /**
      * Instantiates a new matrix dot product.
      *
      * @param m the m
      */
     public MatrixDotProduct(RealMatrix m) {
       this.matrix = m;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
       if (matrix.getColumnDimension() != arg1 || matrix.getRowDimension() != arg0) {
         throw new MatrixDimensionMismatchException(matrix.getRowDimension(),
                 matrix.getColumnDimension(), arg0, arg1);
       }
 
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
 
       return arg2 * matrix.getEntry(arg0, arg1);
     }
 
   }
 
   /**
    * Divide in-place this matrix by another.
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixDotDiv implements RealMatrixChangingVisitor {
 
     /** The matrix. */
     RealMatrix matrix;
 
     /**
      * Instantiates a new matrix dot div.
      *
      * @param m the m
      */
     public MatrixDotDiv(RealMatrix m) {
       this.matrix = m;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
       if (matrix.getColumnDimension() != arg1 || matrix.getRowDimension() != arg0) {
         throw new MatrixDimensionMismatchException(matrix.getRowDimension(),
                 matrix.getColumnDimension(), arg0, arg1);
       }
 
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return arg2 / matrix.getEntry(arg0, arg1);
     }
 
   }
 
   /**
    * Divide matrix by matrix (element by element) setting 0 to result if divider is 0. Prevent NaNs.
    * 
    * <p>It performs the operation:
    * 
    * <pre>
    * <code>
    * A = [....];
    * B = [....];
    * R = A./B;
    * R(isnan(R)) = 0;
    * </code>
    * </pre>
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixDotDivN extends MatrixDotDiv {
 
     /**
      * Instantiates a new matrix dot div N.
      *
      * @param m the m
      */
     public MatrixDotDivN(RealMatrix m) {
       super(m);
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see
      * com.github.celldynamics.quimp.plugin.randomwalk.RandomWalkSegmentation.MatrixDotDiv#visit(
      * int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       double entry = matrix.getEntry(arg0, arg1);
       if (entry != 0.0) {
         return arg2 / matrix.getEntry(arg0, arg1);
       } else {
         return 0.0;
       }
     }
   }
 
   /**
    * Perform operation this = this>d*matrix?1:0.
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixCompareWeighted implements RealMatrixChangingVisitor {
 
     /** The matrix. */
     RealMatrix matrix;
 
     /** The weight. */
     double weight;
 
     /**
      * Instantiates a new matrix compare weighted.
      *
      * @param matrix the matrix
      * @param weight the weight
      */
     public MatrixCompareWeighted(RealMatrix matrix, double weight) {
       this.matrix = matrix;
       this.weight = weight;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
       if (matrix.getColumnDimension() != arg1 || matrix.getRowDimension() != arg0) {
         throw new MatrixDimensionMismatchException(matrix.getRowDimension(),
                 matrix.getColumnDimension(), arg0, arg1);
       }
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       if (arg2 > matrix.getEntry(arg0, arg1) * weight) {
         return 1.0;
       } else {
         return 0.0;
       }
     }
 
   }
 
   /**
    * Add in-place this matrix to another.
    * 
    * <p>src/test/Resources-static/Matlab/rw_laplace4_java_base.m This is source file of segmentation
    * in Matlab that was a base for RandomWalkSegmentation Implementation
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixDotAdd implements RealMatrixChangingVisitor {
 
     /** The matrix. */
     RealMatrix matrix;
 
     /**
      * Instantiates a new matrix dot add.
      *
      * @param m the m
      */
     public MatrixDotAdd(RealMatrix m) {
       this.matrix = m;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
       if (matrix.getColumnDimension() != arg1 || matrix.getRowDimension() != arg0) {
         throw new MatrixDimensionMismatchException(matrix.getRowDimension(),
                 matrix.getColumnDimension(), arg0, arg1);
       }
 
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return arg2 + matrix.getEntry(arg0, arg1);
     }
 
   }
 
   /**
    * Sub in-place this matrix to another.
    * 
    * @author p.baniukiewicz
    *
    */
   static class MatrixDotSub implements RealMatrixChangingVisitor {
 
     /**
      * Example src/test/Resources-static/Matlab/rw_laplace4_java_base.m This is source file of
      * segmentation in Matlab that was a base for RandomWalkSegmentation Implementation.
      */
     RealMatrix matrix;
 
     /**
      * Instantiates a new matrix dot sub.
      *
      * @param m the m
      */
     public MatrixDotSub(RealMatrix m) {
       this.matrix = m;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
       if (matrix.getColumnDimension() != arg1 || matrix.getRowDimension() != arg0) {
         throw new MatrixDimensionMismatchException(matrix.getRowDimension(),
                 matrix.getColumnDimension(), arg0, arg1);
       }
 
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#visit(int, int, double)
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return arg2 - matrix.getEntry(arg0, arg1);
     }
 
   }
 
   /**
    * Perform element-wise power (.^2 in Matlab) and then divide by val. Done in-place.
    * 
    * @author p.baniukiewicz
    */
   static class MatrixElementPowerDiv implements RealMatrixChangingVisitor {
 
     /** The val. */
     double val;
 
     /**
      * Instantiates a new matrix element power div.
      *
      * @param val the val
      */
     public MatrixElementPowerDiv(double val) {
       this.val = val;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#end()
      */
     @Override
     public double end() {
       return 0;
     }
 
     /*
      * (non-Javadoc)
      * 
      * @see org.apache.commons.math3.linear.RealMatrixChangingVisitor#start(int, int, int, int, int,
      * int)
      */
     @Override
     public void start(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5) {
     }
 
     /**
      * Multiply entry by itself.
      *
      * @param arg0 the arg 0
      * @param arg1 the arg 1
      * @param arg2 the arg 2
      * @return the double
      */
     @Override
     public double visit(int arg0, int arg1, double arg2) {
       return (arg2 * arg2) / val;
     }
   }
 
 }