Skip to content

Instantly share code, notes, and snippets.

@ekatrukha
Last active October 31, 2025 12:10
Show Gist options
  • Select an option

  • Save ekatrukha/6b0679cc6822daae32f22219fbfa54b6 to your computer and use it in GitHub Desktop.

Select an option

Save ekatrukha/6b0679cc6822daae32f22219fbfa54b6 to your computer and use it in GitHub Desktop.
Image Cytofluorogram Logarithmic Colocalization
nDefaultBinNumber = 256;
if (nImages<1)
{
exit("This macro needs an input image. Open something");
}
origID = getImageID();
Stack.getDimensions(width, height, channels, slices, frames);
imTitle = getTitle();
if(channels == 1)
{
exit("macro needs an image with at lease two channels");
}
sCh1Names = newArray(channels);
sCh2Names = newArray(channels);
for (i = 0; i < channels; i++)
{
sCh1Names[i] = toString(i+1);
sCh2Names[i] = toString(i+1);
}
bits = bitDepth();
if(bits == 24)
{
exit("RGB images are not supported, please convert to composite 8bit.");
}
nW = getWidth();
nH = getHeight();
sAxisChoice = newArray(2);
sAxisChoice[0] = "Y axis up";
sAxisChoice[1] = "Y axis down";
sScaleChoice = newArray(2);
sScaleChoice[0] = "Log10";
sScaleChoice[1] = "Linear";
Dialog.create("Params cytofluorogram");
Dialog.addMessage("Build a cytofluorogram, where ");
Dialog.addChoice("For X-axis use channel", sCh1Names);
Dialog.addChoice("Fir Y-axis use channel", sCh2Names);
Dialog.addChoice("Y-axis direction", sAxisChoice);
Dialog.addNumber("Bin number", nDefaultBinNumber);
Dialog.addChoice("Axes scale", sScaleChoice);
Dialog.addMessage("Normalization percentiles:");
Dialog.addNumber("Intensity MIN", 0.01);
Dialog.addNumber("Intensity MAX", 0.99);
Dialog.show();
chX = parseInt( Dialog.getChoice() );
chY = parseInt( Dialog.getChoice() );
sYAxisDown = Dialog.getChoice();
bYAxisDown = false;
if(sYAxisDown == "Y axis down")
{
bYAxisDown = true;
}
nBinNumber = Dialog.getNumber();
if(chX == chY)
{
IJ.log("AXIS X CHANNEL IS THE SAME AS AXIS Y CHANNEL");
}
bLogScale = false;
if(Dialog.getChoice() == "Log10")
{
bLogScale = true;
}
//setBatchMode(true);
valsChX = newArray(nW*nH);
valsChY = newArray(nW*nH);
selectImage(origID);
Stack.setChannel(chX);
run("Select All");
//get mode if linear mode
if(!bLogScale)
{
nModeX = getValue("Mode");
}
//store all values
for (i = 0; i < nW; i++)
{
for (j = 0; j < nH; j++)
{
valsChX[i*nH+j] = getPixel(i, j);
}
}
Stack.setChannel(chY);
run("Select All");
//get mode if linear mode
if(!bLogScale)
{
nModeY = getValue("Mode");
}
//store all values
for (i = 0; i < nW; i++)
{
for (j = 0; j < nH; j++)
{
valsChY[i*nH+j] = getPixel(i, j);
}
}
//arrays to sort values to determine percentiles
valsChXsorted = newArray(nW*nH);
valsChYsorted = newArray(nW*nH);
for (i = 0; i < nW; i++)
{
for (j = 0; j < nH; j++)
{
valsChXsorted[i*nH+j] = valsChX[i*nH+j];
valsChYsorted[i*nH+j] = valsChY[i*nH+j];
}
}
Array.sort(valsChXsorted);
Array.sort(valsChYsorted);
nPercMin = Dialog.getNumber();
nPercMax = Dialog.getNumber();
nMinIntX = valsChXsorted[Math.round(nPercMin*nW*nH)];
nMinIntY = valsChYsorted[Math.round(nPercMin*nW*nH)];
nMaxIntX = valsChXsorted[Math.round(nPercMax*nW*nH)];
nMaxIntY = valsChYsorted[Math.round(nPercMax*nW*nH)];
selectImage(origID);
if(bLogScale)
{
nMinIntX = Math.log10(nMinIntX);
nMinIntY = Math.log10(nMinIntY);
nMaxIntX = Math.log10(nMaxIntX);
nMaxIntY = Math.log10(nMaxIntY);
}
nBinSizeX = (nMaxIntX - nMinIntX+0.5) / nBinNumber;
nBinSizeY = (nMaxIntY - nMinIntY+0.5) / nBinNumber;
//get MODEs estimates (in pixel coordinates) in case of Log scale
if(bLogScale)
{
selectImage(origID);
Stack.setChannel(chX);
run("Select All");
twoGaussX = newArray(2);
nModeX = getModeLogScaleValue(nBinNumber, nBinSizeX, nMinIntX, nMaxIntX, true, twoGaussX);
selectImage(origID);
Stack.setChannel(chY);
run("Select All");
twoGaussY = newArray(2);
nModeY = getModeLogScaleValue(nBinNumber, nBinSizeY, nMinIntY, nMaxIntY,true, twoGaussY);
}
//get MODEs estimates (in pixel coordinates) in case of Linear scale
else
{
nModeX = Math.floor((nModeX-nMinIntX)/nBinSizeX);
nModeY = Math.floor((nModeY-nMinIntY)/nBinSizeY);
}
selectImage(origID);
if(bLogScale)
{
print("Intensity X MIN: " + toString(Math.pow(10,nMinIntX))+ " (log10 "+toString(nMinIntX)+")");
print("Intensity X MAX: " + toString(Math.pow(10,nMaxIntX))+ " (log10 "+toString(nMaxIntX)+")");
print("Intensity Y MIN: " + toString(Math.pow(10,nMinIntY)) + " (log10 "+toString(nMinIntY)+")");
print("Intensity Y MAX: " + toString(Math.pow(10,nMaxIntY))+ " (log10 "+toString(nMaxIntY)+")");
}
else
{
print("Intensity X MIN: " + toString(nMinIntX));
print("Intensity X MAX: " + toString(nMaxIntX));
print("Intensity Y MIN: " + toString(nMinIntY));
print("Intensity Y MAX: " + toString(nMaxIntY));
}
print("Intensity X BIN SIZE: " + toString(nBinSizeX));
print("Intensity Y BIN SIZE: " + toString(nBinSizeY));
//make an output image
nWout = nBinNumber;
nHout = nBinNumber;
print("Building CytoFluorogram between channels "+toString(chX)+" (X-axis) and " +toString(chY)+" (Y-axis)");
//print("Intensity MIN: " + toString(nMinInt));
//print("Intensity MAX: " + toString(nMaxInt));
//print("Intensity BIN SIZE: " + toString(nBinSize));
newImage("cytofluo_chx"+toString(chX)+"_vs_chy"+toString(chY)+"_"+imTitle, "16-bit black", nWout, nHout, 1);
outID = getImageID();
nCount = 0;
nMaxPixN = (nW*nH);
selectImage(outID);
for (i = 0; i < nW; i++)
{
for (j = 0; j < nH; j++)
{
valx = valsChX[i*nH+j];
valy = valsChY[i*nH+j];
if(bLogScale)
{
xout = Math.floor((Math.log10(valx)-nMinIntX)/nBinSizeX);
yout = Math.floor((Math.log10(valy)-nMinIntY)/nBinSizeY);
}
else
{
xout = Math.floor((valx-nMinIntX)/nBinSizeX);
yout = Math.floor((valy-nMinIntY)/nBinSizeY);
}
valcurr = getPixel(xout, yout);
setPixel(xout, yout, valcurr+1);
nCount++;
}
showProgress(nCount, nMaxPixN);
}
//setBatchMode(false);
selectImage(outID);
if(!bYAxisDown)
{
run("Flip Vertically");
nModeY = nHout - nModeY - 1.0;
if(bLogScale)
{
twoGaussY[0] = nHout - twoGaussY[0] - 1.0;
twoGaussY[1] = nHout - twoGaussY[1] - 1.0;
}
}
run("Fire");
run("Enhance Contrast", "saturated=0.35");
//makePoint(nModeX, nModeY, "large yellow hybrid");
//roiManager("add");
nPercCircle = 0.15;
if(bLogScale)
{
makeLine(twoGaussX[0], 0, twoGaussX[0], nH, 2);
//roiManager("add");
run("Add Selection...");
makeLine(twoGaussX[1], 0, twoGaussX[1], nH, 2);
//roiManager("add");
run("Add Selection...");
makeLine(0, twoGaussY[0], nW, twoGaussY[0], 2);
//roiManager("add");
run("Add Selection...");
makeLine(0, twoGaussY[1], nW, twoGaussY[1], 2);
//roiManager("add");
run("Add Selection...");
}
else
{
makeRectangle(nModeX-nBinNumber*0.5*nPercCircle, nModeY-nBinNumber*0.5*nPercCircle, nBinNumber*nPercCircle, nBinNumber*nPercCircle);
run("Add Selection...");
}
print("Done");
/** returns nMode already in pixels **/
function getModeLogScaleValue(nBinNumber, nBinSize, nMinInt, nMaxInt, bDoFit, fitResults)
{
//let's build Log histogram
hist = newArray(nBinNumber);
histBins = newArray(nBinNumber);
nW = getWidth();
nH = getHeight();
for (i = 0; i < nW; i++)
{
for (j = 0; j < nH; j++)
{
val = Math.log10(getPixel(i, j));
if(val>nMinInt && val<nMaxInt)
{
ind = Math.floor((val-nMinInt)/nBinSize);
hist[ind]++;
}
}
}
//find a maximum
dMax = -10000000;
dMin = 10000000;
dInd = 0;
for (i = 0; i < hist.length; i++)
{
histBins[i] = nMinInt + i*nBinSize;
if(hist[i]>dMax)
{
dMax = hist[i];
dInd = i;
}
if(hist[i]<dMin)
{
dMin = hist[i];
}
}
if(bDoFit)
{
//prepare fit
TwoGaussEquation = "y = a*Math.exp(-(x-b)*(x-b)/(2*c*c))+d*Math.exp(-(x-e)*(x-e)/(2*f*f))";
nRange = nMaxInt - nMinInt;
initialGuesses = newArray((dMax-dMin)*0.75, nMinInt+0.33*nRange, nRange/12, (dMax-dMin)*0.75, nMinInt+0.66*nRange, nRange/12) ;
Fit.doFit(TwoGaussEquation, histBins, hist, initialGuesses);
//Fit.plot;
fitResults[0] = Math.floor((Fit.p(1)-nMinInt)/nBinSize);
fitResults[1] = Math.floor((Fit.p(4)-nMinInt)/nBinSize);
}
//Plot.create("LogHist", "bins", "counts", histBins, hist);
//print(dInd);
//print(histBins[dInd]);
return dInd;
}
@ekatrukha
Copy link
Author

ekatrukha commented Sep 23, 2025

An ImageJ macro that takes as input a multi-channel image
and builds (Log or Linear) cytofluorogram between two selected channels.
The output is usually "normalized" between user-specified percentiles.

If logarithmic output is selected, it performs "sum-of-two-gaussians" fitting
on the log histogram in an attempt to find background and signal.
It can be used to visualize and quantify colocalizations between different structures.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment