// Calculate_GP_02.txt
// https://sils.fnwi.uva.nl/bcb/objectj/examples/CalculateGP/MD/gp.html
// N. Vischer 10-jun-2016   23:55



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() {
	reds = newArray(255);
	greens = newArray(255);
	blues = newArray(255);
	for (jj = 1; jj < 128; jj++) {//start at 1 to keep endpoints black
		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];
	}            
	setLut(reds, greens, blues);
 }