//****************************************************************
// Create Dialog to select the appropriate Stacks
//****************************************************************
Dialog.create("DistanceCorrelation");
if(nImages==0)
  {
  Dialog.addMessage("No Images Open!");
  Dialog.show();
  exit();
  }
ImageList=newArray(nImages);
for (i=1; i<=nImages; i++) {
        selectImage(i);
        ImageList[i-1]=getTitle;
     }
Dialog.addChoice("Fluorescence Image",ImageList);
Dialog.addCheckbox("Calculate Diference?",false);
Dialog.addNumber("relativ position before ROIj",1);
Dialog.addNumber("relativ position after ROIj",1);
Dialog.addChoice("Distance Map", ImageList);
Dialog.addNumber("extend border by Pixels",5);
ImageList=Array.concat("*none*",ImageList);
Dialog.addChoice("Mask", ImageList);
Dialog.addString("exclude the following Mask Values from analysis (e.g. 0-5,22)","0",50);
Dialog.show();

FluorescenceImageOriginal=Dialog.getChoice();
DistanceMapOriginal=Dialog.getChoice();
MaskOriginal=Dialog.getChoice();
CalculateDifference=Dialog.getCheckbox();
PosBefore=Dialog.getNumber();
PosAfter=Dialog.getNumber();
ROIjWidth=Dialog.getNumber();
getDimensions(xmax,ymax,channels,slices,frames);


//****************************************************************
// Create Binary Mapping Function to excluded particular Mask Values from further analysis
//****************************************************************
valid=newArray(65536);  // By Default all values bigger than 0 in the Mask image are considered as valid 
for(i=0;i<65536;i++)
  {
  valid[i]=true;
  }
excludedValues=Dialog.getString()+",";
do
  {
  indexC=indexOf(excludedValues,",");
  subStringC=substring(excludedValues,0,indexC);
  excludedValues=substring(excludedValues,indexC+1);
  indexR=indexOf(subStringC,"-");
  if (indexR==-1)
    {
    low=parseInt(subStringC);
    high=low;
    }
  else
    {
    low=parseInt(substring(subStringC,0,indexR));
    high=parseInt(substring(subStringC,indexR+1,indexC));
    }
  for(i=low;i<=high;i++)
    {
    valid[i]=false;
    }
  } while (lengthOf(excludedValues)>0)

//******************************************************************
/*print("test");
for(i=0;i<65536;i++)
  {
  if (!valid[i])
    {
    print("not Valid "+i);
    }
  }*/

//****************************************************************
// Make Copy of the stacks and combine them for faster analysis
//****************************************************************
//setBatchMode(true);

FluorescenceImage="CopyFluorescence";
selectImage(FluorescenceImageOriginal);
if (CalculateDifference)
  {
print(" "+slices-PosBefore-PosAfter);
  run("Make Substack...", "  slices="+PosBefore+PosAfter+2+"-"+slices);
  rename(FluorescenceImage);
  setSlice(nSlices);
  for(i=0;i<PosBefore+PosAfter+1;i++)
    {  
    run("Add Slice");
    }
  imageCalculator("difference stack",FluorescenceImage,FluorescenceImageOriginal);
  }
else
  {
  newImage(FluorescenceImage,"16-bit",xmax,ymax,slices);
  imageCalculator("copy stack",FluorescenceImage,FluorescenceImageOriginal);
  }
//selectImage(FluorescenceImage);
run("16-bit");

DistanceMap="CopyDistance";
newImage(DistanceMap,"16-bit",xmax,ymax,slices);
imageCalculator("copy stack",DistanceMap,DistanceMapOriginal);
//selectImage(DistanceMap);
run("16-bit");

Mask="CopyMask";
newImage(Mask,"16-bit",xmax,ymax,slices);
if(MaskOriginal=="*none*")
  {
  run("Add...", "value=65535 stack");  
  }
else
  {
  imageCalculator("copy stack",Mask,MaskOriginal);
  }
//selectImage(Mask);
run("16-bit");

//slices=3;

for(s=1;s<slices+1;s++)
  {
  setSlice(s);
  for(x=0;x<xmax;x++)
    {
    for(y=0;y<ymax;y++)
      {
      if (valid[getPixel(x,y)]) 
        {
        setPixel(x,y,65535);
        }
      else
        {  
        setPixel(x,y,0);
        }
      }
    }
  } // of s

run("Combine...","stack1=["+FluorescenceImage+"] stack2=["+DistanceMap+"]");
run("Combine...","stack1=[Combined Stacks] stack2=["+Mask+"]");

//****************************************************************
// Analyse the Combined stacks
//****************************************************************
//TODO: Chekcen if Table exists

f="[Intensity]";
run("Table...", "name="+f+" width=250 height=600");
print(f, "\\Headings:slice\tmean\tSD");

for(s=1;s<slices+1;s++)
  {
  setSlice(s);
  sum=0;
  n=0;
  
  for(y=0;y<ymax;y++)
    {
    for(x=0;x<xmax;x++)
      {
      intensity=getPixel(x,y);
      distance=getPixel(x+xmax,y);
      mask=getPixel(x+2*xmax,y);
      if ((mask!=0)&&(distance<=ROIjWidth))
        {
        sum+=intensity;
        n++;
        }
      }
    }
  mean=sum/n;
  SD=0;
  for(y=0;y<ymax;y++)
    {
    for(x=0;x<xmax;x++)
      {
      intensity=getPixel(x,y);
      distance=getPixel(x+xmax,y);
      mask=getPixel(x+2*xmax,y);
      if ((mask!=0)&&(distance<=ROIjWidth))
        {
        SD+=(intensity-mean)*(intensity-mean);
        }
      }
    }
    SD=sqrt(SD/(n-1));
    print(f,s+"\t"+mean+"\t"+SD);
  } // of s
selectWindow("Combined Stacks");
//run("Close");
setBatchMode(false);