//****************************************************************
// 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.addChoice("Distance Map", ImageList);
ImageList=Array.concat("*none*",ImageList);
Dialog.addChoice("Mask", ImageList);
Dialog.addString("exclude the following Mask Values from analysis (e.g. 0-5,22)","0",55);
Dialog.show();

FluorescenceImageOriginal=Dialog.getChoice();
DistanceMapOriginal=Dialog.getChoice();
MaskOriginal=Dialog.getChoice();
getDimensions(xmax,ymax,channels,slices,frames);
selectWindow(DistanceMapOriginal);
resetMinAndMax();
selectWindow(FluorescenceImageOriginal);
resetMinAndMax();

//****************************************************************
// 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,",");
  /*if (indexC==-1)
    {
    indexC=lengthOf(excludedValues)-1;
    }*/
    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)



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

FluorescenceImage="CopyFluorescence";
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="[Distribution]";
run("Table...", "name=[Distribution] width=250 height=600");
print(f, "\\Headings:slice\tdistance\tmean\tSD");

for(s=1;s<slices+1;s++)
  {
  setSlice(s);
  sum=newArray(256);
  n=newArray(256);
  mean=newArray(256);
  SD=newArray(256);
  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)
        {
        sum[distance]+=intensity;
        n[distance]++;
        }
      }
    }
  for(i=0;i<255;i++)
    {
    mean[i]=sum[i]/n[i];
    }
  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)
        {
        SD[distance]+=(intensity-mean[distance])*(intensity-mean[distance]);
        }
      }
    }


  for(i=0;i<255;i++)
    {
      if (n[i]>0)
        {
        SD[i]=sqrt(SD[i]/(n[i]-1));
        print(f,s+"\t"+i+"\t"+mean[i]+"\t"+SD[i]);
        }
    }
  } // of s
selectWindow("Combined Stacks");
run("Close");
setBatchMode(false);