package stainingestimation;

import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.WindowManager;
import ij.gui.GenericDialog;
import ij.gui.ImageCanvas;
import ij.gui.ImageWindow;
import ij.gui.NewImage;
import ij.plugin.PlugIn;
import ij.process.ImageConverter;
import ij.process.ImageProcessor;
import java.awt.Color;
import java.awt.Font;
import java.awt.Point;
import java.awt.image.IndexColorModel;
import java.util.ArrayList;
import javax.swing.BoxLayout;
import org.jdesktop.swingx.JXColorSelectionButton;
import tmarker.TMAspot.TMALabel;
import tmarker.TMAspot.TMAspot;


public class Colour_Deconvolution implements PlugIn {

// G.Landini at bham ac uk
// 30/Mar/2004 released
// 03/Apr/2004 resolved ROI exiting
// 07/Apr/2004 added Methyl Green DAB vectors
// 08/Jul/2004 shortened the code
// 01/Aug/2005 added fast red/blue/DAB vectors
// 02/Nov/2005 changed code to work with image stacks (DLC - dchao at fhcrc org)
// 02/Nov/2005 changed field names so user-defined colours can be set within 
//             macros (DLC - dchao at fhcrc org)
// 04/Feb/2007 1.3 disable popup menu when right clicking
// 23/May/2009 added Feulgen-light green vectors
//14/Apr/2010 v 1.4 added Giemsa vector (Methylene blue & eosin) 
//           the images are now names "title"-(Colour_1) etc so there are not clash of names when using [ ]
//           the log window now prints the java code of the translation matrix to include new vectors in the plugin.
//           added "Hide legend" option
//22/Jun/2010 v 1.5 added Masson Trichrome vector (Methyl blue & Ponceau Fuchsin only (this does not have Iron Haematoxylin vector!)
//           fixed bug: check for 0 components before hiding legend (otherwise there was no image shown if legent hidden) 
//
// This plugin implements stain separation using the colour deconvolution
// method described in:
//
//     Ruifrok AC, Johnston DA. Quantification of histochemical
//     staining by color deconvolution. Analytical & Quantitative
//     Cytology & Histology 2001; 23: 291-299.
//
// The code is based on "Color separation-30", a macro for NIH Image kindly provided
// by A.C. Ruifrok. Thanks Arnout!
//
// The plugin assumes images generated by color subtraction (i.e. light-absorbing dyes
// such as those used in bright field histology or ink on printed paper) but the dyes
// should not be neutral grey.
//
// I strongly suggest to read the paper reference above to understand how to determine
// new vectors and how the whole procedure works.
//
// The plugin works correctly when the background is neutral (white to light grey), 
// so background subtraction and colour correction must be applied to the images before 
// processing.
//
// The plugin provides a number of "built in" stain vectors some of which were determined
// experimentally in our lab (marked GL), but you may have to determine your own vectors to
// provide a more accurate stain separation, depending on the stains and methods you use.
// Ideally, vector determination should be done on slides stained with only one colour
// at a time (using the "From ROI" interactive option).
//
// The plugin takes an RGB image and returns three 8-bit images. If the specimen is
// stained with a 2 colour scheme (such as H & E) the 3rd image represents the
// complimentary of the first two colours (i.e. green).
//
// Please be *very* careful about how to interpret the results of colour deconvolution
// when analysing histological images.
// Most staining methods are not stochiometric and so optical density of the chromogen
// may not correlate well with the *quantity* of the reactants.
// This means that optical density of the colour may not be a good indicator of
// the amount of material stained.
//
// Read the paper!
//

        /**
         * Calculates and returns the channel images for a given color image.
         * @param imp The image to be deconvolved.
         * @param se The StainingEstimation instance, used for user value input, if myStaing = "User values".
         * @param ts The currently visible TMAspots, used for user value input, if myStaing = "User values".
         * @param hideLegend If true, the legend of the color deconvolution method is hidden.
         * @param myStain The staining protocol string (e.g. "H&E" or "H DAB").
         * @param substractChannels If true, channel 2 will be substracted from channel 1.
         * @param invert_ch1 If true, channel 1 gets inverted.
         * @param invert_ch2 If true, channel 2 gets inverted.
         * @param invert_ch3 If true, channel 3 gets inverted.
         * @return The three calculated channel images of the deconvolution.
         */
	public static java.util.List<ImagePlus> get_deconvolution_images(ImagePlus imp, StainingEstimation se, TMAspot ts, boolean hideLegend, String myStain, boolean substractChannels, boolean invert_ch1, boolean invert_ch2, boolean invert_ch3) {
		if (imp==null){
			IJ.error("No image!");
			return null;
		}
                
                if (imp.getBitDepth()!=24){
			ImageConverter ic = new ImageConverter(imp);
                        ic.convertToRGB();
                        //IJ.error("RGB image needed.");
			//return null;
		}
                
                ImagePlus imp_tmp = imp.duplicate();
                // normalize the orginal Image
                //StainingEstimation.normalize_image_rgb(imp_tmp);
                
                java.util.List<ImagePlus> d_imps = new ArrayList<>(3);
                String title = imp_tmp.getTitle();
                // "Black / White" does not mean deconvolution with black and white color vectors (which will result in monotone gray images).
                // But it means "just treat the images as black white images".
                if (myStain.equals("Black / White")){
                    d_imps.add(new ImagePlus(title+"-(Colour_1)",imp_tmp.getImage()));
                    // inverted Image
                    ImagePlus ip_inv = imp_tmp.duplicate();
                    ip_inv.getProcessor().invert();
                    d_imps.add(ip_inv);
                    // original image as third channel
                    d_imps.add(new ImagePlus(title+"-(Colour_3)",imp_tmp.getImage()));
                    return d_imps;
                }
                
		ImageStack stack = imp_tmp.getStack();
		int width = stack.getWidth();
		int height = stack.getHeight();
                
		//GenericDialog gd = new GenericDialog("Colour Deconvolution 1.5", parent);
		//gd.addMessage("Select subtractive colour model");
		//String [] stain={"From ROI", "H&E", "H&E 2","H DAB", "Feulgen Light Green", "Giemsa", "FastRed FastBlue DAB", "Methyl Green DAB", "H&E DAB", "H AEC","Azan-Mallory","Masson Trichrome","Alcian blue & H","H PAS","RGB","CMY", "User values"};
		//String [] stain={              "H&E", "H&E 2","H DAB", "Feulgen Light Green", "Giemsa", "FastRed FastBlue DAB", "Methyl Green DAB", "H&E DAB", "H AEC","Azan-Mallory","Masson Trichrome","Alcian blue & H","H PAS","RGB","CMY", "User values"};
		//gd.addChoice("Vectors", stain, stain[0]);
		//gd.addCheckbox("Show matrices",false);
		//gd.addCheckbox("Hide legend",true);

		//gd.addCheckbox("Linear image",false);

		//gd.showDialog();
		//if (gd.wasCanceled())
		//	return null;
		boolean doIshow = false;//gd.getNextBoolean();
		
		double leng, A, V, C, log255=Math.log(255.0);
		int i,j;
		double [] MODx = new double[3];
		double [] MODy = new double[3];
		double [] MODz = new double[3];
		double [] cosx = new double[3];
		double [] cosy = new double[3];
		double [] cosz  = new double[3];
		double [] len = new double[3];
		double [] q = new double[9];
		byte [] rLUT = new byte[256];
		byte [] gLUT = new byte[256];
		byte [] bLUT = new byte[256];
                
                getColorVector(myStain, MODx, MODy, MODz, se, ts);
                
                boolean contine = false;
                for (i=0; i<3; i++) {
                    contine = contine || MODx[i]!=0;
                    contine = contine || MODy[i]!=0;
                    contine = contine || MODz[i]!=0;                    
                }
                if (!contine) {
                    return d_imps;
                }
                
                // start
		for (i=0; i<3; i++){
			//normalise vector length
			cosx[i]=cosy[i]=cosz[i]=0.0;
			len[i]=Math.sqrt(MODx[i]*MODx[i] + MODy[i]*MODy[i] + MODz[i]*MODz[i]);
			if (len[i] != 0.0){
				cosx[i]= MODx[i]/len[i];
				cosy[i]= MODy[i]/len[i];
				cosz[i]= MODz[i]/len[i];
			}
		}


		// translation matrix
		if (cosx[1]==0.0){ //2nd colour is unspecified
			if (cosy[1]==0.0){
				if (cosz[1]==0.0){
					cosx[1]=cosz[0];
					cosy[1]=cosx[0];
					cosz[1]=cosy[0];
				}
			}
		}

		if (cosx[2]==0.0){ // 3rd colour is unspecified
			if (cosy[2]==0.0){
				if (cosz[2]==0.0){
					if ((cosx[0]*cosx[0] + cosx[1]*cosx[1])> 1){
						if (doIshow)
							IJ.log("Colour_3 has a negative R component.");
						cosx[2]=0.0;
					}
					else {
						cosx[2]=Math.sqrt(1.0-(cosx[0]*cosx[0])-(cosx[1]*cosx[1]));
					}

					if ((cosy[0]*cosy[0] + cosy[1]*cosy[1])> 1){
						if (doIshow)
							IJ.log("Colour_3 has a negative G component.");
						cosy[2]=0.0;
					}
					else {
						cosy[2]=Math.sqrt(1.0-(cosy[0]*cosy[0])-(cosy[1]*cosy[1]));
					}

					if ((cosz[0]*cosz[0] + cosz[1]*cosz[1])> 1){
						if (doIshow)
							IJ.log("Colour_3 has a negative B component.");
						cosz[2]=0.0;
					}
					else {
						cosz[2]=Math.sqrt(1.0-(cosz[0]*cosz[0])-(cosz[1]*cosz[1]));
					}
				}
			}
		}

		leng=Math.sqrt(cosx[2]*cosx[2] + cosy[2]*cosy[2] + cosz[2]*cosz[2]);

		cosx[2]= cosx[2]/leng;
		cosy[2]= cosy[2]/leng;
		cosz[2]= cosz[2]/leng;

		for (i=0; i<3; i++){
			if (cosx[i] == 0.0) cosx[i] = 0.001;
			if (cosy[i] == 0.0) cosy[i] = 0.001;
			if (cosz[i] == 0.0) cosz[i] = 0.001;
		}

		if (!hideLegend) {
			ImagePlus imp0 = NewImage.createRGBImage("Colour Deconvolution", 350, 65, 1, 0);
			ImageProcessor ip0 = imp0.getProcessor();
			ip0.setFont(new Font("Monospaced", Font.BOLD, 11));
			ip0.setAntialiasedText(true);
			ip0.setColor(Color.black);
			ip0.moveTo(10,15);
			ip0.drawString("Colour deconvolution: "+myStain);
			ip0.setFont(new Font("Monospaced", Font.PLAIN, 10));

			for (i=0; i<3; i++){
				ip0.setRoi(10,18+ i*15, 14, 14); 
				ip0.setColor( 
					(((255 -(int)(255.0* cosx[i])) & 0xff)<<16)+
					(((255 -(int)(255.0* cosy[i])) & 0xff)<<8 )+
					(((255 -(int)(255.0* cosz[i])) & 0xff)    ));
				ip0.fill();
				ip0.setFont(new Font("Monospaced", Font.PLAIN, 10));
				ip0.setAntialiasedText(true);
				ip0.setColor(Color.black);
				ip0.moveTo(27,32+ i*15);
				ip0.drawString("Colour_"+(i+1)+" R:"+(float)cosx[i]+", G:"+(float)cosy[i]+", B:"+(float)cosz[i] );
			}
			imp0.show();
			imp0.updateAndDraw();
		}

		if (myStain.equals("From ROI")){
			IJ.showMessage("Vectors","Colour[1]:\n"+
			"  R: "+(float) cosx[0] +"\n"+
			"  G: "+(float) cosy[0] +"\n"+
			"  B: "+(float) cosz[0] +"\n \n"+
			"Colour[2]:\n"+
			"  R: "+ (float) cosx[1] +"\n"+
			"  G: "+ (float) cosy[1] +"\n"+
			"  B: "+ (float) cosz[1] +"\n \n"+
			"Colour[3]:\n"+
			"  R: "+ (float) cosx[2] +"\n"+
			"  G: "+ (float) cosy[2] +"\n"+
			"  B: "+ (float) cosz[2]);
		}


		if (doIshow){
			IJ.log( myStain +" Vector Matrix ---");
			for (i=0; i<3; i++){
				IJ.log("Colour["+(i+1)+"]:\n"+
				"  R"+(i+1)+": "+ (float) MODx[i] +"\n"+
				"  G"+(i+1)+": "+ (float) MODy[i] +"\n"+
				"  B"+(i+1)+": "+ (float) MODz[i] +"\n \n");
			}

/*			IJ.log( myStain +" Translation Matrix ---");
			for (i=0; i<3; i++){
				IJ.log("Colour["+(i+1)+"]:\n"+
				"  R"+(i+1)+": "+ (float) cosx[i] +"\n"+
				"  G"+(i+1)+": "+ (float) cosy[i] +"\n"+
				"  B"+(i+1)+": "+ (float) cosz[i] +"\n \n");
			}
*/
			IJ.log( myStain +" Java code ---");
			IJ.log("\t\tif (myStain.equals(\"New_Stain\")){");
			IJ.log("\t\t// This is the New_Stain");
			for (i=0; i<3; i++){
				IJ.log("\t\t\tMODx["+i+"]="+ (float) cosx[i] +";\n"+
					"\t\t\tMODy["+i+"]="+ (float) cosy[i] +";\n"+
					"\t\t\tMODz["+i+"]="+ (float) cosz[i] +";\n\n");
			}
			IJ.log("}");
		}

		//matrix inversion
		A = cosy[1] - cosx[1] * cosy[0] / cosx[0];
		V = cosz[1] - cosx[1] * cosz[0] / cosx[0];
		C = cosz[2] - cosy[2] * V/A + cosx[2] * (V/A * cosy[0] / cosx[0] - cosz[0] / cosx[0]);
		q[2] = (-cosx[2] / cosx[0] - cosx[2] / A * cosx[1] / cosx[0] * cosy[0] / cosx[0] + cosy[2] / A * cosx[1] / cosx[0]) / C;
		q[1] = -q[2] * V / A - cosx[1] / (cosx[0] * A);
		q[0] = 1.0 / cosx[0] - q[1] * cosy[0] / cosx[0] - q[2] * cosz[0] / cosx[0];
		q[5] = (-cosy[2] / A + cosx[2] / A * cosy[0] / cosx[0]) / C;
		q[4] = -q[5] * V / A + 1.0 / A;
		q[3] = -q[4] * cosy[0] / cosx[0] - q[5] * cosz[0] / cosx[0];
		q[8] = 1.0 / C;
		q[7] = -q[8] * V / A;
		q[6] = -q[7] * cosy[0] / cosx[0] - q[8] * cosz[0] / cosx[0];

                // initialize 3 output colour stacks
                ImageStack[] outputstack = new ImageStack[3];
		for (i=0; i<3; i++){
			for (j=0; j<256; j++) { //LUT[1]
				//if (cosx[i] < 0)
				//	rLUT[255-j]=(byte)(255.0 + (double)j * cosx[i]);
				//else
					rLUT[255-j]=(byte)(255.0 - (double)j * cosx[i]);

				//if (cosy[i] < 0)
				//	gLUT[255-j]=(byte)(255.0 + (double)j * cosy[i]);
				//else
					gLUT[255-j]=(byte)(255.0 - (double)j * cosy[i]);

				//if (cosz[i] < 0)
				//	bLUT[255-j]=(byte)(255.0 + (double)j * cosz[i]);
				///else
					bLUT[255-j]=(byte)(255.0 - (double)j * cosz[i]);
			}
                        IndexColorModel cm = new IndexColorModel(8, 256, rLUT, gLUT, bLUT);
                        outputstack[i] = new ImageStack(width, height, cm);
                }

		// translate ------------------
		int imagesize = width * height;
		for (int imagenum=1; imagenum<=stack.getSize(); imagenum++) {
			int[] pixels = (int[])stack.getPixels(imagenum);
			String label = stack.getSliceLabel(imagenum);
			byte[][] newpixels = new byte[3][];
			newpixels[0] = new byte[imagesize];
			newpixels[1] = new byte[imagesize];
			newpixels[2] = new byte[imagesize];

			for (j=0;j<imagesize;j++){
				 // log transform the RGB data
				int R = (pixels[j] & 0xff0000)>>16;
				int G = (pixels[j] & 0x00ff00)>>8 ;
				int B = (pixels[j] & 0x0000ff);
				double Rlog = -((255.0*Math.log(((double)R+1)/255.0))/log255);
				double Glog = -((255.0*Math.log(((double)G+1)/255.0))/log255);
				double Blog = -((255.0*Math.log(((double)B+1)/255.0))/log255);
				for (i=0; i<3; i++){
                                    // rescale to match original paper values
                                    double Rscaled = Rlog * q[i*3];
                                    double Gscaled = Glog * q[i*3+1];
                                    double Bscaled = Blog * q[i*3+2];
                                    double output = Math.exp(-((Rscaled + Gscaled + Bscaled) - 255.0) * log255 / 255.0);
                                    if(output>255) {
                                        output=255;
                                    }
                                    newpixels[i][j]=(byte)(0xff&(int)(Math.floor(output+.5)));
				}
			}
                        
                        // substract channel 2 from channel 1 for H DAB staining
                        if (substractChannels) {
                            for (j=0;j<imagesize;j++) {
                                byte c1 = newpixels[0][j];
                                byte c2 = newpixels[1][j];
                                int ch1i = 0x000000ff&c1;
                                int ch2i = 0x000000ff&c2;
                                int ch1ii = 255-ch1i;
                                int ch2ii = 255-ch2i;
                                int d = ch1ii-ch2ii;
                                int df = Math.min(255, Math.max(0, d));
                                int dfi = 255-Math.min(255, Math.max(0, df));
                                byte db = (byte) (0xff&(int)(dfi));
                                newpixels[0][j] = db;
                            }
                        }
                
			 // add new values to output images
			outputstack[0].addSlice(label,newpixels[0]);
			outputstack[1].addSlice(label,newpixels[1]);
			outputstack[2].addSlice(label,newpixels[2]);
		}
		//new ImagePlus(title+"-(Colour_1)",outputstack[0]).show();
		//new ImagePlus(title+"-(Colour_2)",outputstack[1]).show();
		//new ImagePlus(title+"-(Colour_3)",outputstack[2]).show();
                
                d_imps.add(new ImagePlus(title+"-(Colour_1)",outputstack[0]));
                d_imps.add(new ImagePlus(title+"-(Colour_2)",outputstack[1]));
                d_imps.add(new ImagePlus(title+"-(Colour_3)",outputstack[2]));
                
                if (invert_ch1) d_imps.get(0).getProcessor().invert();
                if (invert_ch2) d_imps.get(1).getProcessor().invert();
                if (invert_ch3) d_imps.get(2).getProcessor().invert();
                
                return d_imps;
	}
        
        /**
         * Generates the color vectors based on the staining protocol myStain.
         * @param myStain The staining protocol description, e.g. "H&E", "H&E 2", "H DAB", "Feulgen Light Green" or "Giemsa" "FastRed FastBlue DAB" (or others).
         * @param MODx The RGB color array which is beeing filled for channel 1, size 3.
         * @param MODy The RGB color array which is beeing filled for channel 2, size 3.
         * @param MODz The RGB color array which is beeing filled for channel 3, size 3.
         * @param se The StainingEstimation instance, used for user value input, if myStaing = "User values".
         * @param ts The currently visible TMAspots, used for user value input, if myStaing = "User values".
         */
        public static void getColorVector(String myStain, double[] MODx, double[] MODy, double[] MODz, StainingEstimation se, TMAspot ts) {
            // stains are defined after this line

		if (myStain.equals("H&E")){
			// GL Haem matrix
			MODx[0]= 0.644211; //0.650;
			MODy[0]= 0.716556; //0.704;
			MODz[0]= 0.266844; //0.286;
			// GL Eos matrix
			MODx[1]= 0.092789; //0.072;
			MODy[1]= 0.954111; //0.990;
			MODz[1]= 0.283111; //0.105;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H&E 2")){
			// GL Haem matrix
			MODx[0]= 0.49015734;
			MODy[0]= 0.76897085;
			MODz[0]= 0.41040173;
			// GL Eos matrix
			MODx[1]= 0.04615336;
			MODy[1]= 0.8420684;
			MODz[1]= 0.5373925;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H DAB")){
			// 3,3-diamino-benzidine tetrahydrochloride
			// Haem matrix
			MODx[0]= 0.650;
			MODy[0]= 0.704;
			MODz[0]= 0.286;
			// DAB matrix
			MODx[1]= 0.268;
			MODy[1]= 0.570;
			MODz[1]= 0.776;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("Feulgen Light Green")){
			//GL Feulgen & light green
			//Feulgen
			MODx[0]= 0.46420921;
			MODy[0]= 0.83008335;
			MODz[0]= 0.30827187;
			// light green
			MODx[1]= 0.94705542;
			MODy[1]= 0.25373821;
			MODz[1]= 0.19650764;
			// Zero matrix
			MODx[2]= 0.0; // 0.0010000
			MODy[2]= 0.0; // 0.47027777
			MODz[2]= 0.0; //0.88235928
		}

		if (myStain.equals("Giemsa")){
			// GL  Methylene Blue and Eosin
			MODx[0]= 0.834750233;
			MODy[0]= 0.513556283;
			MODz[0]= 0.196330403;
			// GL Eos matrix
			MODx[1]= 0.092789;
			MODy[1]= 0.954111;
			MODz[1]= 0.283111;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("FastRed FastBlue DAB")){
			//fast red
			MODx[0]= 0.21393921;
			MODy[0]= 0.85112669;
			MODz[0]= 0.47794022;
			// fast blue
			MODx[1]= 0.74890292;
			MODy[1]= 0.60624161;
			MODz[1]= 0.26731082;
			// dab
			MODx[2]= 0.268;
			MODy[2]= 0.570;
			MODz[2]= 0.776;
		}

		if (myStain.equals("Methyl Green DAB")){
			// MG matrix (GL)
			MODx[0]= 0.98003;
			MODy[0]= 0.144316;
			MODz[0]= 0.133146;
			// DAB matrix
			MODx[1]= 0.268;
			MODy[1]= 0.570;
			MODz[1]= 0.776;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H&E DAB")){
			// Haem matrix
			MODx[0]= 0.650;
			MODy[0]= 0.704;
			MODz[0]= 0.286;
			// Eos matrix
			MODx[1]= 0.072;
			MODy[1]= 0.990;
			MODz[1]= 0.105;
			// DAB matrix
			MODx[2]= 0.268;
			MODy[2]= 0.570;
			MODz[2]= 0.776;
		}

		if (myStain.equals("H AEC")){
			// 3-amino-9-ethylcarbazole
			// Haem matrix
			MODx[0]= 0.650;
			MODy[0]= 0.704;
			MODz[0]= 0.286;
			// AEC matrix
			MODx[1]= 0.2743;
			MODy[1]= 0.6796;
			MODz[1]= 0.6803;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("Azan-Mallory")){
			//Azocarmine and Aniline Blue (AZAN)
			// GL Blue matrix Anilline Blue
			MODx[0]= .853033;
			MODy[0]= .508733;
			MODz[0]= .112656;
			// GL Red matrix Azocarmine
			MODx[1]=0.09289875;
			MODy[1]=0.8662008;
			MODz[1]=0.49098468;
			//GL  Orange matrix Orange-G
			MODx[2]=0.10732849;
			MODy[2]=0.36765403;
			MODz[2]=0.9237484;
		}

		if (myStain.equals("Masson Trichrome")){
			// GL Methyl blue
			MODx[0]=0.7995107;
			MODy[0]=0.5913521;
			MODz[0]=0.10528667;
			// GL Ponceau Fuchsin has 2 hues, really this is only approximate
			MODx[1]=0.09997159;
			MODy[1]=0.73738605;
			MODz[1]=0.6680326;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
			// GL Iron Haematoxylin, but this does not seem to work well because it gets confused with the other 2 components
//			MODx[2]=0.6588232;
//			MODy[2]=0.66414213;
//			MODz[2]=0.3533655;
		}

		if (myStain.equals("Alcian blue & H")){
			// GL Alcian Blue matrix
			MODx[0]= 0.874622;
			MODy[0]= 0.457711;
			MODz[0]= 0.158256;
			// GL Haematox after PAS matrix
			MODx[1]= 0.552556;
			MODy[1]= 0.7544;
			MODz[1]= 0.353744;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H PAS")){
			// GL Haem matrix
			MODx[0]= 0.644211; //0.650;
			MODy[0]= 0.716556; //0.704;
			MODz[0]= 0.266844; //0.286;
			// GL PAS matrix
			MODx[1]= 0.175411;
			MODy[1]= 0.972178;
			MODz[1]= 0.154589;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("RGB")){
			//R
			MODx[0]= 0.0;
			MODy[0]= 1.0;
			MODz[0]= 1.0;
			//G
			MODx[1]= 1.0;
			MODy[1]= 0.0;
			MODz[1]= 1.0;
			//B
			MODx[2]= 1.0;
			MODy[2]= 1.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("CMY")){
			//C
			MODx[0]= 1.0;
			MODy[0]= 0.0;
			MODz[0]= 0.0;
			//M
			MODx[1]= 0.0;
			MODy[1]= 1.0;
			MODz[1]= 0.0;
			//Y
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 1.0;
		}

		if (myStain.equals("User values")){
                    if (se != null) {
                        Color c = se.getParam_UserColor1();
                        MODx[0] = 1.0 - c.getRed()/255.0;
                        MODy[0] = 1.0 - c.getGreen()/255.0;
                        MODz[0] = 1.0 - c.getBlue()/255.0;

                        c = se.getParam_UserColor2();
                        MODx[1] = 1.0 - c.getRed()/255.0;
                        MODy[1] = 1.0 - c.getGreen()/255.0;
                        MODz[1] = 1.0 - c.getBlue()/255.0;

                        c = se.getParam_UserColor3();
                        MODx[2] = 1.0 - c.getRed()/255.0;
                        MODy[2] = 1.0 - c.getGreen()/255.0;
                        MODz[2] = 1.0 - c.getBlue()/255.0;
                    } else {
			ArrayList<double[]> uc = getUserColors(se, ts);
                        if (uc!=null) {
                            for (int i=0; i<3; i++) {
                                MODx[i]=1-uc.get(0)[i]/255.0;
                                MODy[i]=1-uc.get(1)[i]/255.0;
                                MODz[i]=1-uc.get(2)[i]/255.0;
                            }
                        }
                    }
		}
                
                // Black / White is handled in function get_deconvolution_images()
                if (myStain.equals("Black / White")){
			//B
			MODx[0]= 0.0;
			MODy[0]= 0.0;
			MODz[0]= 0.0;
			//W
			MODx[1]= 1.0;
			MODy[1]= 1.0;
			MODz[1]= 1.0;
			//Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		/*if (myStain.equals("From ROI")){
			IJ.runMacro("setOption('DisablePopupMenu', true)");
			// imp.getCanvas().disablePopupMenu(true);
			double [] rgbOD = new double[3];
			for (int i=0; i<3; i++){
				getmeanRGBODfromROI(i, rgbOD, imp);
				MODx[i]= rgbOD[0];
				MODy[i]= rgbOD[1];
				MODz[i]= rgbOD[2];
			}

			IJ.runMacro("setOption('DisablePopupMenu', false)");
			// imp.getCanvas().disablePopupMenu(false);

		}*/
        }
        
        @Override
        public void run(String arg) {
		ImagePlus imp = WindowManager.getCurrentImage();
		if (imp==null){
			IJ.error("No image!");
			return;
		}
		if (imp.getBitDepth()!=24){
			IJ.error("RGB image needed.");
			return;
		}
		ImageStack stack = imp.getStack();
		int width = stack.getWidth();
		int height = stack.getHeight();
                String title = imp.getTitle();

		GenericDialog gd = new GenericDialog("Colour Deconvolution 1.5", IJ.getInstance());
		//gd.addMessage("Select subtractive colour model");
		String [] stain={"From ROI", "H&E", "H&E 2","H DAB", "Feulgen Light Green", "Giemsa", "FastRed FastBlue DAB", "Methyl Green DAB", "H&E DAB", "H AEC","Azan-Mallory","Masson Trichrome","Alcian blue & H","H PAS","RGB","CMY", "User values"};
		gd.addChoice("Vectors", stain, stain[0]);
		gd.addCheckbox("Show matrices",false);
		gd.addCheckbox("Hide legend",false);

		//gd.addCheckbox("Linear image",false);

		gd.showDialog();
		if (gd.wasCanceled())
			return;
		String myStain = gd.getNextChoice();
		boolean doIshow = gd.getNextBoolean();
		boolean hideLegend = gd.getNextBoolean();

		double leng, A, V, C, log255=Math.log(255.0);
		int i,j;
		double [] MODx = new double[3];
		double [] MODy = new double[3];
		double [] MODz = new double[3];
		double [] cosx = new double[3];
		double [] cosy = new double[3];
		double [] cosz  = new double[3];
		double [] len = new double[3];
		double [] q = new double[9];
		byte [] rLUT = new byte[256];
		byte [] gLUT = new byte[256];
		byte [] bLUT = new byte[256];

// stains are defined after this line

		if (myStain.equals("H&E")){
			// GL Haem matrix
			MODx[0]= 0.644211; //0.650;
			MODy[0]= 0.716556; //0.704;
			MODz[0]= 0.266844; //0.286;
			// GL Eos matrix
			MODx[1]= 0.092789; //0.072;
			MODy[1]= 0.954111; //0.990;
			MODz[1]= 0.283111; //0.105;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H&E 2")){
			// GL Haem matrix
			MODx[0]= 0.49015734;
			MODy[0]= 0.76897085;
			MODz[0]= 0.41040173;
			// GL Eos matrix
			MODx[1]= 0.04615336;
			MODy[1]= 0.8420684;
			MODz[1]= 0.5373925;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H DAB")){
			// 3,3-diamino-benzidine tetrahydrochloride
			// Haem matrix
			MODx[0]= 0.650;
			MODy[0]= 0.704;
			MODz[0]= 0.286;
			// DAB matrix
			MODx[1]= 0.268;
			MODy[1]= 0.570;
			MODz[1]= 0.776;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("Feulgen Light Green")){
			//GL Feulgen & light green
			//Feulgen
			MODx[0]= 0.46420921;
			MODy[0]= 0.83008335;
			MODz[0]= 0.30827187;
			// light green
			MODx[1]= 0.94705542;
			MODy[1]= 0.25373821;
			MODz[1]= 0.19650764;
			// Zero matrix
			MODx[2]= 0.0; // 0.0010000
			MODy[2]= 0.0; // 0.47027777
			MODz[2]= 0.0; //0.88235928
		}

		if (myStain.equals("Giemsa")){
			// GL  Methylene Blue and Eosin
			MODx[0]= 0.834750233;
			MODy[0]= 0.513556283;
			MODz[0]= 0.196330403;
			// GL Eos matrix
			MODx[1]= 0.092789;
			MODy[1]= 0.954111;
			MODz[1]= 0.283111;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("FastRed FastBlue DAB")){
			//fast red
			MODx[0]= 0.21393921;
			MODy[0]= 0.85112669;
			MODz[0]= 0.47794022;
			// fast blue
			MODx[1]= 0.74890292;
			MODy[1]= 0.60624161;
			MODz[1]= 0.26731082;
			// dab
			MODx[2]= 0.268;
			MODy[2]= 0.570;
			MODz[2]= 0.776;
		}

		if (myStain.equals("Methyl Green DAB")){
			// MG matrix (GL)
			MODx[0]= 0.98003;
			MODy[0]= 0.144316;
			MODz[0]= 0.133146;
			// DAB matrix
			MODx[1]= 0.268;
			MODy[1]= 0.570;
			MODz[1]= 0.776;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H&E DAB")){
			// Haem matrix
			MODx[0]= 0.650;
			MODy[0]= 0.704;
			MODz[0]= 0.286;
			// Eos matrix
			MODx[1]= 0.072;
			MODy[1]= 0.990;
			MODz[1]= 0.105;
			// DAB matrix
			MODx[2]= 0.268;
			MODy[2]= 0.570;
			MODz[2]= 0.776;
		}

		if (myStain.equals("H AEC")){
			// 3-amino-9-ethylcarbazole
			// Haem matrix
			MODx[0]= 0.650;
			MODy[0]= 0.704;
			MODz[0]= 0.286;
			// AEC matrix
			MODx[1]= 0.2743;
			MODy[1]= 0.6796;
			MODz[1]= 0.6803;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("Azan-Mallory")){
			//Azocarmine and Aniline Blue (AZAN)
			// GL Blue matrix Anilline Blue
			MODx[0]= .853033;
			MODy[0]= .508733;
			MODz[0]= .112656;
			// GL Red matrix Azocarmine
			MODx[1]=0.09289875;
			MODy[1]=0.8662008;
			MODz[1]=0.49098468;
			//GL  Orange matrix Orange-G
			MODx[2]=0.10732849;
			MODy[2]=0.36765403;
			MODz[2]=0.9237484;
		}

		if (myStain.equals("Masson Trichrome")){
			// GL Methyl blue
			MODx[0]=0.7995107;
			MODy[0]=0.5913521;
			MODz[0]=0.10528667;
			// GL Ponceau Fuchsin has 2 hues, really this is only approximate
			MODx[1]=0.09997159;
			MODy[1]=0.73738605;
			MODz[1]=0.6680326;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
			// GL Iron Haematoxylin, but this does not seem to work well because it gets confused with the other 2 components
//			MODx[2]=0.6588232;
//			MODy[2]=0.66414213;
//			MODz[2]=0.3533655;
		}

		if (myStain.equals("Alcian blue & H")){
			// GL Alcian Blue matrix
			MODx[0]= 0.874622;
			MODy[0]= 0.457711;
			MODz[0]= 0.158256;
			// GL Haematox after PAS matrix
			MODx[1]= 0.552556;
			MODy[1]= 0.7544;
			MODz[1]= 0.353744;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("H PAS")){
			// GL Haem matrix
			MODx[0]= 0.644211; //0.650;
			MODy[0]= 0.716556; //0.704;
			MODz[0]= 0.266844; //0.286;
			// GL PAS matrix
			MODx[1]= 0.175411;
			MODy[1]= 0.972178;
			MODz[1]= 0.154589;
			// Zero matrix
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("RGB")){
			//R
			MODx[0]= 0.0;
			MODy[0]= 1.0;
			MODz[0]= 1.0;
			//G
			MODx[1]= 1.0;
			MODy[1]= 0.0;
			MODz[1]= 1.0;
			//B
			MODx[2]= 1.0;
			MODy[2]= 1.0;
			MODz[2]= 0.0;
		}

		if (myStain.equals("CMY")){
			//C
			MODx[0]= 1.0;
			MODy[0]= 0.0;
			MODz[0]= 0.0;
			//M
			MODx[1]= 0.0;
			MODy[1]= 1.0;
			MODz[1]= 0.0;
			//Y
			MODx[2]= 0.0;
			MODy[2]= 0.0;
			MODz[2]= 1.0;
		}

		if (myStain.equals("User values")){
			GenericDialog gd2 = new GenericDialog("User values", IJ.getInstance());
			gd2.addMessage("Colour[1]");
			gd2.addNumericField("[R1]", 0, 5);
			gd2.addNumericField("[G1]", 0, 5);
			gd2.addNumericField("[B1]", 0, 5);
			gd2.addMessage("Colour[2]");
			gd2.addNumericField("[R2]", 0, 5);
			gd2.addNumericField("[G2]", 0, 5);
			gd2.addNumericField("[B2]", 0, 5);
			gd2.addMessage("Colour[3]");
			gd2.addNumericField("[R3]", 0, 5);
			gd2.addNumericField("[G3]", 0, 5);
			gd2.addNumericField("[B3]", 0, 5);

			gd2.showDialog();
			if (gd2.wasCanceled())
				return;

			MODx[0]= gd2.getNextNumber();
			MODy[0]= gd2.getNextNumber();
			MODz[0]= gd2.getNextNumber();
			MODx[1]= gd2.getNextNumber();
			MODy[1]= gd2.getNextNumber();
			MODz[1]= gd2.getNextNumber();
			MODx[2]= gd2.getNextNumber();
			MODy[2]= gd2.getNextNumber();
			MODz[2]= gd2.getNextNumber();
		}

		if (myStain.equals("From ROI")){
			IJ.runMacro("setOption('DisablePopupMenu', true)");
			// imp.getCanvas().disablePopupMenu(true);
			double [] rgbOD = new double[3];
			for (i=0; i<3; i++){
				getmeanRGBODfromROI(i, rgbOD, imp);
				MODx[i]= rgbOD[0];
				MODy[i]= rgbOD[1];
				MODz[i]= rgbOD[2];
			}

			IJ.runMacro("setOption('DisablePopupMenu', false)");
			// imp.getCanvas().disablePopupMenu(false);

		}
		// start
		for (i=0; i<3; i++){
			//normalise vector length
			cosx[i]=cosy[i]=cosz[i]=0.0;
			len[i]=Math.sqrt(MODx[i]*MODx[i] + MODy[i]*MODy[i] + MODz[i]*MODz[i]);
			if (len[i] != 0.0){
				cosx[i]= MODx[i]/len[i];
				cosy[i]= MODy[i]/len[i];
				cosz[i]= MODz[i]/len[i];
			}
		}


		// translation matrix
		if (cosx[1]==0.0){ //2nd colour is unspecified
			if (cosy[1]==0.0){
				if (cosz[1]==0.0){
					cosx[1]=cosz[0];
					cosy[1]=cosx[0];
					cosz[1]=cosy[0];
				}
			}
		}

		if (cosx[2]==0.0){ // 3rd colour is unspecified
			if (cosy[2]==0.0){
				if (cosz[2]==0.0){
					if ((cosx[0]*cosx[0] + cosx[1]*cosx[1])> 1){
						if (doIshow)
							IJ.log("Colour_3 has a negative R component.");
						cosx[2]=0.0;
					}
					else {
						cosx[2]=Math.sqrt(1.0-(cosx[0]*cosx[0])-(cosx[1]*cosx[1]));
					}

					if ((cosy[0]*cosy[0] + cosy[1]*cosy[1])> 1){
						if (doIshow)
							IJ.log("Colour_3 has a negative G component.");
						cosy[2]=0.0;
					}
					else {
						cosy[2]=Math.sqrt(1.0-(cosy[0]*cosy[0])-(cosy[1]*cosy[1]));
					}

					if ((cosz[0]*cosz[0] + cosz[1]*cosz[1])> 1){
						if (doIshow)
							IJ.log("Colour_3 has a negative B component.");
						cosz[2]=0.0;
					}
					else {
						cosz[2]=Math.sqrt(1.0-(cosz[0]*cosz[0])-(cosz[1]*cosz[1]));
					}
				}
			}
		}

		leng=Math.sqrt(cosx[2]*cosx[2] + cosy[2]*cosy[2] + cosz[2]*cosz[2]);

		cosx[2]= cosx[2]/leng;
		cosy[2]= cosy[2]/leng;
		cosz[2]= cosz[2]/leng;

		for (i=0; i<3; i++){
			if (cosx[i] == 0.0) cosx[i] = 0.001;
			if (cosy[i] == 0.0) cosy[i] = 0.001;
			if (cosz[i] == 0.0) cosz[i] = 0.001;
		}

		if (!hideLegend) {
			ImagePlus imp0 = NewImage.createRGBImage("Colour Deconvolution", 350, 65, 1, 0);
			ImageProcessor ip0 = imp0.getProcessor();
			ip0.setFont(new Font("Monospaced", Font.BOLD, 11));
			ip0.setAntialiasedText(true);
			ip0.setColor(Color.black);
			ip0.moveTo(10,15);
			ip0.drawString("Colour deconvolution: "+myStain);
			ip0.setFont(new Font("Monospaced", Font.PLAIN, 10));

			for (i=0; i<3; i++){
				ip0.setRoi(10,18+ i*15, 14, 14); 
				ip0.setColor( 
					(((255 -(int)(255.0* cosx[i])) & 0xff)<<16)+
					(((255 -(int)(255.0* cosy[i])) & 0xff)<<8 )+
					(((255 -(int)(255.0* cosz[i])) & 0xff)    ));
				ip0.fill();
				ip0.setFont(new Font("Monospaced", Font.PLAIN, 10));
				ip0.setAntialiasedText(true);
				ip0.setColor(Color.black);
				ip0.moveTo(27,32+ i*15);
				ip0.drawString("Colour_"+(i+1)+" R:"+(float)cosx[i]+", G:"+(float)cosy[i]+", B:"+(float)cosz[i] );
			}
			imp0.show();
			imp0.updateAndDraw();
		}

/*		if (myStain.equals("From ROI")){
			IJ.showMessage("Vectors","Colour[1]:\n"+
			"  R: "+(float) cosx[0] +"\n"+
			"  G: "+(float) cosy[0] +"\n"+
			"  B: "+(float) cosz[0] +"\n \n"+
			"Colour[2]:\n"+
			"  R: "+ (float) cosx[1] +"\n"+
			"  G: "+ (float) cosy[1] +"\n"+
			"  B: "+ (float) cosz[1] +"\n \n"+
			"Colour[3]:\n"+
			"  R: "+ (float) cosx[2] +"\n"+
			"  G: "+ (float) cosy[2] +"\n"+
			"  B: "+ (float) cosz[2]);
		}
*/

		if (doIshow){
			IJ.log( myStain +" Vector Matrix ---");
			for (i=0; i<3; i++){
				IJ.log("Colour["+(i+1)+"]:\n"+
				"  R"+(i+1)+": "+ (float) MODx[i] +"\n"+
				"  G"+(i+1)+": "+ (float) MODy[i] +"\n"+
				"  B"+(i+1)+": "+ (float) MODz[i] +"\n \n");
			}

/*			IJ.log( myStain +" Translation Matrix ---");
			for (i=0; i<3; i++){
				IJ.log("Colour["+(i+1)+"]:\n"+
				"  R"+(i+1)+": "+ (float) cosx[i] +"\n"+
				"  G"+(i+1)+": "+ (float) cosy[i] +"\n"+
				"  B"+(i+1)+": "+ (float) cosz[i] +"\n \n");
			}
*/
			IJ.log( myStain +" Java code ---");
			IJ.log("\t\tif (myStain.equals(\"New_Stain\")){");
			IJ.log("\t\t// This is the New_Stain");
			for (i=0; i<3; i++){
				IJ.log("\t\t\tMODx["+i+"]="+ (float) cosx[i] +";\n"+
					"\t\t\tMODy["+i+"]="+ (float) cosy[i] +";\n"+
					"\t\t\tMODz["+i+"]="+ (float) cosz[i] +";\n\n");
			}
			IJ.log("}");
		}

		//matrix inversion
		A = cosy[1] - cosx[1] * cosy[0] / cosx[0];
		V = cosz[1] - cosx[1] * cosz[0] / cosx[0];
		C = cosz[2] - cosy[2] * V/A + cosx[2] * (V/A * cosy[0] / cosx[0] - cosz[0] / cosx[0]);
		q[2] = (-cosx[2] / cosx[0] - cosx[2] / A * cosx[1] / cosx[0] * cosy[0] / cosx[0] + cosy[2] / A * cosx[1] / cosx[0]) / C;
		q[1] = -q[2] * V / A - cosx[1] / (cosx[0] * A);
		q[0] = 1.0 / cosx[0] - q[1] * cosy[0] / cosx[0] - q[2] * cosz[0] / cosx[0];
		q[5] = (-cosy[2] / A + cosx[2] / A * cosy[0] / cosx[0]) / C;
		q[4] = -q[5] * V / A + 1.0 / A;
		q[3] = -q[4] * cosy[0] / cosx[0] - q[5] * cosz[0] / cosx[0];
		q[8] = 1.0 / C;
		q[7] = -q[8] * V / A;
		q[6] = -q[7] * cosy[0] / cosx[0] - q[8] * cosz[0] / cosx[0];

                // initialize 3 output colour stacks
                ImageStack[] outputstack = new ImageStack[3];
		for (i=0; i<3; i++){
			for (j=0; j<256; j++) { //LUT[1]
				//if (cosx[i] < 0)
				//	rLUT[255-j]=(byte)(255.0 + (double)j * cosx[i]);
				//else
					rLUT[255-j]=(byte)(255.0 - (double)j * cosx[i]);

				//if (cosy[i] < 0)
				//	gLUT[255-j]=(byte)(255.0 + (double)j * cosy[i]);
				//else
					gLUT[255-j]=(byte)(255.0 - (double)j * cosy[i]);

				//if (cosz[i] < 0)
				//	bLUT[255-j]=(byte)(255.0 + (double)j * cosz[i]);
				///else
					bLUT[255-j]=(byte)(255.0 - (double)j * cosz[i]);
			}
			IndexColorModel cm = new IndexColorModel(8, 256, rLUT, gLUT, bLUT);
                        outputstack[i] = new ImageStack(width, height, cm);
                }

		// translate ------------------
		int imagesize = width * height;
		for (int imagenum=1; imagenum<=stack.getSize(); imagenum++) {
			int[] pixels = (int[])stack.getPixels(imagenum);
			String label = stack.getSliceLabel(imagenum);
			byte[][] newpixels = new byte[3][];
			newpixels[0] = new byte[imagesize];
			newpixels[1] = new byte[imagesize];
			newpixels[2] = new byte[imagesize];

			for (j=0;j<imagesize;j++){
				 // log transform the RGB data
				int R = (pixels[j] & 0xff0000)>>16;
				int G = (pixels[j] & 0x00ff00)>>8 ;
				int B = (pixels[j] & 0x0000ff);
				double Rlog = -((255.0*Math.log(((double)R+1)/255.0))/log255);
				double Glog = -((255.0*Math.log(((double)G+1)/255.0))/log255);
				double Blog = -((255.0*Math.log(((double)B+1)/255.0))/log255);
				for (i=0; i<3; i++){
                                    // rescale to match original paper values
                                    double Rscaled = Rlog * q[i*3];
                                    double Gscaled = Glog * q[i*3+1];
                                    double Bscaled = Blog * q[i*3+2];
                                    double output = Math.exp(-((Rscaled + Gscaled + Bscaled) - 255.0) * log255 / 255.0);
                                    if(output>255) {
                                        output=255;
                                    }
                                    newpixels[i][j]=(byte)(0xff&(int)(Math.floor(output+.5)));
				}
			}
			 // add new values to output images
			outputstack[0].addSlice(label,newpixels[0]);
			outputstack[1].addSlice(label,newpixels[1]);
			outputstack[2].addSlice(label,newpixels[2]);
		}
		new ImagePlus(title+"-(Colour_1)",outputstack[0]).show();
		new ImagePlus(title+"-(Colour_2)",outputstack[1]).show();
		new ImagePlus(title+"-(Colour_3)",outputstack[2]).show();
	}

	static void getmeanRGBODfromROI(int i, double [] rgbOD, ImagePlus imp){
		//get a ROI and its mean optical density. GL
		int [] xyzf = new int [4]; //[0]=x, [1]=y, [2]=z, [3]=flags
		int x1, y1, x2, y2, h=0, w=0, px=0, py=0, x, y,p;
		double log255=Math.log(255.0);
		ImageProcessor ip = imp.getProcessor();
		int mw = ip.getWidth()-1;
		int mh = ip.getHeight()-1;

		IJ.showMessage("Select ROI for Colour_"+(i+1)+".\n \n(Right-click to end)");
		getCursorLoc( xyzf, imp );
		while ((xyzf[3] & 4) !=0){  //trap until right released
			getCursorLoc( xyzf, imp );
			IJ.wait(20);
		}

		while (((xyzf[3] & 16) == 0) && ((xyzf[3] & 4) ==0)) { //trap until one is pressed
			getCursorLoc( xyzf, imp );
			IJ.wait(20);
		}

		rgbOD[0]=0;
		rgbOD[1]=0;
		rgbOD[2]=0;

		if ((xyzf[3] & 4) == 0){// right was not pressed, but left (ROI) was
			x1=xyzf[0];
			y1=xyzf[1];
			//IJ.write("first point x:" + x1 + "  y:" + y1);
			x2=x1;  y2=y1;
			while ((xyzf[3] & 4) == 0){//until right pressed
				getCursorLoc( xyzf, imp );
				if (xyzf[0]!=x2 || xyzf[1]!=y2) {
					if (xyzf[0]<0) {
                                        xyzf[0]=0;
                                    }
					if (xyzf[1]<0) {
                                        xyzf[1]=0;
                                    }
					if (xyzf[0]>mw) {
                                        xyzf[0]=mw;
                                    }
					if (xyzf[1]>mh) {
                                        xyzf[1]=mh;
                                    }
					x2=xyzf[0]; y2=xyzf[1];
					w=x2-x1+1;
					h=y2-y1+1;
					if (x2<x1) {px=x2;  w=(x1-x2)+1;} else {
                                        px=x1;
                                    }
					if (y2<y1) {py=y2;  h=(y1-y2)+1;} else {
                                        py=y1;
                                    }
					IJ.makeRectangle(px, py, w, h);
					//IJ.write("Box x:" + x2 +"  y:" + y2+" w:"+w+" h:"+h);
				}
				IJ.wait(20);
			}
			while ((xyzf[3] & 16) !=0){  //trap until left released
				getCursorLoc( xyzf, imp );
				IJ.wait(20);
			}

			for (x=px;x<(px+w);x++){
				for(y=py;y<(py+h);y++){
					p=ip.getPixel(x,y);
					// rescale to match original paper values
					rgbOD[0] = rgbOD[0]+ (-((255.0*Math.log(((double)((p & 0xff0000)>>16)+1)/255.0))/log255));
					rgbOD[1] = rgbOD[1]+ (-((255.0*Math.log(((double)((p & 0x00ff00)>> 8) +1)/255.0))/log255));
					rgbOD[2] = rgbOD[2]+ (-((255.0*Math.log(((double)((p & 0x0000ff))        +1)/255.0))/log255));
				}
			}
			rgbOD[0] = rgbOD[0] / (w*h);
			rgbOD[1] = rgbOD[1] / (w*h);
			rgbOD[2] = rgbOD[2] / (w*h);
		}
		IJ.run("Select None");
	}


	static void getCursorLoc(int [] xyzf, ImagePlus imp ) {
		ImageWindow win = imp.getWindow();
		ImageCanvas ic = win.getCanvas();
		Point p = ic.getCursorLoc();
		xyzf[0]=p.x;
		xyzf[1]=p.y;
		xyzf[2]=imp.getCurrentSlice()-1;
		xyzf[3]=ic.getModifiers();
	}
        
        /**
         * Lets the user indicate colors which should be used for the different channels.
         * @param se The parent TMARKER Frame on which messages are shown and the JXColorSelectionButton is shown.
         * @return A list with three colors. Each color is represented as double[] of size 3 (R G B).
         */
        static ArrayList<double[]> getUserColors(StainingEstimation se, TMAspot ts) {
            double [] MODx = new double[3];
	    double [] MODy = new double[3];
	    double [] MODz = new double[3];
            
            GenericDialog gd2 = new GenericDialog("User values", se);
            gd2.setLayout(new BoxLayout(gd2, BoxLayout.PAGE_AXIS));
            gd2.addMessage("Please select the colors to be separated.\nFor a two-channel separation, leave channel 3 white.");
            gd2.addMessage("Channel 1");
            Color c = (se!=null && se.getParam_UserColor1()!=null) ? se.getParam_UserColor1() : ts.getAverageColor(TMALabel.LABEL_NEG, false);
            JXColorSelectionButton c1 = new JXColorSelectionButton(c);
            gd2.add(c1);
            gd2.addMessage("Channel 2");
            c = (se!=null && se.getParam_UserColor2()!=null) ? se.getParam_UserColor2() : ts.getAverageColor(TMALabel.LABEL_POS, false);
            JXColorSelectionButton c2 = new JXColorSelectionButton(c);
            gd2.add(c2);
            gd2.addMessage("Channel 3");
            c = (se!=null && se.getParam_UserColor3()!=null) ? se.getParam_UserColor3() : Color.WHITE;
            JXColorSelectionButton c3 = new JXColorSelectionButton(c);
            gd2.add(c3);
            
            gd2.showDialog();
            if (gd2.wasCanceled()) {
                return null;
            }

            if (se!=null) {
                se.setParam_UserColor1(((JXColorSelectionButton) (gd2.getComponent(2))).getBackground());
                se.setParam_UserColor2(((JXColorSelectionButton) (gd2.getComponent(4))).getBackground());
                se.setParam_UserColor3(((JXColorSelectionButton) (gd2.getComponent(6))).getBackground());
            }
            
            MODx[0]= ((JXColorSelectionButton) (gd2.getComponent(2))).getBackground().getRed();
            MODy[0]= ((JXColorSelectionButton) (gd2.getComponent(2))).getBackground().getGreen();
            MODz[0]= ((JXColorSelectionButton) (gd2.getComponent(2))).getBackground().getBlue();
            MODx[1]= ((JXColorSelectionButton) (gd2.getComponent(4))).getBackground().getRed();
            MODy[1]= ((JXColorSelectionButton) (gd2.getComponent(4))).getBackground().getGreen();
            MODz[1]= ((JXColorSelectionButton) (gd2.getComponent(4))).getBackground().getBlue();
            MODx[2]= ((JXColorSelectionButton) (gd2.getComponent(6))).getBackground().getRed();
            MODy[2]= ((JXColorSelectionButton) (gd2.getComponent(6))).getBackground().getGreen();
            MODz[2]= ((JXColorSelectionButton) (gd2.getComponent(6))).getBackground().getBlue();
            
            ArrayList<double[]> uc = new ArrayList<>(3);
            uc.add(MODx);
            uc.add(MODy);
            uc.add(MODz);
            
            return uc;
	}

}

