001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 * image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause).
006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009
010package imagingbook.common.histogram;
011
012import static imagingbook.common.math.Arithmetic.sqr;
013
014/**
015 * <p>
016 * This class defines static methods related to histograms. See Ch. 2 of [1] for additional details.
017 * </p>
018 * <p>
019 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
020 * (2022).
021 * </p>
022 *
023 * @author WB
024 * @version 2022/08/24
025 */
026public abstract class HistogramUtils {
027        
028        private HistogramUtils() { }
029
030        /**
031         * Calculates and returns the total population (sum of all bin counts) of a histogram.
032         *
033         * @param h a histogram
034         * @return the histogram's total count
035         */
036        public static int count(int[] h) {
037                return count(h, 0, h.length - 1);
038        }
039
040        /**
041         * Calculates and returns the population (sum of bin counts) of a histogram over the specified range of indexes. The
042         * range is automatically clipped.
043         *
044         * @param h a histogram
045         * @param lo the lower index (inclusive)
046         * @param hi the upper index (inclusive)
047         * @return the population count
048         */
049        public static int count(int[] h, int lo, int hi) {
050                if (lo < 0) lo = 0;
051                if (hi >= h.length) hi = h.length-1;
052                int cnt = 0;
053                for (int i = lo; i <= hi; i++) {
054                        cnt += h[i];
055                }
056                return cnt;
057        }
058        
059        // -----------------------------------------------------------
060        
061        /**
062         * Returns the maximum bin value (count) of the given histogram.
063         * 
064         * @param h a histogram
065         * @return the maximum bin value
066         */
067        public static int max(int[] h) {
068                int hmax = -1;
069                for (int i = 0; i < h.length; i++) {
070                        if (h[i] > hmax)
071                                hmax = h[i];
072                }
073                return hmax;
074        }
075
076        /**
077         * Returns the maximum bin value (count) of the given frequency distribution (histogram).
078         *
079         * @param h a histogram
080         * @return the maximum bin value
081         */
082        public static double max(double[] h) {
083                double hmax = Double.NEGATIVE_INFINITY;
084                for (int i = 0; i < h.length; i++) {
085                        if (h[i] > hmax)
086                                hmax = h[i];
087                }
088                return hmax;
089        }
090        
091        // -----------------------------------------------------------
092
093        /**
094         * Calculates and returns the cumulative histogram from a given histogram.
095         *
096         * @param h a histogram
097         * @return the cumulative histogram
098         */
099        public static int[] cumulate(int[] h) {
100                final int K = h.length;
101                int[] C = new int[K];
102                C[0] = h[0];
103                for (int i = 1; i < K; i++) {
104                C[i] = C[i-1] + h[i];
105            }
106            return C;
107        }
108        
109        // -----------------------------------------------------------
110
111        /**
112         * Calculates and returns the <em>probability distribution function</em> (pdf) for the given histogram. The
113         * resulting {@code double} array has the same length as the original histogram. Its values sum to 1.
114         *
115         * @param h a histogram
116         * @return the probability distribution function
117         */
118        public static double[] pdf(int[] h) {
119                final int K = h.length;
120                final int n = count(h);                 // sum all histogram values     
121                double[] p = new double[K];
122                for (int i = 0; i < h.length; i++) {
123                        p[i] =  (double) h[i] / n;
124                }
125                return p;
126        }
127
128        /**
129         * Calculates and returns the <em>cumulative distribution function</em> (cdf) for the given histogram. The resulting
130         * {@code double} array has the same length as the original histogram. Its maximum value is 1.
131         *
132         * @param h a histogram
133         * @return the cumulative distribution function
134         */
135        public static double[] cdf(int[] h) {
136                // returns the cumul. probability distribution function (cdf) for histogram h
137                final int K = h.length;
138                final int n = count(h);         // sum all histogram values             
139                double[] P = new double[K];
140                long c = h[0];
141                P[0] = (double) c / n;
142                for (int i = 1; i < K; i++) {
143                c = c + h[i];
144                P[i] = (double) c / n;
145            }
146            return P;
147        }
148        
149        // -----------------------------------------------------------
150
151        /**
152         * Calculates and returns the intensity mean (average) of the distribution represented by the given histogram.
153         *
154         * @param h a histogram
155         * @return the mean intensity
156         */
157        public static double mean(int[] h) {
158                return mean(h, 0, h.length - 1);
159        }
160
161        /**
162         * Calculates and returns the intensity mean (average) of the distribution represented by the given histogram,
163         * limited to the specified intensity range. The range is automatically clipped.
164         *
165         * @param h a histogram
166         * @param lo the lower index (inclusive)
167         * @param hi the upper index (inclusive)
168         * @return the mean intensity
169         */
170        public static double mean(int[] h, int lo, int hi) {
171                if (lo < 0)
172                        lo = 0;
173                if (hi >= h.length)
174                        hi = h.length - 1;
175                long cnt = 0;
176                long sum = 0;
177                for (int i = lo; i <= hi; i++) {
178                        cnt = cnt + h[i];
179                        sum = sum + i * h[i];
180                }
181                if (cnt > 0)
182                        return ((double) sum) / cnt;
183                else
184                        return 0;
185        }
186        
187        // -----------------------------------------------------------
188
189        /**
190         * Calculates and returns the intensity variance (&sigma;<sup>2</sup>) of the distribution represented by the given
191         * histogram.
192         *
193         * @param h a histogram
194         * @return the intensity variance
195         */
196        public static double variance(int[] h) {
197                return variance(h, 0, h.length-1);
198        }
199
200        /**
201         * Calculates and returns the intensity variance (&sigma;<sup>2</sup>) of the distribution represented by the given
202         * histogram, limited to the specified intensity range (fast version). The range is automatically clipped.
203         *
204         * @param h a histogram
205         * @param lo the lower index (inclusive)
206         * @param hi the upper index (inclusive)
207         * @return the intensity variance
208         */
209        public static double variance(int[] h, int lo, int hi) {
210                if (lo < 0)
211                        lo = 0;
212                if (hi >= h.length)
213                        hi = h.length - 1;
214                long A = 0;
215                long B = 0;
216                int N = 0;
217                for (int i = lo; i <= hi; i++) {
218                        int ni = h[i];
219                        A = A + (long) i * ni;
220                        B = B + (long) i * i * ni;
221                        N = N + ni;
222                }
223                
224                if (N == 0) {
225                        throw new IllegalArgumentException("empty histogram or range");
226                }
227                
228                return (B - (double) (A * A) / N) / N;
229        }
230        
231        // This is a naive (slow) version, for testing only:
232        public static double varianceSlow(int[] h, int lo, int hi) {
233                if (lo < 0)
234                        lo = 0;
235                if (hi >= h.length)
236                        hi = h.length - 1;
237                final double mu = mean(h, lo, hi);
238                int N = 0;
239                double sum = 0;
240                for (int i = lo; i <= hi; i++) {
241                        N = N + h[i];
242                        sum = sum + sqr(i - mu) * h[i];
243                }
244                
245                if (N == 0) {
246                        throw new IllegalArgumentException("empty histogram or range");
247                }
248
249                return sum / N;
250        }
251        
252        // -----------------------------------------------------------
253
254        /**
255         * Calculates and returns the intensity <em>median</em> of the distribution represented by the given histogram.
256         *
257         * @param h a histogram
258         * @return the intensity median
259         */
260        public static int median(int[] h) {
261                final int K = h.length;
262                final int N = count(h);
263                final int m = N / 2;
264                int i = 0;
265                int sum = h[0];
266                while (sum <= m && i < K) {
267                        i++;
268                        sum += h[i];
269                }
270                return i;
271        }
272        
273        // -----------------------------------------------------------
274
275        /**
276         * Returns a normalized frequency distribution for the given histogram whose maximum entry is 1 ({@code int}
277         * version). Mainly intended for displaying histograms.
278         *
279         * @param h a histogram
280         * @return the max-normalized frequency distribution
281         */
282        public static double[] normalizeMax(int[] h) {
283                // find the max histogram entry
284                final int max = max(h);
285                if (max == 0) {
286                        throw new IllegalArgumentException("empty histogram");
287                }
288                double[] hn = new double[h.length];
289                double s = 1.0 / max;
290                for (int i = 0; i < h.length; i++) {
291                        hn[i] = s * h[i];
292                }
293                return hn;
294        }
295
296        /**
297         * Returns a normalized frequency distribution for the given histogram whose maximum entry is 1 ({@code double}
298         * version). Mainly intended for displaying histograms.
299         *
300         * @param h a histogram
301         * @return the max-normalized frequency distribution
302         */
303        public static double[] normalizeMax(double[] h) {
304                // find max histogram entry
305                final double max = max(h);
306                if (max == 0) {
307                        throw new IllegalArgumentException("empty histogram");
308                }
309                double[] hn = new double[h.length];
310                double s = 1.0 / max;
311                for (int i = 0; i < h.length; i++) {
312                        hn[i] = s * h[i];
313                }
314                return hn;
315        }
316        
317        // -----------------------------------------------------------------
318
319        /**
320         * Histogram matching. Given are two histograms: the histogram hA of the target image IA and a reference histogram
321         * hR, both of size K. The result is a discrete mapping f which, when applied to the target image, produces a new
322         * image with a distribution function similar to the reference histogram.
323         *
324         * @param hA histogram of the target image
325         * @param hR reference histogram (the same size as hA)
326         * @return a discrete mapping f to be applied to the values of a target image
327         */
328        public static int[] matchHistograms (int[] hA, int[] hR) {
329                int K = hA.length;
330                double[] PA = HistogramUtils.cdf(hA); // get CDF of histogram hA
331                double[] PR = HistogramUtils.cdf(hR); // get CDF of histogram hR
332                int[] f = new int[K]; // pixel mapping function f()
333
334                // compute pixel mapping function f():
335                for (int a = 0; a < K; a++) {
336                        int j = K - 1;
337                        do {
338                                f[a] = j;
339                                j--;
340                        } while (j >= 0 && PA[a] <= PR[j]);
341                }
342                return f;
343        }
344
345        /**
346         * Histogram matching to a reference cumulative distribution function that is piecewise linear.
347         *
348         * @param hA histogram of the target image
349         * @param PR a piecewise linear reference cumulative distribution function ({@link PiecewiseLinearCdf})
350         * @return a discrete mapping f to be applied to the values of a target image
351         * @see PiecewiseLinearCdf
352         * @see #matchHistograms(int[], int[])
353         */
354        public static int[] matchHistograms(int[] hA, PiecewiseLinearCdf PR) {
355                int K = hA.length;
356                double[] PA = HistogramUtils.cdf(hA); // get p.d.f. of histogram Ha
357                int[] f = new int[K];                                   // pixel mapping function f()
358
359                // compute pixel mapping function f():
360                for (int a = 0; a < K; a++) {
361                        double b = PA[a];
362                        f[a] = (int) Math.round(PR.getInverseCdf(b));
363                }
364                return f;
365        }
366        
367}