001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.image.processing.threshold;
031
032import org.openimaj.citation.annotation.Reference;
033import org.openimaj.citation.annotation.ReferenceType;
034import org.openimaj.image.FImage;
035import org.openimaj.image.processor.ImageProcessor;
036import org.openimaj.util.array.ArrayUtils;
037import org.openimaj.util.pair.FloatFloatPair;
038
039/**
040 * Otsu's adaptive thresholding algorithm.
041 * 
042 * @see <a
043 *      href="http://en.wikipedia.org/wiki/Otsu's_method">http://en.wikipedia.org/wiki/Otsu&apos;s_method</a>
044 * 
045 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
046 */
047@Reference(
048                type = ReferenceType.Article,
049                author = { "Nobuyuki Otsu" },
050                title = "A Threshold Selection Method from Gray-Level Histograms",
051                year = "1979",
052                journal = "Systems, Man and Cybernetics, IEEE Transactions on",
053                pages = { "62", "66" },
054                number = "1",
055                volume = "9",
056                customData = {
057                                "keywords",
058                                "Displays;Gaussian distribution;Histograms;Least squares approximation;Marine vehicles;Q measurement;Radar tracking;Sea measurements;Surveillance;Target tracking",
059                                "doi", "10.1109/TSMC.1979.4310076",
060                                "ISSN", "0018-9472"
061                })
062public class OtsuThreshold implements ImageProcessor<FImage> {
063        private static final int DEFAULT_NUM_BINS = 256;
064        int numBins = DEFAULT_NUM_BINS;
065
066        /**
067         * Default constructor
068         */
069        public OtsuThreshold() {
070
071        }
072
073        /**
074         * Construct with the given number of histogram bins
075         * 
076         * @param numBins
077         *            the number of histogram bins
078         */
079        public OtsuThreshold(int numBins) {
080                this.numBins = numBins;
081        }
082
083        protected static int[] makeHistogram(FImage fimg, int numBins) {
084                final int[] histData = new int[numBins];
085
086                // Calculate histogram
087                for (int r = 0; r < fimg.height; r++) {
088                        for (int c = 0; c < fimg.width; c++) {
089                                final int h = (int) (fimg.pixels[r][c] * (numBins - 1));
090                                histData[h]++;
091                        }
092                }
093
094                return histData;
095        }
096
097        protected static int[] makeHistogram(float[] data, int numBins, float min, float max) {
098                final int[] histData = new int[numBins];
099
100                // Calculate histogram
101                for (int c = 0; c < data.length; c++) {
102                        final float d = (data[c] - min) / (max - min);
103                        final int h = (int) (d * (numBins - 1));
104                        histData[h]++;
105                }
106
107                return histData;
108        }
109
110        /**
111         * Estimate the threshold for the given image.
112         * 
113         * @param img
114         *            the image
115         * @param numBins
116         *            the number of histogram bins
117         * @return the estimated threshold
118         */
119        public static float calculateThreshold(FImage img, int numBins) {
120                final int[] histData = makeHistogram(img, numBins);
121
122                // Total number of pixels
123                final int total = img.getWidth() * img.getHeight();
124
125                return computeThresholdFromHistogram(histData, total);
126        }
127
128        /**
129         * Estimate the threshold for the given data.
130         * <p>
131         * Internally, the data will be min-max normalised before the histogram is
132         * built, and the specified number of bins will cover the entire
133         * <code>max-min</code> range. The returned threshold will have
134         * <code>min</code> added to it to return it to the original range.
135         * 
136         * @param data
137         *            the data
138         * @param numBins
139         *            the number of histogram bins
140         * @return the estimated threshold
141         */
142        public static float calculateThreshold(float[] data, int numBins) {
143                final float min = ArrayUtils.minValue(data);
144                final float max = ArrayUtils.maxValue(data);
145                final int[] histData = makeHistogram(data, numBins, min, max);
146
147                return computeThresholdFromHistogram(histData, data.length) + min;
148        }
149
150        /**
151         * Estimate the threshold and inter-class variance for the given data.
152         * <p>
153         * Internally, the data will be min-max normalised before the histogram is
154         * built, and the specified number of bins will cover the entire
155         * <code>max-min</code> range. The returned threshold will have
156         * <code>min</code> added to it to return it to the original range.
157         * 
158         * @param data
159         *            the data
160         * @param numBins
161         *            the number of histogram bins
162         * @return the estimated threshold and variance
163         */
164        public static FloatFloatPair calculateThresholdAndVariance(float[] data, int numBins) {
165                final float min = ArrayUtils.minValue(data);
166                final float max = ArrayUtils.maxValue(data);
167                final int[] histData = makeHistogram(data, numBins, min, max);
168
169                final FloatFloatPair result = computeThresholdAndVarianceFromHistogram(histData, data.length);
170                result.first += min;
171                return result;
172        }
173
174        /**
175         * Estimate the threshold for the given histogram.
176         * 
177         * @param histData
178         *            the histogram
179         * @param total
180         *            the total number of items in the histogram
181         * @return the estimated threshold
182         */
183        public static float computeThresholdFromHistogram(int[] histData, int total) {
184                return computeThresholdAndVarianceFromHistogram(histData, total).first;
185        }
186
187        /**
188         * Estimate the threshold and inter-class variance for the given histogram.
189         * 
190         * @param histData
191         *            the histogram
192         * @param total
193         *            the total number of items in the histogram
194         * @return the estimated threshold and variance
195         */
196        public static FloatFloatPair computeThresholdAndVarianceFromHistogram(int[] histData, int total) {
197                final int numBins = histData.length;
198                float sum = 0;
199                for (int t = 0; t < numBins; t++)
200                        sum += t * histData[t];
201
202                float sumB = 0;
203                int wB = 0;
204                int wF = 0;
205
206                float varMax = 0;
207                float threshold = 0;
208
209                for (int t = 0; t < numBins; t++) {
210                        wB += histData[t]; // Weight Background
211                        if (wB == 0)
212                                continue;
213
214                        wF = total - wB; // Weight Foreground
215                        if (wF == 0)
216                                break;
217
218                        sumB += (t * histData[t]);
219
220                        final float mB = sumB / wB; // Mean Background
221                        final float mF = (sum - sumB) / wF; // Mean Foreground
222
223                        // Calculate Between Class Variance
224                        final float varBetween = (float) wB * (float) wF * (mB - mF) * (mB - mF);
225
226                        // Check if new maximum found
227                        if (varBetween > varMax) {
228                                varMax = varBetween;
229                                threshold = t;
230                        }
231                }
232
233                return new FloatFloatPair(threshold / (numBins - 1), varMax / total / total);
234        }
235
236        @Override
237        public void processImage(FImage image) {
238                final float threshold = calculateThreshold(image, numBins);
239
240                image.threshold(threshold);
241        }
242}