Determining changes in soil core porosity with depth

Hello, I’m am new to ImageJ and this is my first forum post.

As part of a project I am trying to determine the porosity of a soil core and how it varies with depth. I have micro X-Ray CT images of soil cores. A re-slice of the raw image from one core is below.

Reslice of Core7 - for forum post.tif (892.3 KB)

I have applied the Otsu thresholding plug-in which provides a good differentiation between the air above the soil (black), the soil (white) and the soil pores (black) as shown in the image below.

Otsu threshold of Reslice of Core7 for forum post.tif (892.2 KB)

What I am trying to achieve is a calculation of the average porosity with depth below the soil surface. The reason for this is that I am studying the impact of soil crusts (which have lower porosity) upon infiltration of water into the soil profile.

As the soil surface is uneven I think that I need to flatten the soil surface in order to be able to use the analyse-measure command to calculate the porosity. Another method could be to use the colour difference between the air above the soil and the soil surface to create a line that follows the surface and then create regions of interest based on this line at regular depths below the soil surface to provide a region of interest to enable calculation of porosity.

I have not been able to work out how to achieve either of these methods and would welcome any feedback on how to do so.


Cameron from Australia

I would try to create a different segmentation with the aim to segment the soil from the air above. For example by applying a larger Gaussian Blur of lets say 10 px you blur the original image. Then you threshold that. You can clean up the smaller holes using Process > Binary > Fill holes…

You could also filter out the largest air region by first inverting the image. Then applying the Analyze Particles with a size filter with the Show: Masks option. Which then gives back a cleaned up Mask.

From this you can calculate the euclidean distance map using Process > Binary > Distance Map…

Each foreground pixel in the binary image is replaced with a gray value equal to that pixel’s distance from the nearest background pixel (for background pixels the EDM is 0)

Make sure the foreground and background color is specified correctly:
Edit > Options > Colors
Process > Binary > Options

Plus make sure you have EDM Output set to 16-bit under Process > Binary > Options.

This gives you the distance from the soil surface. Using thresholding with fixed upper and lower values you can define new ROIs with regular distances from that.

Then you can perform measurements within that ROI.



Brilliant, thanks Christopher. I will give this a go!



Hello Christopher,

I have tried following the second approach you detailed but have not come up with the final result that you have.

Firstly I inverted the image (image below).
[Reslice of Core7 - inverted.tif|attachment]
(upload://pp5tlkrSSxjQieezjHgwDBvzf4e.tif) (583.8 KB)

I then attempted to use the Analyse Particles (including the Show: Masks option). It said that a thresholded image or 8 bit binary image was required. I applied the Otsu filter (image below).

Reslice of Core7 - inverted with otsu filter applied.tif (583.7 KB)

To this image I then applied Analyse Particles resulting with the cleaned up mask (image below).

Mask of Reslice of Core7 - inverted with otsu filter applied.tif (585.3 KB)

I checked the Edit>Options>Colours settings which had foreground (black), background(blue) and selection (yellow). I was not sure what these should be set to so I left them as they were.

Under Process>Binary>Options Black background was selected, EDM output set at 16-bit and Do set to nothing.

I then created the distance map. The distance map when applied to the mask does not look like the distance map you created (attached below).
EDM of Mask Core7.tif (1.1 MB)

So I applied the EDM to the inverted and binarised original image. It looked more like the EDM image you attached but not quite.
EDM of Reslice Core7.tif (1.1 MB)

I appear to have missed a step/have a wrong setting here of which I am unsure of. As a result I cannot get to the final step of applying the thresholding.

I also don’t quite understand how to apply the final step of “Using thresholding with fixed upper and lower values you can define new ROIs with regular distances from that.” It would be greatly appreciated if you could expand upon that for me.

Any additional assistance would be greatly appreciated.

Thank you


Hi Cameron,

It really looks like you have somewhere a black background option that is (or not) ticked. It’s strange !

The distance map computes how far a pixel is from the closest object.
When you look at the image, you can see it makes a ramp (growing intensities) from the edges of the soil.
So if you use a threshold, e.g. 1 to 1, you will select all pixels that are 1 pixel away from the edge. with 1 to 8, you will get every pixels that are less than 9 pixels from the edge. with 13 to 21, you will get pixels that are at least 13 pixels away, but max 21 pixels away. etc.
So you can play with the lower and upper values to get the ROI you desire.


Hi Cameron,

the first step that you are missing it the large Gaussian Blur. Process > Filters > Gaussian Blur…
The aim is to really get rid of the holes in the soil sample to create a mask of the entire soil region. Also choose a threshold that eliminates the holes as well.

If that segmentation is not accurate enough for you, than you can always erode a bit the resulting mask. But be sure to set the binary options to pad edges when eroding, otherwise you get black edges and that messes up the EDM as well.

The settings for how the particle analyzer and the binary processing are interpreting a binary image are a bit messed up. It is confusing but that is not your fault. If you make sure that your foreground is always white and has 255 as pixel value and the background is black and has 0 as pixel value than the processing should work as indicated. Just try it out a bit.

For the last part I mean that you do not use a threshold method. You set fixed pixel values for the upper and lower bound of the threshold. This allows you to extract the part of the EDM that you require. You can then do this in intervals as you wish.

I wrote small macro to exemplify the workflow steps:

// set the options to:
// foreground = white
// background = black
run("Options...", "iterations=1 count=1 black edm=16-bit do=Nothing");
run("Colors...", "foreground=white background=black selection=yellow");

// select image and duplicate
selectWindow("Reslice of Core7 - for forum post.tif");
run("Duplicate...", " ");
selectWindow("Reslice of Core7 - for forum post-1.tif");

// needs lots of blur to get rid of holes
run("Gaussian Blur...", "sigma=10");

// threshold
// aim to get as a complete hole free soil as possible
setAutoThreshold("Default dark");
setOption("BlackBackground", true);
run("Convert to Mask");

// fill any remaining gaps
run("Fill Holes");

// Clean up any larger objects at the edge of the image
// invert image to get holes at the edge as objects

// size restricted analyze particles with masks as output
run("Analyze Particles...", "size=100-Infinity show=Masks");

// reset the LUT and binary to:
// foreground = white
// background = black
run("Invert LUT");

// create EDM
run("Distance Map");

// duplicate EDM
run("Duplicate...", " ");
selectWindow("EDM of Mask-1");

// set a manual threshold
// fixed lower and upper threshold value
// i.e. size in pixels of submask
setThreshold(1, 100);
setOption("BlackBackground", true);
run("Convert to Mask");
1 Like

Thanks Nico, that helps me understand how to manipulate the ROI.



Thanks so much Christopher.

I have used your code and it is achieving what I am after (and helping me to learn more about ImageJ).

Thanks again!


1 Like