Extended Depth of Field

File

Microscopes have an extremely narrow depth of field, with a 100x objective it may be as shallow as 1 micron. Unfortunately this makes capturing a 2D image of a 3D object at high magnification tricky! This macro is a basic example of an image stack focuser, it takes a stack of images from a range of focal depths and builds a 2D image from it using only in focus regions of the images, i.e. it generates an extended depth of field.

Extended_Depth_of_Field.ijm

//Generates a focused image from a z stack by selecting in-focus regions and automatically
//stitching them together to generate a single focused image.
//
//Required input:
//A greyscale (8, 16 or 32bit) z stack (with no sample movement between slices).
//For correct heightmap generation the slices must be in focus order.
//This macro assumes all pixels have positive values in 32bit images.
//
//Output:
//A 32bit greyscale heightmap based on sample focus
//A 32bit greyscale focused image generated from the most in focus regions
//
//Parameters:
//Radius - the radius of region to consisder for focusing.
//Large values are good to avoid the influence of noise.
//Small values are good for detecting and preserving fine detail.
//
//Principle:
//In focus regions have sharper detail, therefore have a stronger edge detection result.
//The slice with the maximum edge detection result is most in focus.
//Build up the image from the original stack in patches according to the most in focus slice.
//Use gaussian blending to reduce the appearance of sharp edges.
//This is similar to the plugin by Michael Umorin:
//http://rsbweb.nih.gov/ij/plugins/stack-focuser.html
//
//Tips:
//This works well even with relatively few z slices, and will process much faster.

//User variables
radius=3;

//Get start image properties
w=getWidth();
h=getHeight();
d=nSlices();
source=getImageID();
origtitle=getTitle();
rename("tempnameforprocessing");
sourcetitle=getTitle();
setBatchMode(true);

//Generate edge-detected image for detecting focus
run("Duplicate...", "title=["+sourcetitle+"_Heightmap] duplicate range=1-"+d);
heightmap=getImageID();
heightmaptitle=getTitle();
run("Find Edges", "stack");
run("Maximum...", "radius="+radius+" stack");

//Alter edge detected image to desired structure
run("32-bit");
for (x=0; x<w; x++) {
    showStatus("Finding Z depths");
    showProgress(x/w);
    for (y=0; y<h; y++) {
        slice=0;
        max=0;
        for (z=0; z<d; z++) {
            setZCoordinate(z);
            v=getPixel(x,y);
            if (v>=max) {
                max=v;
                slice=z;
            }
        }
        for (z=0; z<d; z++) {
            setZCoordinate(z);
            if (z==slice) {
                setPixel(x,y,1);
            } else {
                setPixel(x,y,0);
            }
        }
    }
}
run("Gaussian Blur...", "sigma="+radius+" stack");

//Generation of the final image
setBatchMode(false);
//Multiply modified edge detect (the depth map) with the source image
run("Image Calculator...", "image1="+sourcetitle+" operation=Multiply image2="+heightmaptitle+" create 32-bit stack");multiplication=getImageID();
//Z project the multiplication result
run("Z Project...", "start=1 stop="+d+" projection=[Sum Slices]");

//Some tidying
setBatchMode(true);
rename(origtitle+"_focused");
selectImage(heightmap);
close();
selectImage(multiplication);
close();
selectImage(source);
rename(origtitle);
setBatchMode(false);

ImageJ Macros Gallery