//Make a list of open images
function getOpenima(id){
openim=newArray(nImages);
if(nImages==0) exit("no image open");
for(i=1;i<=nImages;i++) { 
	selectImage(i);
          if(id==true) openim[i-1]=getImageID();
          else openim[i-1]=getTitle;
      } 
return openim;
}



//Close windows except those spécified
function closeWindowexept(openim,exept){
ima=getOpenima(true);
for(i=0;i<ima.length;i++){
	closed=1;
	for(j=0;j<exept.length;j++){
			if(ima[i]==exept[j]) closed=0;
		}
	for(k=0;k<openim.length;k++){
			if(ima[i]==openim[k]) closed=0;
		}
	if(closed==1){
		selectImage(ima[i]);
		close();
		}
	}
}


//Create a dialogue box to choose the images and parametres used for treatement 
function getChoiceim(){
list=getOpenima(false);
Dialog.create("choix des images");
Dialog.addChoice("mask",list);
Dialog.addChoice("image à traiter",list);
Dialog.addNumber("coefficient pour le traitement", 0);
Dialog.addNumber("rayon pour le filtre médian",0);
Dialog.addChoice("traitement continuité", newArray("non","oui"));
Dialog.show;
cmask=Dialog.getChoice(); 
selectImage(cmask);
imask=getImageID();
cima=Dialog.getChoice(); 
selectImage(cima);
iima=getImageID();
coeff=Dialog.getNumber();
med=Dialog.getNumber();
cont=Dialog.getChoice(); 
if(cont=="oui") conti=1;
else conti=0;
choice=newArray(imask,iima,coeff,med,conti);
return choice;
}


//Evaluate type of the image to be treated 
function getType(choice){
type=newArray(2);
selectImage(choice[0]);
type[0]=bitDepth();
selectImage(choice[1]);
type[1]=bitDepth();
return type;
}


//Show an error message if image to be treated are not in te corrected type
function showErrorchoice(type){
if(type[0]!=8){
	exit("mask (black background) image must be 8-bit type"); 
	}
if(type[1]!=16){
	exit("image to be treated must be 16-bit type"); 
	}
}


//Calculate minimal and maximal grey values in a stack serie
function getstackMinandMax(){
stackMinandMax=newArray(65535,0,0,0);
for(i=0;i<nSlices;i++){
	setSlice(i+1);
	getRawStatistics(npix,mean,min,max,std,hist);
	if (min<stackMinandMax[0]) {
		stackMinandMax[0]=min;
		stackMinandMax[2]=i+1;
		}
	if (max>stackMinandMax[1]) {
		stackMinandMax[1]=max;
		stackMinandMax[3]=i+1;
		}
	}
return stackMinandMax;
}


//Apply by subtraction the mask to the image to be treated 
function applyMask(mask,im,openim){
selectImage(mask);
run("Duplicate...", "title=inter-mask");
inter=getImageID();
run("16-bit");
run("Invert");
run("Multiply...", "value=16");
selectImage(im);
if(nSlices!=1){ 
	m=getstackMinandMax();
	setSlice(m[3]);
	}
run("Duplicate...", "title=slicecalculator");
slicecalc=getImageID(); 
imageCalculator("Subtract create",slicecalc,inter); 
rename("result");
run("32-bit"); 
setThreshold(1,4095);
run("NaN Background");
slicecalcres=getImageID();
exept=newArray(1); exept[0]=slicecalcres;
closeWindowexept(openim, exept);
return slicecalcres;
}

//Calculate cumulate sum of an array  
function getnPix(hist){
npix=0;
for(i=0;i<hist.length;i++){
npix=npix+hist[i];
}
return npix;
}


//calculate median value of an array (count: cumulate sum of the array)
function getMedian(hist, count) {
      n = hist.length;
      sum = 0;
      i = -1;
      count2 = count/2;
      do {
          sum += hist[++i];
      } while (sum<=count2);
      return i;
  }


//Calculate the window of intensities used for the image treatment (based on histogramm analysis and choosen coefficent)
function getWindowanalyse(slicecalc,coeff){
selectImage(slicecalc);
getRawStatistics(n,meanc,minc,maxc,stdc);
getHistogram(values,histc,4096,0,4096);
npixc=getnPix(histc);	
medianc=getMedian(histc, npixc);
newmin1=medianc+(coeff-2)*stdc; if(newmin1<minc) newmin1=minc;
newmax1=1.5*maxc; if(newmax1>4095) newmax1=4095;
newmin2=medianc+coeff*stdc; if(newmin2<minc) newmin2=minc;	
newmax2=1.2*maxc; if(newmax2>4095) newmax2=4095;
selectImage(slicecalc); close();
windowa=newArray(newmin1,newmax1,newmin2,newmax2);
return windowa;
}


//Apply image processing
function applyTraitement(im,windowa,coeff,med,conti){
selectImage(im);
run("Duplicate...", "title=inter-1 duplicate");
im1=getImageID();
run("Duplicate...", "title=inter-2 duplicate");
im2=getImageID();

selectImage(im1);
setSlice(1);
setMinAndMax(windowa[0],windowa[1]);
run("8-bit");
run("Convolve...", "text1=[-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 24 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ] normalize stack");
if(conti!=0) {setMinAndMax(0,128);}
if(med!=0){run("Median...", "radius="+med+" stack");}
	
selectImage(im2);
setMinAndMax(windowa[2],windowa[3]);
run("8-bit");

imageCalculator("Add create stack", im1,im2);
selectImage(im1); close();
selectImage(im2);close();
rename("traitement_coefficient"+coeff+"_median"+med+"_continuité"+conti);
}


macro "image processing"{
openim=getOpenima(true);
choice=getChoiceim();
showErrorchoice(getType(choice));
slicecal=applyMask(choice[0],choice[1],openim);
windowa=getWindowanalyse(slicecal,choice[2]);
applyTraitement(choice[1],windowa,choice[2],choice[3],choice[4]);
run("Tile");

}

macro"help"{
message1= "Image informations:\n \t-Image must be 16-bit type \n \t-Mask is binary black background image to specify pixels that belong to the area (cell) to be treated";
message2="Processing informations:\n \t -Background coefficient is required to specify the histogramm window to be analysed depending on object occupancy within your area: \n \t \t For large occupancy coefficient should be small (-2-0) whereas for small occupancy it should be higher (1-4)" ";
message3="\t - Median filter radius must be specificied: for small object choose small radius (1) whereas for big objects higher radius (2) is recommended";
message4= "\t -In some case some object discontinuity can appear after processing: to avoid this you can use the continuity option";
message=message1+"\n"+message2+ "\n"+ message3+ "\n"+ message4;
showMessage("help",message);
}