/*JACoP: "Just Another Colocalization Plugin..." v1, 13/02/06
Fabrice P Cordelières, fabrice.cordelieres at curie.u-psud.fr
Susanne Bolte, Susanne.bolte@isv.cnrs-gif.fr
Copyright (C) 2006 Susanne Bolte & Fabrice P. Cordelières
License:
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*
*
*/
import ij.*;
import ij.ImagePlus.*;
import ij.plugin.*;
import ij.plugin.filter.*;
import ij.process.*;
import ij.gui.*;
import ij.measure.*;
import ij.util.*;
import java.awt.*;
import java.awt.event.*;
import java.util.*;
public class JACoP_ implements PlugIn, AdjustmentListener, TextListener {
//Load preferences
private static boolean PearsonBool=Prefs.get("JACoP_Pearson.boolean", true);
private static boolean OverlapBool=Prefs.get("JACoP_Overlap.boolean", true);
private static boolean MMBool=Prefs.get("JACoP_MM.boolean", true);
private static boolean CostesThrBool=Prefs.get("JACoP_CostesThr.boolean", true);
private static boolean CCFBool=Prefs.get("JACoP_CCF.boolean", true);
private static int CCFx=(int)Prefs.get("JACoP_CCFx.double", 20);
private static boolean CytoBool=Prefs.get("JACoP_Cyto.boolean", true);
private static boolean ICABool=Prefs.get("JACoP_ICA.boolean", true);
private static boolean CostesRandBool=Prefs.get("JACoP_CostesRand.boolean", true);
private static boolean DistBool=Prefs.get("JACoP_Dist.boolean", false);
private static boolean SpaPearBool=Prefs.get("JACoP_SpaPear.boolean", false);
private static boolean lineBool=Prefs.get("JACoP_Line.boolean", true);
private static int MicroscopeType=(int) Prefs.get("JACoP_MicroscopeType.double",0);
private static int xycal=(int) Prefs.get("JACoP_xycal.double",67);
private static int zcal=(int) Prefs.get("JACoP_zcal.double",200);
private static int wa=(int) Prefs.get("JACoP_wa.double",519);
private static int wb=(int) Prefs.get("JACoP_wb.double",565);
private static double NA=Prefs.get("JACoP_NA.double",1.4);
private static double IR=Prefs.get("JACoP_IR.double",1.518);
private static int xyBlock=(int) Prefs.get("JACoP_xyBlock.double",3);
private static int zBlock=(int) Prefs.get("JACoP_zBlock.double",3);
private static int nbRand=(int) Prefs.get("JACoP_nbRand.double", 200);
private static double binWidth=Prefs.get("JACoP_binWidth.double",0.0005);
private static int fillMeth=(int) Prefs.get("JACoP_fillMeth.double", 0);
private static boolean xyRand=Prefs.get("JACoP_xyRand.boolean", true);
private static boolean zRand=Prefs.get("JACoP_zRand.boolean", true);
private static boolean showRand=Prefs.get("JACoP_showRand.boolean", true);
String[] fill={"Shrink", "Pad w/black px"};
int widthCostes;
int heightCostes;
int nbsliceCostes;
//Variables
int [] imgIDList;
String [] imgTitleList;
String title;
ImagePlus imgA;
String titleA;
int depthA;
int widthA;
int heightA;
int nbsliceA;
double[] A;
int thrA;
double Amean=0;
double Amin;
double Amax=0;
double[] ACostes;
ImagePlus imgB;
String titleB;
int depthB;
int widthB;
int heightB;
int nbsliceB;
double[] B;
int thrB;
double Bmean=0;
double Bmin;
double Bmax=0;
double[] BCostes;
double[] BRandCostes;
Vector sliders;
Vector value;
//Values for stats
boolean doThat;
double sumA;
double sumB;
double sumAB;
double sumsqrA;
double Aarraymean;
double Barraymean;
//Parameters
String [] MicroscopeTypeList={"WideField","Confocal"};
String [] Chromophore={"Other","Alexa 350","DAPI","eCFP","eGFP","FITC","Alexa 488","YFP","CY3","Alexa 555","Alexa 546","EtBr","DsRed","Alexa 633","Alexa 647","CY5"};
int[] Excitation={0,346,358,436,488,494,495,514,548,555,556,518,558,632,650,650};
int[] Emission={0,442,461,474,507,518,519,527,562,565,573,605,583,647,668,670};
double resxy;
int pixxy;
double resz;
int pixz;
//int CCFy;
//int CCFz;
int a;
int b;
int c;
int d;
int e;
int f;
int i;
int j;
int k;
int l;
int m;
public void run(String arg) {
//Check that at least two images are opened
if (WindowManager.getImageCount()<2){
IJ.showMessage("Error", "Man,\n"+"You're in deep troubles:\n"+"you need at least two images...");
return;
}
//Get the list of currently opened windows
imgIDList= new int[WindowManager.getImageCount()];
imgIDList=WindowManager.getIDList();
imgTitleList=new String [WindowManager.getImageCount()];
for (i=0;i0) sumAcoloc+=A[i];
if (B[i]>thrB) sumAcolocThr+=A[i];
if (A[i]>0) sumBcoloc+=B[i];
if (A[i]>thrA) sumBcolocThr+=B[i];
sumA+=A[i];
sumB+=B[i];
}
double M1=sumAcoloc/sumA;
double M1Thr=sumAcolocThr/sumA;
double M2=sumBcoloc/sumB;
double M2Thr=sumBcolocThr/sumB;
IJ.log("\nManders' Coefficients (original):\nM1="+round(M1,3)+" (fraction of A overlapping B)\nM2="+round(M2,3)+" (fraction of B overlapping A)");
IJ.log("\nManders' Coefficients (using thresholds):\nM1="+round(M1Thr,3)+" (fraction of A overlapping B)\nM2="+round(M2Thr,3)+" (fraction of B overlapping A)");
}
public void CostesThr() {
int CostesThrA=(int)Amax;
int CostesThrB=(int)Bmax;
double CostesSumAThr=0;
double CostesSumA=0;
double CostesSumBThr=0;
double CostesSumB=0;
double CostesPearson=1;
double [] rx= new double[(int)(Amax-Amin+1)];
double [] ry= new double[(int)(Amax-Amin+1)];
double rmax=0;
double rmin=1;
doThat=true;
int count=0;
//First Step: define line equation
doThat=true;
double[] tmp=linreg(A,B,0,0);
double a=tmp[0];
double b=tmp[1];
double CoeffCorr=tmp[2];
doThat=false;
int LoopMin= (int) Math.max(Amin, (Bmin-b)/a);
int LoopMax= (int) Math.min(Amax, (Bmax-b)/a);
//Minimize r of points below (thrA,a*thrA+b)
for (i=LoopMax;i>=LoopMin;i--){
IJ.showStatus("Costes' threshold calculation in progress : "+(int)(100*(LoopMax-i)/(LoopMax-LoopMin))+"% done");
IJ.showProgress((int)LoopMax-i, (int)(LoopMax-LoopMin));
if (IJ.escapePressed()) {
IJ.showStatus("Task canceled by user");
IJ.showProgress(2,1);
return;
}
CostesPearson=linregCostes(A,B,i,a*i+b)[2];
rx[count]=i;
ry[count]=CostesPearson;
rmax=Math.max(rmax,CostesPearson);
rmin=Math.min(rmin,CostesPearson);
count++;
if (CostesPearson<=0){
CostesThrA=i;
CostesThrB=(int)(a*i+b);
i=(int)Amin-1;
}
}
for (i=0; iCostesThrA) CostesSumAThr+=A[i];
CostesSumB+=B[i];
if (B[i]>CostesThrB) CostesSumBThr+=B[i];
}
Plot plot=new Plot("Costes' threshold "+titleA+" and "+titleB,"ThrA", "Pearson's coefficient below",rx,ry);
plot.setLimits(LoopMax, CostesThrA, CostesPearson, rmax);
plot.setColor(Color.black);
plot.draw();
if (lineBool){
double[] xline={LoopMax, CostesThrA};
double[] yline={0, 0};
plot.setColor(Color.red);
plot.addPoints(xline, yline, 2);
}
plot.show();
ImagePlus CostesMask=NewImage.createRGBImage("Costes' mask",widthA,heightA,nbsliceA,0);
CostesMask.getProcessor().setValue(Math.pow(2, depthA));
for (k=1; k<=nbsliceA; k++){
CostesMask.setSlice(k);
for (j=0; jCostesThrA && color[1]>CostesThrB){
//CostesMask.getProcessor().setValue(((A[position]-CostesThrA)/(LoopMax-CostesThrA))*Math.pow(2, depthA));
//CostesMask.getProcessor().drawPixel(i,j);
for (l=0; l<=2; l++) color[l]=255;
}
CostesMask.getProcessor().putPixel(i,j,color);
}
}
}
CostesMask.setSlice(1);
CostesMask.show();
//IJ.setMinAndMax(0,Math.pow(2, depthA));
//IJ.run("Invert LUT");
IJ.showStatus("");
IJ.showProgress(2,1);
IJ.log("\nCostes' automatic threshold set to "+CostesThrA+" for imgA & "+CostesThrB+" for imgB");
IJ.log("Pearson's Coefficient:\nr="+round(linreg(A,B,CostesThrA,CostesThrB)[2],3)+" ("+round(CostesPearson,3)+" below thresholds)");
IJ.log("M1="+round(CostesSumAThr/CostesSumA,3)+" & M2="+round(CostesSumBThr/CostesSumB,3));
}
public void CCF(){
double num;
double den1;
double den2;
double CCF0=0;
double CCFmin=0;
int lmin=-CCFx;
double CCFmax=0;
int lmax=-CCFx;
double [] CCFarray=new double[2*CCFx+1];
double [] x=new double[2*CCFx+1];
int count=0;
IJ.log("\nVan Steensel's Cross-correlation Coefficient between "+titleA+" and "+titleB+":");
for (l=-CCFx; l<=CCFx; l++){
IJ.showStatus("CCF calculation in progress: "+(count+1)+"/"+(2*CCFx+1));
IJ.showProgress(count+1, 2*CCFx+1);
if (IJ.escapePressed()) {
IJ.showStatus("Task canceled by user");
IJ.showProgress(2,1);
return;
}
num=0;
den1=0;
den2=0;
for (k=1; k<=nbsliceA; k++){
for (j=0; j0 && i+lCCFmax){
CCFmax=CCF;
lmax=l;
}
}
x[count]=l;
CCFarray[count]=CCF;
count++;
}
IJ.log ("CCF min.: "+round(CCFmin,3)+" (obtained for dx="+lmin+") CCF max.: "+round(CCFmax,3)+" (obtained for dx="+lmax+")");
Plot plot=new Plot("Van Steensel's CCF between "+titleA+" and "+titleB,"dx", "CCF",x,CCFarray);
plot.setLimits(-CCFx, CCFx, CCFmin, CCFmax);
plot.setColor(Color.black);
plot.draw();
if (lineBool){
double[] xline={0,0};
double[] yline={CCFmin,CCFmax};
plot.setColor(Color.red);
plot.addPoints(xline, yline, 2);
}
IJ.showStatus("");
IJ.showProgress(2,1);
plot.show();
}
public void CytoFluo(){
//ImagePlus cyto=NewImage.createRGBImage("Cytofluo. "+titleB+"=f("+titleA+")", (int) Math.pow(2,depthA), (int) Math.pow(2,depthA), 1, 1);
Plot plot = new Plot("Cytofluorogram between "+titleA+" and "+titleB, titleA, titleB, A, B);
double limHigh=Math.max(Amax, Bmax);
double limLow=Math.min(Amin, Bmin);
plot.setLimits(limLow, limHigh, limLow, limHigh);
plot.setColor(Color.white);
doThat=true;
double[] tmp=linreg(A,B,0,0);
double a=tmp[0];
double b=tmp[1];
double CoeffCorr=tmp[2];
plot.draw();
plot.setColor(Color.black);
plot.addPoints(A, B, 6);
if (lineBool){
double[] xline={limLow,limHigh};
double[] yline={a*limLow+b,a*limHigh+b};
plot.setColor(Color.red);
plot.addPoints(xline, yline, 2);
}
//cyto.show();
plot.show();
IJ.log("\nCytofluorogram's parameters:\na: "+round(a,3)+"\nb: "+round(b,3)+"\nCorrelation coefficient: "+round(CoeffCorr,3));
}
public void ICA(){
double[] array=new double[1];
double[] Anorm=new double[A.length];
double[] Bnorm=new double[A.length];
double AnormMean=0;
double BnormMean=0;
double prodMin=0;
double prodMax=0;
double lim=0;
double[] x= new double[A.length];
double ICQ=0;
//Intensities are normalized to range from 0 to 1
for (i=0; iprodMax) prodMax=x[i];
if (x[i]0) ICQ++;
}
if (Math.abs(prodMin)>Math.abs(prodMax)){
lim=Math.abs(prodMin);
}else{
lim=Math.abs(prodMax);
}
ICQ=ICQ/A.length-0.5;
Plot plotA = new Plot("ICA A ("+titleA+")", "(Ai-a)(Bi-b)", titleA, x, Anorm);
plotA.setColor(Color.white);
plotA.setLimits(-lim, lim, 0, 1);
plotA.draw();
plotA.setColor(Color.black);
plotA.addPoints(x, Anorm, 6);
if (lineBool){
double[] xline={0,0};
double[] yline={0,1};
plotA.setColor(Color.red);
plotA.addPoints(xline, yline, 2);
}
plotA.show();
Plot plotB = new Plot("ICA B ("+titleB+")", "(Ai-a)(Bi-b)", titleB, x, Bnorm);
plotB.setColor(Color.white);
plotB.setLimits(-lim, lim, 0, 1);
plotB.draw();
plotB.setColor(Color.black);
plotB.addPoints(x, Bnorm, 6);
if (lineBool){
double[] xline={0,0};
double[] yline={0,1};
plotB.setColor(Color.red);
plotB.addPoints(xline, yline, 2);
}
plotB.show();
IJ.log("\nLi's Intensity correlation coefficient:\nICQ: "+ICQ);
}
public void CostesRand(){
double direction;
int shift;
int newposition;
if (xyRand || nbsliceCostes==1){
//If slices independent 2D there is no need to take into account the z thickness and ranndomization along z axis should not be done
zBlock=1;
zRand=false;
}
doThat=true;
double r2test=linreg(ACostes, BCostes, 0, 0)[2];
doThat=false;
double[] arrayR= new double[nbRand];
double mean=0;
double SD=0;
double Pval=0;
double[] arrayDistribR= new double[(int)(2/binWidth+1)];
for (f=0; f=widthCostes) newposition-=widthCostes;
if (newposition<0) newposition+=widthCostes;
BRandCostes[offsetCostes(newposition,b,c)]=BCostes[offsetCostes(a,b,c)];
}
}
}
}
}
for (i=0; i=heightCostes) newposition-=heightCostes;
if (newposition<0) newposition+=heightCostes;
BRandCostes[offsetCostes(b,newposition,c)]=BCostes[offsetCostes(b,a,c)];
}
}
}
}
}
for (i=0; inbsliceCostes) newposition-=nbsliceCostes;
if (newposition<1) newposition+=nbsliceCostes;
BRandCostes[offsetCostes(b,c,newposition)]=BCostes[offsetCostes(b,c,a)];
}
}
}
}
}
for (i=0; ithrA){
for (c=k-pixz; c<=k+pixz; c++){
for (b=j-pixxy; b<=j+pixxy; b++){
for (a=i-pixxy; athrB){
if (!alreadydone){
DistA.setSlice(k);
DistA.getProcessor().setValue(A[offset(i,j,k)]);
DistA.getProcessor().drawPixel(i,j);
InterA.setSlice(k);
InterA.getProcessor().setValue(A[offset(i,j,k)]);
InterA.getProcessor().drawPixel(i,j);
alreadydone=true;
}
if (i==a && j==b && k==c){
InterB.setSlice(k);
InterB.getProcessor().setValue(B[offset(i,j,k)]);
InterB.getProcessor().drawPixel(i,j);
}
DistB.setSlice(c);
DistB.getProcessor().setValue(B[offset(a,b,c)]);
DistB.getProcessor().drawPixel(a,b);
}
}
}
}
}
}
}
}
for (k=1; k<=nbsliceA; k++){
for (j=0; jthrA) thrpixA++;
if (B[offset(i,j,k)]>thrB) thrpixB++;
DistA.setSlice(k);
DistB.setSlice(k);
if (DistA.getProcessor().getPixel(i,j)!=0) countA++;
if (DistB.getProcessor().getPixel(i,j)!=0) countB++;
}
}
}
DistA.setSlice(1);
DistA.show();
IJ.setMinAndMax(imgA.getProcessor().getMin(),imgA.getProcessor().getMax());
IJ.run("Invert LUT");
DistB.setSlice(1);
DistB.show();
IJ.setMinAndMax(imgB.getProcessor().getMin(),imgB.getProcessor().getMax());
IJ.run("Invert LUT");
InterA.setSlice(1);
InterA.show();
IJ.setMinAndMax(imgA.getProcessor().getMin(),imgA.getProcessor().getMax());
IJ.run("Invert LUT");
InterB.setSlice(1);
InterB.show();
IJ.setMinAndMax(imgB.getProcessor().getMin(),imgB.getProcessor().getMax());
IJ.run("Invert LUT");
IJ.log("\nDistance based colocalization:\n% of positive A thresholded pixels="+100*countA/thrpixA+"\n% of positive B thresholded pixels="+100*countB/thrpixB);
}
public void SpatialPearson() {
ImagePlus SpaPear=NewImage.createFloatImage("Spatial Pearson of "+titleA+" & "+titleB,widthA,heightA,nbsliceA,0);
ImageProcessor ip;
double locMeanA=0;
double locMeanB=0;
double num=0;
double den1=0;
double den2=0;
for (k=1; k<=nbsliceA; k++){
SpaPear.setSlice(k);
ip=SpaPear.getProcessor();
for (j=0; jnbsliceA){
((Scrollbar)sliders.elementAt(2)).setValue(nbsliceA);
((TextField)value.elementAt(2)).setText(""+nbsliceA);
}
if ((int) Tools.parseDouble(((TextField)value.elementAt(2)).getText())<1){
((Scrollbar)sliders.elementAt(2)).setValue(1);
((TextField)value.elementAt(2)).setText("1");
}
thrA=((Scrollbar)sliders.elementAt(0)).getValue();
imgA.getProcessor().setThreshold(thrA,Math.pow(2,16),ImageProcessor.RED_LUT);
imgA.updateAndDraw();
thrB=((Scrollbar)sliders.elementAt(1)).getValue();
imgB.getProcessor().setThreshold(thrB,Math.pow(2,16),ImageProcessor.RED_LUT);
imgB.updateAndDraw();
imgA.setSlice(((Scrollbar)sliders.elementAt(2)).getValue());
imgB.setSlice(((Scrollbar)sliders.elementAt(2)).getValue());
}
public double round(double y, int z){
//Special tip to round numbers to 10^-2
y*=Math.pow(10,z);
y=(int) y;
y/=Math.pow(10,z);
return y;
}
public int offset(int m,int n,int o){
if (m+n*widthA+(o-1)*widthA*heightA>=widthA*heightA*nbsliceA){
return widthA*heightA*nbsliceA-1;
}else{
if (m+n*widthA+(o-1)*widthA*heightA<0){
return 0;
}else{
return m+n*widthA+(o-1)*widthA*heightA;
}
}
}
public int offsetCostes(int m,int n,int o){
if (m+n*widthCostes+(o-1)*widthCostes*heightCostes>=widthCostes*heightCostes*nbsliceCostes){
return widthCostes*heightCostes*nbsliceCostes-1;
}else{
if (m+n*widthCostes+(o-1)*widthCostes*heightCostes<0){
return 0;
}else{
return m+n*widthCostes+(o-1)*widthCostes*heightCostes;
}
}
}
public double[] linreg(double[] Aarray, double[] Barray, double TA, double TB){
double num=0;
double den1=0;
double den2=0;
double[] coeff=new double[6];
int count=0;
if (doThat){
sumA=0;
sumB=0;
sumAB=0;
sumsqrA=0;
Aarraymean=0;
Barraymean=0;
for (m=0; m=TA && Barray[m]>=TB){
sumA+=Aarray[m];
sumB+=Barray[m];
sumAB+=Aarray[m]*Barray[m];
sumsqrA+=Math.pow(Aarray[m],2);
count++;
}
}
Aarraymean=sumA/count;
Barraymean=sumB/count;
}
for (m=0; m=TA && Barray[m]>=TB){
num+=(Aarray[m]-Aarraymean)*(Barray[m]-Barraymean);
den1+=Math.pow((Aarray[m]-Aarraymean), 2);
den2+=Math.pow((Barray[m]-Barraymean), 2);
}
}
//0:a, 1:b, 2:corr coeff, 3: num, 4: den1, 5: den2
coeff[0]=(count*sumAB-sumA*sumB)/(count*sumsqrA-Math.pow(sumA,2));
coeff[1]=(sumsqrA*sumB-sumA*sumAB)/(count*sumsqrA-Math.pow(sumA,2));
coeff[2]=num/(Math.sqrt(den1*den2));
coeff[3]=num;
coeff[4]=den1;
coeff[5]=den2;
return coeff;
}
public double[] linregCostes(double[] Aarray, double[] Barray, double TA, double TB){
double num=0;
double den1=0;
double den2=0;
double[] coeff=new double[3];
int count=0;
sumA=0;
sumB=0;
sumAB=0;
sumsqrA=0;
Aarraymean=0;
Barraymean=0;
for (m=0; m