// Calculate_GP_02.txt
// http://www.fnwi.uva.nl/objectj/examples/CalculateGP/
// N. Vischer 10-jun-2016



macro "Calculate GP [F1]"{
	backgrounds = newArray(102, 103);//will be subtracted from images
	significants = newArray(10, 12);//smaller intensities will be invalid

	if (nImages < 2)
		exit("Exactly 2 images must be open: 1=ImageA, 2=ImageB");

	titles = newArray(2);
	selectImage(1);
	titles[0] = getTitle;
	selectImage(2);
	titles[1] = getTitle;
	


	if (nImages > 2){
		showMessageWithCancel("Delete all images except these two:\n \n" + titles[0] + "\n" + titles[1]);
	for (img = nImages; img >2; img--){
		selectImage(img);
		close;
	}
	Dialog.create("Enter values");
   	Dialog.addNumber("Subtract Backgnd A:", backgrounds[0]);
   	Dialog.addNumber("Subtract Backgnd B:", backgrounds[1]);
  	Dialog.addNumber("Significance A:", significants[0]);
   	Dialog.addNumber("Significance B:", significants[1]);
	Dialog.show();
  	backgrounds[0] = Dialog.getNumber();
  	backgrounds[1] = Dialog.getNumber();
  	significants[0] = Dialog.getNumber();
  	significants[1] = Dialog.getNumber();
 


	selectImage(1);
	run("Duplicate...", "title=A");
	selectImage(2);
	run("Duplicate...", "title=B");

	selectImage("A");
	run("Subtract...", "value=" + backgrounds[0]);
	changeValues(0, significants[0], 0);	
	run("Enhance Contrast", "saturated=0.35");
	
	selectImage("B");
	run("Subtract...", "value=" + backgrounds[1]);
	changeValues(0, significants[1], 0);	
	run("Enhance Contrast", "saturated=0.35");
	

	imageCalculator("Subtract create 32-bit", "A", "B");
	rename("A-B");
	imageCalculator("Add create 32-bit", "A", "B");
	rename("A+B");
	imageCalculator("Divide create 32-bit", "A-B","A+B");
	rename("Ratio");
	changeValues(-1, -1, NaN);//make extrema invalid	
	changeValues(1, 1, NaN);	
	selectImage("Ratio");
	setLutRedWhiteBlue();
	setMinAndMax(-1, 1);		
}

macro 	"Set LUT  example"{
	setLutRedWhiteBlue();
	setMinAndMax(-0.4, 0.8);
}



function setLutRedWhiteBlue() {
	getLut(reds, greens, blues);
	for (jj = 0; jj < 128; jj++) {
		blues[jj] = round(64.0 + jj / 128.0 * (255 - 64));
		greens[jj] = blues[jj];
		reds[jj] =  255;

		blues[255 - jj] = reds[jj];
		reds[255 - jj] = blues[jj];
		greens[255 - jj] = greens[jj];
	}
	//set either end to black:
	reds[0] = 0;
	greens[0] = 0;
	blues[0] = 0;
	reds[255] = 0;
	greens[255] = 0;
	blues[255] = 0;
            
	setLut(reds, greens, blues);
 }