# Efficient way to perform local Otsu thresholding on image

I have images of cell nuclei that have locally varying intensity distributions. Therefore, I need a robust local threshold algorithm. LocalAdaptiveBinarize works quite sufficient for many cases, but still gives a lot of segmented objects in the background. Binarize can be used with Otsu's cluster variance maximization algorithm. I try to implement this algorithm using a local neighborhood in Mathematica. Here is a straightforward approach I came up with using ImageFilter and FindThreshold:

localOtsu[img_Image, r_Integer: 7] :=
ImageFilter[Boole[#[[r + 1, r + 1]] > FindThreshold[Flatten[#]]] &,
img, r]


Applied to a test image gives the following result and timing:

img = Import["...\\Test.tif"];
{img,AbsoluteTiming[Binarize[img]], AbsoluteTiming[localOtsu[img]]}


{,{0.040000,},{4.876000,}}

As you can see, this approach is terribly slow and I do not understand why this is the case. I appreciate any suggestions for improvement. Is ImageFilter the best way to do this?

With local Otsu thresholding I obtain much better results in many areas of my images than with LocalAdaptiveBinarize, no matter what parameters I use. That is why I am really looking for a way to implement it in Mathematica.

• as a workaround, you could calculate thresholds for every nth pixel, and interpolate in between. – Niki Estner Apr 15 '15 at 15:52

I would use some edge-preserving filter (like CurvatureFlowFilter) to get a smooth image before feeding to LocalAdaptiveBinarize:

img       = Import["Test.tif"]
smoothimg = CurvatureFlowFilter[img, 2]
biimg     = LocalAdaptiveBinarize[smoothimg, 10, {1, 0, .05}]
FlipView[Image[#, ImageSize -> 400] & /@ {img, cleanimg, biimg}]


• Thanks Silvia for this interesting approach. The problem with the CurvatureFlowFilter is that I also have objects in the image that are not circular. The structure of those gets distorted by the filter. As I said, I know about LocalAdaptiveBinarize, but local Otsu gives in many cases much better results (I tried this on images with ImageJ). That is why I am really looking for a implementation of it in Mathematica. – g3kk0 Apr 27 '15 at 6:21
• @g3kk0 Thanks. would it be possible to see the image with non-circular objects? There is no universal filter/algorithm in image processing, it is common (and often easier) to develop a special method for a specific problem. – Silvia Apr 27 '15 at 13:12

Update: Current implementation about twice as fast as OP, though there are some caveats.

Since you seem to be worried more about the speed of your method, maybe this could help. In the days of version 9 (which for me is still today) I implemented a local adaptive method similar to LocalAdaptiveBinarize, based only on the local mean of pixel values. I will look for the paper it came from later.

In this case it was done by computing the integrated image (This can be done with ImageAccumulate). This can be used to compute local mean in any rectangular region using four calls per pixel. This makes the problem O(N), where N is total number of pixels.

The code below, where the integrated image is computed explicitly, with For loops.

adaptiveAverageThres[img_, s_, t_, variance_: 0.0] :=
Module[{Intimg, xmin, xmax, ymin, ymax, dat1, sum, count, result},

(* Building an integrated image O(N),
this could all be replaced with ImageAccumulate for a huge speed boost,
though you may need to change the ordering, my summing is top-left to bottom right,
and ImageAccumulate is bottom-left to top-right *)

dat1 = ImageData@ColorNegate@ColorConvert[img, "Grayscale"];

Intimg = ConstantArray[0.0, Dimensions[dat1]];
Intimg[[1, 1 ;;]] = Accumulate[dat1[[1, 1 ;;]]];
Intimg[[1 ;;, 1]] = Accumulate[dat1[[1 ;;, 1]]];

For[i = 2, i <= First@Dimensions[dat1], i++,
For [j = 2, j <= Last@Dimensions[dat1], j++,
Intimg[[i, j]] = dat1[[i, j]] + Intimg[[i, j - 1]] +
Intimg[[i - 1, j]]-Intimg[[i - 1, j - 1]];
];
];

(* Computing average using the integrated image for the local sum
and comparing the threshold % values, O(N) *)

result = MapIndexed[
(xmin = Max[1, First@#2 - Ceiling[s/2.]];
xmax = Min[First@Dimensions[dat1], First@#2 + Ceiling[s/2.]];
ymin = Max[1, Last@#2 - Ceiling[s/2.]];
ymax = Min[Last@Dimensions[dat1], Last@#2 + Ceiling[s/2.]];
count = (xmax - xmin)*(ymax - min); (* number of pixels in region *)
sum = Intimg[[xmax, ymax]] - Intimg[[xmax, ymin]] -
Intimg[[xmin, ymax]] + Intimg[[xmin, ymin]];
(* implementing threshold based on the percentage away from mean*)
If[((#1*count) <= (sum*(100. - t)/100.)), 0.0,1.0]) &, dat1, {2}];
Return[result];
];


You can see how you could also adapt this to include the variance etc by computing the integrated square of pixel value. The idea would be for you to keep track of a quantity that can used to compute Otsu's method. Otsu's uses the histogram of pixel values, so that is a possibility. Since counting within each bin can be integrated with a simple sum.

## Local Otsu's

How can we change the implementation above to do Otsu's method? First problem is I don't know exactly how Otsu's method is implemented in FindThreshold, it is definitely histogramming the pixel values, but the binning is most likely depending on the distribution of values, so it may vary region by region.

Second problem is that Otsu's method requires us to go through each bin in the histogrammed pixel values in the specific region. So even if we have a O(N) algorithm pixel wise, an Otsu's version would have to be O(bins*N), which for small images, could be the leading speed.

Choosing too few bins, may merge two separate pixel values ranges into one, making the threshold useless. Too many and we can slow down considerably. The most bins we could have would be the encoded range of values, which for the tiff above (encoded in 16-bit), that's almost $10^5$ bins!

If we fix the binning to say 100 values, it should give us a reasonable amount of possible threshold values, one could increase to 8-bit range of 256 bins.

We proceed by computing an "Integrated Image" of the Cumulative Pixel values Histogram. This should allow us to compute the cumulative pixel histogram at any rectangular region with only 4 calls.

To speed up the Otsu's algorithm, we will also compute "Integrated Image" of the Cumulative Pixel histogram times the bin position also.

 dat1 = ImageData@img;

Intcumhist = ConstantArray[0.0,{First@Dimensions[dat1],Last@Dimensions[dat1], 100}];

bins=Table[i, {i, 0., 1., 1.0/100.}]
binsmidvalues = MovingAverage[bins, 2];
binpartitions = Partition[bins, 2, 1];

cumupixelhistogram[pix_] := cumupixelhistogram[pix] =
Accumulate@Boole[Thread[# <= pix < #2] & @@@ partitions];
cumupixeltimesbinvaluehistogram[pix_] := cumupixeltimesbinvaluehistogram[pix]=
(Accumulate@(#*binsmidvalues)&)@Boole[Thread[# <= pix < #2] & @@@ partitions];

Intcumhist[[1, 1 ;;]] = Accumulate[(cumpixelhistogram[#] &) /@ (dat1[[1, 1 ;;]])];
Intcumhist[[1 ;;, 1]] = Accumulate[(cumpixelhistogram[#] &) /@ (dat1[[1 ;;, 1]])];

Intcumhist2[[1, 1 ;;]] = Accumulate[(cumupixeltimesbinvaluehistogram[#] &) /@ (dat1[[1, 1 ;;]])];
Intcumhist2[[1 ;;, 1]] = Accumulate[(cumupixeltimesbinvaluehistogram[#] &) /@ (dat1[[1 ;;, 1]])];

For[i = 2, i <= First@Dimensions[dat1], i++,
For[j = 2, j <= Last@Dimensions[dat1], j++,
Intcumhist[[i, j]] = cumpixelhistogram[dat1[[i, j]]] +
Intcumhist[[i - 1, j]] + Intcumhist[[i, j - 1]] - Intcumhist[[i - 1, j - 1]];
Intcumhist2[[i, j]] = cumupixeltimesbinvaluehistogram[dat1[[i, j]]] +
Intcumhist2[[i - 1, j]] + Intcumhist2[[i, j - 1]] - Intcumhist2[[i - 1, j - 1]];
] ;
];


I computed the summation in the first row and column with Accumulate. I separated the first row and column, because one) it is really fast to write in mathematica, two) I didn't want to introduce if statements into the for loops. This is not the bottleneck of the algorithm yet, so I didn't study the timing. Right now its about 10% of Timing.

Now we need go to each pixel and compute the cumulative histogram for a window of size s*s=2*r+1 $\times$ 2*r +1. Then compute the Otsu's threshold value associated with each pixel window. Note that for your image having local size that doesn't include an object completely misses the point of Otsu's and will pick-up "noise". Also this part can be parallelize as we map at each pixel, and will give you the main speed boost.

For the example image, if you run more the 2 parallel kernels, starting up those kernels will actually make the problem slower! This is because the image is small. Once the kernels have started you get that speed boost regardless so it may be useful for large images and/or applying to a large set of images.

 r=16
newimgthres = ParallelTable[

(* computing window corners, note no padding atm *)

xmin = Max[1, row - r];
xmax = Min[First@Dimensions[dat1], row + r];
ymin = Max[1, column - r];
ymax = Min[Last@Dimensions[dat1], column + r];

cumul = Intcumhist[[xmax, ymax]] - Intcumhist[[xmax, ymin]] -
Intcumhist[[xmin, ymax]] + Intcumhist[[xmin, ymin]];

cumul2 = Intcumhist2[[xmax, ymax]] - Intcumhist2[[xmax, ymin]] -
Intcumhist2[[xmin, ymax]] + Intcumhist2[[xmin, ymin]];

cumulprob =N[cumul/(cumul[[-1]])]; (* cummulative probability distribution *)
cumul2prob = N[cumul2/(cumul2[[-1]])];

(* Otsu's Algorithm, Main bottleneck atm *)
max = 0.0;
binpos = 1;

Catch[Scan[
(ql = cumprob[[#]];
qh = cumprob[[-1]] - ql;
If[qh == 0.0, Throw[0.0]]; (* This is true if we reach the max nonzero bin in the region,
break out scan *)
If[ql != 0.0,
res = ql*(cum2prob[[#]]*cumprob[[-1]]/ql - cum2prob[[-1]])^2/qh;
If[max < res, binpos = #; max = res;];
]; &), numberofbins]];

binsmidvalues[[binpos]],
{column, 1, Last@Dimensions[dat1], 1}, {row, 1, First@Dimensions[dat1], 1}] ];


From time to time I will come back to it. At the moment with region size r=16

AbsoluteTiming[yourresult=localOtsu[imgtest, 16]]


2.655135

while the timing for the above for 100 bins, and with 4 Parallel Kernels is (not including the time taken to Launch them).

1.219444

Just to give you an idea, building the integrated image of both cumulative objects takes about (AbsoluteTiming)

~0.196480 seconds

while applying Otsu at each pixel,

~1.018954 seconds

result is similar for the same values of r

  GraphicsGrid[{{yourresult,
Image@MapThread[Boole[#2 > #1] &, {newimgthres, ImageData@imgtest}, 2]}}]


The differences are mainly due to the Standard padding differences.

## Some Possible Improvements and Suggestions.

I didn't include padding but you could apply a padding with ArrayPad, compute the integrated images while computing the threshold only in the non-padded regions.

Note that for any local thresholding method: It is rare that there is a reason thresholds should vary too much pixel to pixel, so it may be useful to apply a gaussian filter (or some other smoothing filter) on the image that contains only the threshold values.

Another possibility is to do what @nikie suggested, only apply the threshold method in a few regions and then Interpolate the result.

A third possibility is to apply these at the GPU level with CUDA/OpenCL since the second part is trivially parallelized.