| /*====================================================================* |
| - Copyright (C) 2001 Leptonica. All rights reserved. |
| - This software is distributed in the hope that it will be |
| - useful, but with NO WARRANTY OF ANY KIND. |
| - No author or distributor accepts responsibility to anyone for the |
| - consequences of using this software, or for whether it serves any |
| - particular purpose or works at all, unless he or she says so in |
| - writing. Everyone is granted permission to copy, modify and |
| - redistribute this source code, for commercial or non-commercial |
| - purposes, with the following restrictions: (1) the origin of this |
| - source code must not be misrepresented; (2) modified versions must |
| - be plainly marked as such; and (3) this notice may not be removed |
| - or altered from any source or modified source distribution. |
| *====================================================================*/ |
| |
| /* |
| * numafunc2.c |
| * |
| * Transformations |
| * NUMA *numaTransform() |
| * NUMA *numaConvolve() |
| * NUMA *numaConvertToInt() |
| * |
| * Histogram generation and statistics |
| * NUMA *numaMakeHistogram() |
| * NUMA *numaMakeHistogramAuto() |
| * NUMA *numaMakeHistogramClipped() |
| * NUMA *numaRebinHistogram() |
| * NUMA *numaNormalizeHistogram() |
| * l_int32 numaGetStatsUsingHistogram() |
| * l_int32 numaGetHistogramStats() |
| * l_int32 numaGetHistogramStatsOnInterval() |
| * l_int32 numaMakeRankFromHistogram() |
| * l_int32 numaHistogramGetRankFromVal() |
| * l_int32 numaHistogramGetValFromRank() |
| * |
| * Splitting a distribution |
| * l_int32 numaSplitDistribution() |
| * |
| * Extrema finding |
| * NUMA *numaFindPeaks() |
| * NUMA *numaFindExtrema() |
| * |
| * Threshold crossings and frequency analysis |
| * l_int32 numaSelectCrossingThreshold() |
| * NUMA *numaCrossingsByThreshold() |
| * NUMA *numaCrossingsByPeaks() |
| * NUMA *numaEvalBestHaarParameters() |
| * l_int32 numaEvalHaarSum() |
| * |
| * Things to remember when using the Numa: |
| * |
| * (1) The numa is a struct, not an array. Always use accessors |
| * (see numabasic.c), never the fields directly. |
| * |
| * (2) The number array holds l_float32 values. It can also |
| * be used to store l_int32 values. See numabasic.c for |
| * details on using the accessors. |
| * |
| * (3) Occasionally, in the comments we denote the i-th element of a |
| * numa by na[i]. This is conceptual only -- the numa is not an array! |
| * |
| * Some general comments on histograms: |
| * |
| * (1) Histograms are the generic statistical representation of |
| * the data about some attribute. Typically they're not |
| * normalized -- they simply give the number of occurrences |
| * within each range of values of the attribute. This range |
| * of values is referred to as a 'bucket'. For example, |
| * the histogram could specify how many connected components |
| * are found for each value of their width; in that case, |
| * the bucket size is 1. |
| * |
| * (2) In leptonica, all buckets have the same size. Histograms |
| * are therefore specified by a numa of occurrences, along |
| * with two other numbers: the 'value' associated with the |
| * occupants of the first bucket and the size (i.e., 'width') |
| * of each bucket. These two numbers then allow us to calculate |
| * the value associated with the occupants of each bucket. |
| * These numbers are fields in the numa, initialized to |
| * a startx value of 0.0 and a binsize of 1.0. Accessors for |
| * these fields are functions numa*XParameters(). All histograms |
| * must have these two numbers properly set. |
| */ |
| |
| #include <stdio.h> |
| #include <string.h> |
| #include <stdlib.h> |
| #include <math.h> |
| #include "allheaders.h" |
| |
| /* bin sizes in numaMakeHistogram() */ |
| static const l_int32 BinSizeArray[] = {2, 5, 10, 20, 50, 100, 200, 500, 1000,\ |
| 2000, 5000, 10000, 20000, 50000, 100000, 200000,\ |
| 500000, 1000000, 2000000, 5000000, 10000000,\ |
| 200000000, 50000000, 100000000}; |
| static const l_int32 NBinSizes = 24; |
| |
| |
| #ifndef NO_CONSOLE_IO |
| #define DEBUG_HISTO 0 |
| #define DEBUG_CROSSINGS 0 |
| #define DEBUG_FREQUENCY 0 |
| #endif /* ~NO_CONSOLE_IO */ |
| |
| |
| /*----------------------------------------------------------------------* |
| * Transformations * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * numaTransform() |
| * |
| * Input: nas |
| * shift (add this to each number) |
| * scale (multiply each number by this) |
| * Return: na with all values shifted and scaled, or null on error |
| * |
| * Notes: |
| * (1) Each number is shifted before scaling. |
| * (2) The operation sequence is opposite to that for Box and Pta: |
| * scale first, then shift. |
| */ |
| NUMA * |
| numaTransform(NUMA *nas, |
| l_float32 shift, |
| l_float32 scale) |
| { |
| l_int32 i, n; |
| l_float32 val; |
| NUMA *nad; |
| |
| PROCNAME("numaTransform"); |
| |
| if (!nas) |
| return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); |
| n = numaGetCount(nas); |
| if ((nad = numaCreate(n)) == NULL) |
| return (NUMA *)ERROR_PTR("nad not made", procName, NULL); |
| for (i = 0; i < n; i++) { |
| numaGetFValue(nas, i, &val); |
| val = scale * val + shift; |
| numaAddNumber(nad, val); |
| } |
| return nad; |
| } |
| |
| |
| /*! |
| * numaConvolve() |
| * |
| * Input: na |
| * halfwidth (of rectangular filter, minus the center) |
| * Return: na (after low-pass filtering), or null on error |
| * |
| * Notes: |
| * (1) Full convolution takes place only from i = halfwidth to |
| * i = n - halfwidth - 1. We do the end parts using only |
| * the partial array available. We do not pad the ends with 0. |
| * (2) This implementation assumes specific fields in the Numa! |
| */ |
| NUMA * |
| numaConvolve(NUMA *na, |
| l_int32 halfwidth) |
| { |
| l_int32 i, n, rval; |
| l_float32 sum, norm; |
| l_float32 *array, *carray, *sumarray; |
| NUMA *nac; |
| |
| PROCNAME("numaConvolve"); |
| |
| if (!na) |
| return (NUMA *)ERROR_PTR("na not defined", procName, NULL); |
| n = numaGetCount(na); |
| if (2 * halfwidth + 1 > n) |
| L_WARNING("filter wider than input array!", procName); |
| array = na->array; |
| |
| if ((nac = numaCreate(n)) == NULL) |
| return (NUMA *)ERROR_PTR("nac not made", procName, NULL); |
| carray = nac->array; |
| nac->n = n; /* fill with zeroes */ |
| |
| /* Make sum array; note the indexing */ |
| if ((sumarray = (l_float32 *)CALLOC(n + 1, sizeof(l_float32))) == NULL) |
| return (NUMA *)ERROR_PTR("sumarray not made", procName, NULL); |
| sum = 0.0; |
| sumarray[0] = 0.0; |
| for (i = 0; i < n; i++) { |
| sum += array[i]; |
| sumarray[i + 1] = sum; |
| } |
| |
| /* Central part */ |
| norm = 1. / (2 * halfwidth + 1); |
| rval = n - halfwidth; |
| for (i = halfwidth; i < rval; ++i) |
| carray[i] = norm * |
| (sumarray[i + halfwidth + 1] - sumarray[i - halfwidth]); |
| |
| /* Left side */ |
| for (i = 0; i < halfwidth; i++) |
| carray[i] = sumarray[i + halfwidth + 1] / (halfwidth + i + 1); |
| |
| /* Right side */ |
| for (i = rval; i < n; i++) |
| carray[i] = (1. / (n - i + halfwidth)) * |
| (sumarray[n] - sumarray[i - halfwidth]); |
| |
| FREE(sumarray); |
| return nac; |
| } |
| |
| |
| /*! |
| * numaConvertToInt() |
| * |
| * Input: na |
| * Return: na with all values rounded to nearest integer, or |
| * null on error |
| */ |
| NUMA * |
| numaConvertToInt(NUMA *nas) |
| { |
| l_int32 i, n, ival; |
| NUMA *nad; |
| |
| PROCNAME("numaConvertToInt"); |
| |
| if (!nas) |
| return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); |
| |
| n = numaGetCount(nas); |
| if ((nad = numaCreate(n)) == NULL) |
| return (NUMA *)ERROR_PTR("nad not made", procName, NULL); |
| for (i = 0; i < n; i++) { |
| numaGetIValue(nas, i, &ival); |
| numaAddNumber(nad, ival); |
| } |
| return nad; |
| } |
| |
| |
| /*----------------------------------------------------------------------* |
| * Histogram generation and statistics * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * numaMakeHistogram() |
| * |
| * Input: na |
| * maxbins (max number of histogram bins) |
| * &binsize (<return> size of histogram bins) |
| * &binstart (<optional return> start val of minimum bin; |
| * input NULL to force start at 0) |
| * Return: na consisiting of histogram of integerized values, |
| * or null on error. |
| * |
| * Note: |
| * (1) This simple interface is designed for integer data. |
| * The bins are of integer width and start on integer boundaries, |
| * so the results on float data will not have high precision. |
| * (2) Specify the max number of input bins. Then @binsize, |
| * the size of bins necessary to accommodate the input data, |
| * is returned. It is one of the sequence: |
| * {1, 2, 5, 10, 20, 50, ...}. |
| * (3) If &binstart is given, all values are accommodated, |
| * and the min value of the starting bin is returned. |
| * Otherwise, all negative values are discarded and |
| * the histogram bins start at 0. |
| */ |
| NUMA * |
| numaMakeHistogram(NUMA *na, |
| l_int32 maxbins, |
| l_int32 *pbinsize, |
| l_int32 *pbinstart) |
| { |
| l_int32 i, n, ival, hval; |
| l_int32 iminval, imaxval, range, binsize, nbins, ibin; |
| l_float32 val, ratio; |
| NUMA *nai, *nahist; |
| |
| PROCNAME("numaMakeHistogram"); |
| |
| if (!na) |
| return (NUMA *)ERROR_PTR("na not defined", procName, NULL); |
| if (!pbinsize) |
| return (NUMA *)ERROR_PTR("&binsize not defined", procName, NULL); |
| |
| /* Determine input range */ |
| numaGetMin(na, &val, NULL); |
| iminval = (l_int32)(val + 0.5); |
| numaGetMax(na, &val, NULL); |
| imaxval = (l_int32)(val + 0.5); |
| if (pbinstart == NULL) { /* clip negative vals; start from 0 */ |
| iminval = 0; |
| if (imaxval < 0) |
| return (NUMA *)ERROR_PTR("all values < 0", procName, NULL); |
| } |
| |
| /* Determine binsize */ |
| range = imaxval - iminval + 1; |
| if (range > maxbins - 1) { |
| ratio = (l_float64)range / (l_float64)maxbins; |
| binsize = 0; |
| for (i = 0; i < NBinSizes; i++) { |
| if (ratio < BinSizeArray[i]) { |
| binsize = BinSizeArray[i]; |
| break; |
| } |
| } |
| if (binsize == 0) |
| return (NUMA *)ERROR_PTR("numbers too large", procName, NULL); |
| } |
| else |
| binsize = 1; |
| *pbinsize = binsize; |
| nbins = 1 + range / binsize; /* +1 seems to be sufficient */ |
| |
| /* Redetermine iminval */ |
| if (pbinstart && binsize > 1) { |
| if (iminval >= 0) |
| iminval = binsize * (iminval / binsize); |
| else |
| iminval = binsize * ((iminval - binsize + 1) / binsize); |
| } |
| if (pbinstart) |
| *pbinstart = iminval; |
| |
| #if DEBUG_HISTO |
| fprintf(stderr, " imaxval = %d, range = %d, nbins = %d\n", |
| imaxval, range, nbins); |
| #endif /* DEBUG_HISTO */ |
| |
| /* Use integerized data for input */ |
| if ((nai = numaConvertToInt(na)) == NULL) |
| return (NUMA *)ERROR_PTR("nai not made", procName, NULL); |
| n = numaGetCount(nai); |
| |
| /* Make histogram, converting value in input array |
| * into a bin number for this histogram array. */ |
| if ((nahist = numaCreate(nbins)) == NULL) |
| return (NUMA *)ERROR_PTR("nahist not made", procName, NULL); |
| numaSetCount(nahist, nbins); |
| numaSetXParameters(nahist, iminval, binsize); |
| for (i = 0; i < n; i++) { |
| numaGetIValue(nai, i, &ival); |
| ibin = (ival - iminval) / binsize; |
| if (ibin >= 0 && ibin < nbins) { |
| numaGetIValue(nahist, ibin, &hval); |
| numaSetValue(nahist, ibin, hval + 1.0); |
| } |
| } |
| |
| numaDestroy(&nai); |
| return nahist; |
| } |
| |
| |
| /*! |
| * numaMakeHistogramAuto() |
| * |
| * Input: na (numa of floats; these may be integers) |
| * maxbins (max number of histogram bins; >= 1) |
| * Return: na consisiting of histogram of quantized float values, |
| * or null on error. |
| * |
| * Notes: |
| * (1) This simple interface is designed for accurate binning |
| * of both integer and float data. |
| * (2) If the array data is integers, and the range of integers |
| * is smaller than @maxbins, they are binned as they fall, |
| * with binsize = 1. |
| * (3) If the range of data, (maxval - minval), is larger than |
| * @maxbins, or if the data is floats, they are binned into |
| * exactly @maxbins bins. |
| * (4) Unlike numaMakeHistogram(), these bins in general have |
| * non-integer location and width, even for integer data. |
| */ |
| NUMA * |
| numaMakeHistogramAuto(NUMA *na, |
| l_int32 maxbins) |
| { |
| l_int32 i, n, imin, imax, irange, ibin, ival, allints; |
| l_float32 minval, maxval, range, binsize, fval; |
| NUMA *nah; |
| |
| PROCNAME("numaMakeHistogramAuto"); |
| |
| if (!na) |
| return (NUMA *)ERROR_PTR("na not defined", procName, NULL); |
| maxbins = L_MAX(1, maxbins); |
| |
| /* Determine input range */ |
| numaGetMin(na, &minval, NULL); |
| numaGetMax(na, &maxval, NULL); |
| |
| /* Determine if values are all integers */ |
| n = numaGetCount(na); |
| numaHasOnlyIntegers(na, maxbins, &allints); |
| |
| /* Do simple integer binning if possible */ |
| if (allints && (maxval - minval < maxbins)) { |
| imin = (l_int32)minval; |
| imax = (l_int32)maxval; |
| irange = imax - imin + 1; |
| nah = numaCreate(irange); |
| numaSetCount(nah, irange); /* init */ |
| numaSetXParameters(nah, minval, 1.0); |
| for (i = 0; i < n; i++) { |
| numaGetIValue(na, i, &ival); |
| ibin = ival - imin; |
| numaGetIValue(nah, ibin, &ival); |
| numaSetValue(nah, ibin, ival + 1.0); |
| } |
| |
| return nah; |
| } |
| |
| /* Do float binning, even if the data is integers. */ |
| range = maxval - minval; |
| binsize = range / (l_float32)maxbins; |
| if (range == 0.0) { |
| nah = numaCreate(1); |
| numaSetXParameters(nah, minval, binsize); |
| numaAddNumber(nah, n); |
| return nah; |
| } |
| nah = numaCreate(maxbins); |
| numaSetCount(nah, maxbins); |
| numaSetXParameters(nah, minval, binsize); |
| for (i = 0; i < n; i++) { |
| numaGetFValue(na, i, &fval); |
| ibin = (l_int32)((fval - minval) / binsize); |
| numaGetIValue(nah, ibin, &ival); |
| numaSetValue(nah, ibin, ival + 1.0); |
| } |
| |
| return nah; |
| } |
| |
| |
| /*! |
| * numaMakeHistogramClipped() |
| * |
| * Input: na |
| * binsize (typically 1.0) |
| * maxsize (of histogram ordinate) |
| * Return: na (histogram of bins of size @binsize, starting with |
| * the na[0] (x = 0.0) and going up to a maximum of |
| * x = @maxsize, by increments of @binsize), or null on error |
| * |
| * Notes: |
| * (1) This simple function generates a histogram of values |
| * from na, discarding all values < 0.0 or greater than |
| * min(@maxsize, maxval), where maxval is the maximum value in na. |
| * The histogram data is put in bins of size delx = @binsize, |
| * starting at x = 0.0. We use as many bins as are |
| * needed to hold the data. |
| */ |
| NUMA * |
| numaMakeHistogramClipped(NUMA *na, |
| l_float32 binsize, |
| l_float32 maxsize) |
| { |
| l_int32 i, n, nbins, ival, ibin; |
| l_float32 val, maxval; |
| NUMA *nad; |
| |
| PROCNAME("numaMakeHistogramClipped"); |
| |
| if (!na) |
| return (NUMA *)ERROR_PTR("na not defined", procName, NULL); |
| if (binsize <= 0.0) |
| return (NUMA *)ERROR_PTR("binsize must be > 0.0", procName, NULL); |
| if (binsize > maxsize) |
| binsize = maxsize; /* just one bin */ |
| |
| numaGetMax(na, &maxval, NULL); |
| n = numaGetCount(na); |
| maxsize = L_MIN(maxsize, maxval); |
| nbins = (l_int32)(maxsize / binsize) + 1; |
| |
| /* fprintf(stderr, "maxsize = %7.3f, nbins = %d\n", maxsize, nbins); */ |
| |
| if ((nad = numaCreate(nbins)) == NULL) |
| return (NUMA *)ERROR_PTR("nad not made", procName, NULL); |
| numaSetXParameters(nad, 0.0, binsize); |
| numaSetCount(nad, nbins); /* interpret zeroes in bins as data */ |
| for (i = 0; i < n; i++) { |
| numaGetFValue(na, i, &val); |
| ibin = (l_int32)(val / binsize); |
| if (ibin >= 0 && ibin < nbins) { |
| numaGetIValue(nad, ibin, &ival); |
| numaSetValue(nad, ibin, ival + 1.0); |
| } |
| } |
| |
| return nad; |
| } |
| |
| |
| /*! |
| * numaRebinHistogram() |
| * |
| * Input: nas (input histogram) |
| * newsize (number of old bins contained in each new bin) |
| * Return: nad (more coarsely re-binned histogram), or null on error |
| */ |
| NUMA * |
| numaRebinHistogram(NUMA *nas, |
| l_int32 newsize) |
| { |
| l_int32 i, j, ns, nd, index, count, val; |
| l_float32 start, oldsize; |
| NUMA *nad; |
| |
| PROCNAME("numaRebinHistogram"); |
| |
| if (!nas) |
| return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); |
| if (newsize <= 1) |
| return (NUMA *)ERROR_PTR("newsize must be > 1", procName, NULL); |
| if ((ns = numaGetCount(nas)) == 0) |
| return (NUMA *)ERROR_PTR("no bins in nas", procName, NULL); |
| |
| nd = (ns + newsize - 1) / newsize; |
| if ((nad = numaCreate(nd)) == NULL) |
| return (NUMA *)ERROR_PTR("nad not made", procName, NULL); |
| numaGetXParameters(nad, &start, &oldsize); |
| numaSetXParameters(nad, start, oldsize * newsize); |
| |
| for (i = 0; i < nd; i++) { /* new bins */ |
| count = 0; |
| index = i * newsize; |
| for (j = 0; j < newsize; j++) { |
| if (index < ns) { |
| numaGetIValue(nas, index, &val); |
| count += val; |
| index++; |
| } |
| } |
| numaAddNumber(nad, count); |
| } |
| |
| return nad; |
| } |
| |
| |
| /*! |
| * numaNormalizeHistogram() |
| * |
| * Input: nas (input histogram) |
| * area (target sum of all numbers in dest histogram; |
| * e.g., use area = 1.0 if this represents a |
| * probability distribution) |
| * Return: nad (normalized histogram), or null on error |
| */ |
| NUMA * |
| numaNormalizeHistogram(NUMA *nas, |
| l_float32 area) |
| { |
| l_int32 i, ns; |
| l_float32 sum, factor, fval; |
| NUMA *nad; |
| |
| PROCNAME("numaNormalizeHistogram"); |
| |
| if (!nas) |
| return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); |
| if (area <= 0.0) |
| return (NUMA *)ERROR_PTR("area must be > 0.0", procName, NULL); |
| if ((ns = numaGetCount(nas)) == 0) |
| return (NUMA *)ERROR_PTR("no bins in nas", procName, NULL); |
| |
| numaGetSum(nas, &sum); |
| factor = area / sum; |
| |
| if ((nad = numaCreate(ns)) == NULL) |
| return (NUMA *)ERROR_PTR("nad not made", procName, NULL); |
| numaCopyXParameters(nad, nas); |
| |
| for (i = 0; i < ns; i++) { |
| numaGetFValue(nas, i, &fval); |
| fval *= factor; |
| numaAddNumber(nad, fval); |
| } |
| |
| return nad; |
| } |
| |
| |
| /*! |
| * numaGetStatsUsingHistogram() |
| * |
| * Input: na (an arbitrary set of numbers; not ordered and not |
| * a histogram) |
| * maxbins (the maximum number of bins to be allowed in |
| * the histogram; use 0 for consecutive integer bins) |
| * &min (<optional return> min value of set) |
| * &max (<optional return> max value of set) |
| * &mean (<optional return> mean value of set) |
| * &variance (<optional return> variance) |
| * &median (<optional return> median value of set) |
| * rank (in [0.0 ... 1.0]; median has a rank 0.5; ignored |
| * if &rval == NULL) |
| * &rval (<optional return> value in na corresponding to @rank) |
| * &histo (<optional return> Numa histogram; use NULL to prevent) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) This is a simple interface for gathering statistics |
| * from a numa, where a histogram is used 'under the covers' |
| * to avoid sorting if a rank value is requested. In that case, |
| * by using a histogram we are trading speed for accuracy, because |
| * the values in @na are quantized to the center of a set of bins. |
| * (2) If the median, other rank value, or histogram are not requested, |
| * the calculation is all performed on the input Numa. |
| * (3) The variance is the average of the square of the |
| * difference from the mean. The median is the value in na |
| * with rank 0.5. |
| * (4) There are two situations where this gives rank results with |
| * accuracy comparable to computing stastics directly on the input |
| * data, without binning into a histogram: |
| * (a) the data is integers and the range of data is less than |
| * @maxbins, and |
| * (b) the data is floats and the range is small compared to |
| * @maxbins, so that the binsize is much less than 1. |
| * (5) If a histogram is used and the numbers in the Numa extend |
| * over a large range, you can limit the required storage by |
| * specifying the maximum number of bins in the histogram. |
| * Use @maxbins == 0 to force the bin size to be 1. |
| * (6) This optionally returns the median and one arbitrary rank value. |
| * If you need several rank values, return the histogram and use |
| * numaHistogramGetValFromRank(nah, rank, &rval) |
| * multiple times. |
| */ |
| l_int32 |
| numaGetStatsUsingHistogram(NUMA *na, |
| l_int32 maxbins, |
| l_float32 *pmin, |
| l_float32 *pmax, |
| l_float32 *pmean, |
| l_float32 *pvariance, |
| l_float32 *pmedian, |
| l_float32 rank, |
| l_float32 *prval, |
| NUMA **phisto) |
| { |
| l_int32 i, n; |
| l_float32 minval, maxval, fval, mean, sum; |
| NUMA *nah; |
| |
| PROCNAME("numaGetStatsUsingHistogram"); |
| |
| if (pmin) *pmin = 0.0; |
| if (pmax) *pmax = 0.0; |
| if (pmean) *pmean = 0.0; |
| if (pmedian) *pmedian = 0.0; |
| if (pvariance) *pvariance = 0.0; |
| if (!na) |
| return ERROR_INT("na not defined", procName, 1); |
| if ((n = numaGetCount(na)) == 0) |
| return ERROR_INT("numa is empty", procName, 1); |
| |
| numaGetMin(na, &minval, NULL); |
| numaGetMax(na, &maxval, NULL); |
| if (pmin) *pmin = minval; |
| if (pmax) *pmax = maxval; |
| if (pmean || pvariance) { |
| sum = 0.0; |
| for (i = 0; i < n; i++) { |
| numaGetFValue(na, i, &fval); |
| sum += fval; |
| } |
| mean = sum / (l_float32)n; |
| if (pmean) *pmean = mean; |
| } |
| if (pvariance) { |
| sum = 0.0; |
| for (i = 0; i < n; i++) { |
| numaGetFValue(na, i, &fval); |
| sum += fval * fval; |
| } |
| *pvariance = sum / (l_float32)n - mean * mean; |
| } |
| |
| if (!pmedian && !prval && !phisto) |
| return 0; |
| |
| nah = numaMakeHistogramAuto(na, maxbins); |
| if (pmedian) |
| numaHistogramGetValFromRank(nah, 0.5, pmedian); |
| if (prval) |
| numaHistogramGetValFromRank(nah, rank, prval); |
| if (phisto) |
| *phisto = nah; |
| else |
| numaDestroy(&nah); |
| return 0; |
| } |
| |
| |
| /*! |
| * numaGetHistogramStats() |
| * |
| * Input: nahisto (histogram: y(x(i)), i = 0 ... nbins - 1) |
| * startx (x value of first bin: x(0)) |
| * deltax (x increment between bins; the bin size; x(1) - x(0)) |
| * &xmean (<optional return> mean value of histogram) |
| * &xmedian (<optional return> median value of histogram) |
| * &xmode (<optional return> mode value of histogram: |
| * xmode = x(imode), where y(xmode) >= y(x(i)) for |
| * all i != imode) |
| * &xvariance (<optional return> variance of x) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) If the histogram represents the relation y(x), the |
| * computed values that are returned are the x values. |
| * These are NOT the bucket indices i; they are related to the |
| * bucket indices by |
| * x(i) = startx + i * deltax |
| */ |
| l_int32 |
| numaGetHistogramStats(NUMA *nahisto, |
| l_float32 startx, |
| l_float32 deltax, |
| l_float32 *pxmean, |
| l_float32 *pxmedian, |
| l_float32 *pxmode, |
| l_float32 *pxvariance) |
| { |
| PROCNAME("numaGetHistogramStats"); |
| |
| if (pxmean) *pxmean = 0.0; |
| if (pxmedian) *pxmedian = 0.0; |
| if (pxmode) *pxmode = 0.0; |
| if (pxvariance) *pxvariance = 0.0; |
| if (!nahisto) |
| return ERROR_INT("nahisto not defined", procName, 1); |
| |
| return numaGetHistogramStatsOnInterval(nahisto, startx, deltax, 0, 0, |
| pxmean, pxmedian, pxmode, |
| pxvariance); |
| } |
| |
| |
| /*! |
| * numaGetHistogramStatsOnInterval() |
| * |
| * Input: nahisto (histogram: y(x(i)), i = 0 ... nbins - 1) |
| * startx (x value of first bin: x(0)) |
| * deltax (x increment between bins; the bin size; x(1) - x(0)) |
| * ifirst (first bin to use for collecting stats) |
| * ilast (last bin for collecting stats; use 0 to go to the end) |
| * &xmean (<optional return> mean value of histogram) |
| * &xmedian (<optional return> median value of histogram) |
| * &xmode (<optional return> mode value of histogram: |
| * xmode = x(imode), where y(xmode) >= y(x(i)) for |
| * all i != imode) |
| * &xvariance (<optional return> variance of x) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) If the histogram represents the relation y(x), the |
| * computed values that are returned are the x values. |
| * These are NOT the bucket indices i; they are related to the |
| * bucket indices by |
| * x(i) = startx + i * deltax |
| */ |
| l_int32 |
| numaGetHistogramStatsOnInterval(NUMA *nahisto, |
| l_float32 startx, |
| l_float32 deltax, |
| l_int32 ifirst, |
| l_int32 ilast, |
| l_float32 *pxmean, |
| l_float32 *pxmedian, |
| l_float32 *pxmode, |
| l_float32 *pxvariance) |
| { |
| l_int32 i, n, imax; |
| l_float32 sum, sumval, halfsum, moment, var, x, y, ymax; |
| |
| PROCNAME("numaGetHistogramStats"); |
| |
| if (pxmean) *pxmean = 0.0; |
| if (pxmedian) *pxmedian = 0.0; |
| if (pxmode) *pxmode = 0.0; |
| if (pxvariance) *pxvariance = 0.0; |
| if (!nahisto) |
| return ERROR_INT("nahisto not defined", procName, 1); |
| if (!pxmean && !pxmedian && !pxmode && !pxvariance) |
| return ERROR_INT("nothing to compute", procName, 1); |
| |
| n = numaGetCount(nahisto); |
| if (ilast <= 0) ilast = n - 1; |
| if (ifirst < 0) ifirst = 0; |
| if (ifirst > ilast || ifirst > n - 1) |
| return ERROR_INT("ifirst is too large", procName, 1); |
| for (sum = 0.0, moment = 0.0, var = 0.0, i = ifirst; i <= ilast ; i++) { |
| x = startx + i * deltax; |
| numaGetFValue(nahisto, i, &y); |
| sum += y; |
| moment += x * y; |
| var += x * x * y; |
| } |
| if (sum == 0.0) |
| return ERROR_INT("sum is 0", procName, 1); |
| |
| if (pxmean) |
| *pxmean = moment / sum; |
| if (pxvariance) |
| *pxvariance = var / sum - moment * moment / (sum * sum); |
| |
| if (pxmedian) { |
| halfsum = sum / 2.0; |
| for (sumval = 0.0, i = ifirst; i <= ilast; i++) { |
| numaGetFValue(nahisto, i, &y); |
| sumval += y; |
| if (sumval >= halfsum) { |
| *pxmedian = startx + i * deltax; |
| break; |
| } |
| } |
| } |
| |
| if (pxmode) { |
| ymax = -1.0e10; |
| for (i = ifirst; i <= ilast; i++) { |
| numaGetFValue(nahisto, i, &y); |
| if (y > ymax) { |
| ymax = y; |
| imax = i; |
| } |
| } |
| *pxmode = startx + imax * deltax; |
| } |
| |
| return 0; |
| } |
| |
| |
| /*! |
| * numaMakeRankFromHistogram() |
| * |
| * Input: startx (xval corresponding to first element in nay) |
| * deltax (x increment between array elements in nay) |
| * nasy (input histogram, assumed equally spaced) |
| * npts (number of points to evaluate rank function) |
| * &nax (<optional return> array of x values in range) |
| * &nay (<return> rank array of specified npts) |
| * Return: 0 if OK, 1 on error |
| */ |
| l_int32 |
| numaMakeRankFromHistogram(l_float32 startx, |
| l_float32 deltax, |
| NUMA *nasy, |
| l_int32 npts, |
| NUMA **pnax, |
| NUMA **pnay) |
| { |
| l_int32 i, n; |
| l_float32 sum, fval; |
| NUMA *nan, *nar; |
| |
| PROCNAME("numaMakeRankFromHistogram"); |
| |
| if (pnax) *pnax = NULL; |
| if (!pnay) |
| return ERROR_INT("&nay not defined", procName, 1); |
| *pnay = NULL; |
| if (!nasy) |
| return ERROR_INT("nasy not defined", procName, 1); |
| if (!pnay) |
| return ERROR_INT("&nay not defined", procName, 1); |
| if ((n = numaGetCount(nasy)) == 0) |
| return ERROR_INT("no bins in nas", procName, 1); |
| |
| /* Normalize and generate the rank array corresponding to |
| * the binned histogram. */ |
| nan = numaNormalizeHistogram(nasy, 1.0); |
| nar = numaCreate(n + 1); /* rank numa corresponding to nan */ |
| sum = 0.0; |
| numaAddNumber(nar, sum); /* first element is 0.0 */ |
| for (i = 0; i < n; i++) { |
| numaGetFValue(nan, i, &fval); |
| sum += fval; |
| numaAddNumber(nar, sum); |
| } |
| |
| /* Compute rank array on full range with specified |
| * number of points and correspondence to x-values. */ |
| numaInterpolateEqxInterval(startx, deltax, nar, L_LINEAR_INTERP, |
| startx, startx + n * deltax, npts, |
| pnax, pnay); |
| numaDestroy(&nan); |
| numaDestroy(&nar); |
| return 0; |
| } |
| |
| |
| /*! |
| * numaHistogramGetRankFromVal() |
| * |
| * Input: na (histogram) |
| * rval (value of input sample for which we want the rank) |
| * &rank (<return> fraction of total samples below rval) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) If we think of the histogram as a function y(x), normalized |
| * to 1, for a given input value of x, this computes the |
| * rank of x, which is the integral of y(x) from the start |
| * value of x to the input value. |
| * (2) This function only makes sense when applied to a Numa that |
| * is a histogram. The values in the histogram can be ints and |
| * floats, and are computed as floats. The rank is returned |
| * as a float between 0.0 and 1.0. |
| * (3) The numa parameters startx and binsize are used to |
| * compute x from the Numa index i. |
| */ |
| l_int32 |
| numaHistogramGetRankFromVal(NUMA *na, |
| l_float32 rval, |
| l_float32 *prank) |
| { |
| l_int32 i, ibinval, n; |
| l_float32 startval, binsize, binval, maxval, fractval, total, sum, val; |
| |
| PROCNAME("numaHistogramGetRankFromVal"); |
| |
| if (!prank) |
| return ERROR_INT("prank not defined", procName, 1); |
| *prank = 0.0; |
| if (!na) |
| return ERROR_INT("na not defined", procName, 1); |
| numaGetXParameters(na, &startval, &binsize); |
| n = numaGetCount(na); |
| if (rval < startval) |
| return 0; |
| maxval = startval + n * binsize; |
| if (rval > maxval) { |
| *prank = 1.0; |
| return 0; |
| } |
| |
| binval = (rval - startval) / binsize; |
| ibinval = (l_int32)binval; |
| if (ibinval >= n) { |
| *prank = 1.0; |
| return 0; |
| } |
| fractval = binval - (l_float32)ibinval; |
| |
| sum = 0.0; |
| for (i = 0; i < ibinval; i++) { |
| numaGetFValue(na, i, &val); |
| sum += val; |
| } |
| numaGetFValue(na, ibinval, &val); |
| sum += fractval * val; |
| numaGetSum(na, &total); |
| *prank = sum / total; |
| |
| /* fprintf(stderr, "binval = %7.3f, rank = %7.3f\n", binval, *prank); */ |
| |
| return 0; |
| } |
| |
| |
| /*! |
| * numaHistogramGetValFromRank() |
| * |
| * Input: na (histogram) |
| * rank (fraction of total samples) |
| * &rval (<return> approx. to the bin value) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) If we think of the histogram as a function y(x), this returns |
| * the value x such that the integral of y(x) from the start |
| * value to x gives the fraction 'rank' of the integral |
| * of y(x) over all bins. |
| * (2) This function only makes sense when applied to a Numa that |
| * is a histogram. The values in the histogram can be ints and |
| * floats, and are computed as floats. The val is returned |
| * as a float, even though the buckets are of integer width. |
| * (3) The numa parameters startx and binsize are used to |
| * compute x from the Numa index i. |
| */ |
| l_int32 |
| numaHistogramGetValFromRank(NUMA *na, |
| l_float32 rank, |
| l_float32 *prval) |
| { |
| l_int32 i, n; |
| l_float32 startval, binsize, rankcount, total, sum, fract, val; |
| |
| PROCNAME("numaHistogramGetValFromRank"); |
| |
| if (!prval) |
| return ERROR_INT("prval not defined", procName, 1); |
| *prval = 0.0; |
| if (!na) |
| return ERROR_INT("na not defined", procName, 1); |
| if (rank < 0.0) { |
| L_WARNING("rank < 0; setting to 0.0", procName); |
| rank = 0.0; |
| } |
| if (rank > 1.0) { |
| L_WARNING("rank > 1.0; setting to 1.0", procName); |
| rank = 1.0; |
| } |
| |
| n = numaGetCount(na); |
| numaGetXParameters(na, &startval, &binsize); |
| numaGetSum(na, &total); |
| rankcount = rank * total; /* count that corresponds to rank */ |
| sum = 0.0; |
| for (i = 0; i < n; i++) { |
| numaGetFValue(na, i, &val); |
| if (sum + val >= rankcount) |
| break; |
| sum += val; |
| } |
| if (val <= 0.0) /* can == 0 if rank == 0.0 */ |
| fract = 0.0; |
| else /* sum + fract * val = rankcount */ |
| fract = (rankcount - sum) / val; |
| |
| /* The use of the fraction of a bin allows a simple calculation |
| * for the histogram value at the given rank. */ |
| *prval = startval + binsize * ((l_float32)i + fract); |
| |
| /* fprintf(stderr, "rank = %7.3f, val = %7.3f\n", rank, *prval); */ |
| |
| return 0; |
| } |
| |
| |
| /*----------------------------------------------------------------------* |
| * Splitting a distribution * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * numaSplitDistribution() |
| * |
| * Input: na (histogram) |
| * scorefract (fraction of the max score, used to determine |
| * the range over which the histogram min is searched) |
| * &splitindex (<optional return> index for splitting) |
| * &ave1 (<optional return> average of lower distribution) |
| * &ave2 (<optional return> average of upper distribution) |
| * &num1 (<optional return> population of lower distribution) |
| * &num2 (<optional return> population of upper distribution) |
| * &nascore (<optional return> for debugging; otherwise use null) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) This function is intended to be used on a distribution of |
| * values that represent two sets, such as a histogram of |
| * pixel values, and the goal is to determine the means of |
| * the two sets and the best splitting point. |
| * (2) The Otsu method finds a split point that divides the distribution |
| * into two parts by maximizing a score function that is the |
| * product of two terms: |
| * (a) the square of the difference of centroids, (ave1 - ave2)^2 |
| * (b) fract1 * (1 - fract1) |
| * where fract1 is the fraction in the lower distribution. |
| * This biases the split point into the larger "bump" (i.e., toward |
| * the point where the (b) term reaches its maximum of 0.25 at |
| * fract1 = 0.5. To avoid this, we define a range of values near |
| * the maximum of the score function, and choose the value within |
| * this range such that the histogram itself has a minimum value. |
| * The range is determined by scorefract: we include all abscissa |
| * values to the left and right of the value that maximizes the |
| * score, such that the score stays above (1 - scorefract) * maxscore. |
| * (3) We normalize the score so that if the two distributions |
| * were of equal size and at opposite ends of the numa, the |
| * score would be 1.0. |
| */ |
| l_int32 |
| numaSplitDistribution(NUMA *na, |
| l_float32 scorefract, |
| l_int32 *psplitindex, |
| l_float32 *pave1, |
| l_float32 *pave2, |
| l_float32 *pnum1, |
| l_float32 *pnum2, |
| NUMA **pnascore) |
| { |
| l_int32 i, n, bestsplit, minrange, maxrange, maxindex; |
| l_float32 ave1, ave2, ave1prev, ave2prev; |
| l_float32 num1, num2, num1prev, num2prev; |
| l_float32 val, minval, sum, fract1; |
| l_float32 norm, score, minscore, maxscore; |
| NUMA *nascore, *naave1, *naave2, *nanum1, *nanum2; |
| |
| PROCNAME("numaSplitDistribution"); |
| |
| if (!na) |
| return ERROR_INT("na not defined", procName, 1); |
| |
| n = numaGetCount(na); |
| if (n <= 1) |
| return ERROR_INT("n = 1 in histogram", procName, 1); |
| numaGetSum(na, &sum); |
| if (sum <= 0.0) |
| return ERROR_INT("sum <= 0.0", procName, 1); |
| norm = 4.0 / ((n - 1) * (n - 1)); |
| ave1prev = 0.0; |
| numaGetHistogramStats(na, 0.0, 1.0, &ave2prev, NULL, NULL, NULL); |
| num1prev = 0.0; |
| num2prev = sum; |
| maxindex = n / 2; /* initialize with something */ |
| |
| /* Split the histogram with [0 ... i] in the lower part |
| * and [i+1 ... n-1] in upper part. First, compute an otsu |
| * score for each possible splitting. */ |
| nascore = numaCreate(n); |
| if (pave2) naave1 = numaCreate(n); |
| if (pave2) naave2 = numaCreate(n); |
| if (pnum1) nanum1 = numaCreate(n); |
| if (pnum2) nanum2 = numaCreate(n); |
| maxscore = 0.0; |
| for (i = 0; i < n; i++) { |
| numaGetFValue(na, i, &val); |
| num1 = num1prev + val; |
| if (num1 == 0) |
| ave1 = ave1prev; |
| else |
| ave1 = (num1prev * ave1prev + i * val) / num1; |
| num2 = num2prev - val; |
| if (num2 == 0) |
| ave2 = ave2prev; |
| else |
| ave2 = (num2prev * ave2prev - i * val) / num2; |
| fract1 = num1 / sum; |
| score = norm * (fract1 * (1 - fract1)) * (ave2 - ave1) * (ave2 - ave1); |
| numaAddNumber(nascore, score); |
| if (pave1) numaAddNumber(naave1, ave1); |
| if (pave2) numaAddNumber(naave2, ave2); |
| if (pnum1) numaAddNumber(nanum1, num1); |
| if (pnum1) numaAddNumber(nanum2, num2); |
| if (score > maxscore) { |
| maxscore = score; |
| maxindex = i; |
| } |
| num1prev = num1; |
| num2prev = num2; |
| ave1prev = ave1; |
| ave2prev = ave2; |
| } |
| |
| /* Next, for all contiguous scores within a specified fraction |
| * of the max, choose the split point as the value with the |
| * minimum in the histogram. */ |
| minscore = (1. - scorefract) * maxscore; |
| for (i = maxindex - 1; i >= 0; i--) { |
| numaGetFValue(nascore, i, &val); |
| if (val < minscore) |
| break; |
| } |
| minrange = i + 1; |
| for (i = maxindex + 1; i < n; i++) { |
| numaGetFValue(nascore, i, &val); |
| if (val < minscore) |
| break; |
| } |
| maxrange = i - 1; |
| numaGetFValue(na, minrange, &minval); |
| bestsplit = minrange; |
| for (i = minrange + 1; i <= maxrange; i++) { |
| numaGetFValue(na, i, &val); |
| if (val < minval) { |
| minval = val; |
| bestsplit = i; |
| } |
| } |
| |
| if (psplitindex) *psplitindex = bestsplit; |
| if (pave1) numaGetFValue(naave1, bestsplit, pave1); |
| if (pave2) numaGetFValue(naave2, bestsplit, pave2); |
| if (pnum1) numaGetFValue(nanum1, bestsplit, pnum1); |
| if (pnum2) numaGetFValue(nanum2, bestsplit, pnum2); |
| |
| if (pnascore) { /* debug mode */ |
| fprintf(stderr, "minrange = %d, maxrange = %d\n", minrange, maxrange); |
| fprintf(stderr, "minval = %10.0f\n", minval); |
| gplotSimple1(nascore, GPLOT_X11, "junkoutroot", |
| "Score for split distribution"); |
| *pnascore = nascore; |
| } |
| else |
| numaDestroy(&nascore); |
| |
| if (pave1) numaDestroy(&naave1); |
| if (pave2) numaDestroy(&naave2); |
| if (pnum1) numaDestroy(&nanum1); |
| if (pnum2) numaDestroy(&nanum2); |
| return 0; |
| } |
| |
| |
| /*----------------------------------------------------------------------* |
| * Extrema finding * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * numaFindPeaks() |
| * |
| * Input: source na |
| * max number of peaks to be found |
| * fract1 (min fraction of peak value) |
| * fract2 (min slope) |
| * Return: peak na, or null on error. |
| * |
| * Notes: |
| * (1) The returned na consists of sets of four numbers representing |
| * the peak, in the following order: |
| * left edge; peak center; right edge; normalized peak area |
| */ |
| NUMA * |
| numaFindPeaks(NUMA *nas, |
| l_int32 nmax, |
| l_float32 fract1, |
| l_float32 fract2) |
| { |
| l_int32 i, k, n, maxloc, lloc, rloc; |
| l_float32 fmaxval, sum, total, newtotal, val, lastval; |
| l_float32 peakfract; |
| NUMA *na, *napeak; |
| |
| PROCNAME("numaFindPeaks"); |
| |
| if (!nas) |
| return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); |
| n = numaGetCount(nas); |
| numaGetSum(nas, &total); |
| |
| /* We munge this copy */ |
| if ((na = numaCopy(nas)) == NULL) |
| return (NUMA *)ERROR_PTR("na not made", procName, NULL); |
| if ((napeak = numaCreate(4 * nmax)) == NULL) |
| return (NUMA *)ERROR_PTR("napeak not made", procName, NULL); |
| |
| for (k = 0; k < nmax; k++) { |
| numaGetSum(na, &newtotal); |
| if (newtotal == 0.0) /* sanity check */ |
| break; |
| numaGetMax(na, &fmaxval, &maxloc); |
| sum = fmaxval; |
| lastval = fmaxval; |
| lloc = 0; |
| for (i = maxloc - 1; i >= 0; --i) { |
| numaGetFValue(na, i, &val); |
| if (val == 0.0) { |
| lloc = i + 1; |
| break; |
| } |
| if (val > fract1 * fmaxval) { |
| sum += val; |
| lastval = val; |
| continue; |
| } |
| if (lastval - val > fract2 * lastval) { |
| sum += val; |
| lastval = val; |
| continue; |
| } |
| lloc = i; |
| break; |
| } |
| lastval = fmaxval; |
| rloc = n - 1; |
| for (i = maxloc + 1; i < n; ++i) { |
| numaGetFValue(na, i, &val); |
| if (val == 0.0) { |
| rloc = i - 1; |
| break; |
| } |
| if (val > fract1 * fmaxval) { |
| sum += val; |
| lastval = val; |
| continue; |
| } |
| if (lastval - val > fract2 * lastval) { |
| sum += val; |
| lastval = val; |
| continue; |
| } |
| rloc = i; |
| break; |
| } |
| peakfract = sum / total; |
| numaAddNumber(napeak, lloc); |
| numaAddNumber(napeak, maxloc); |
| numaAddNumber(napeak, rloc); |
| numaAddNumber(napeak, peakfract); |
| |
| for (i = lloc; i <= rloc; i++) |
| numaSetValue(na, i, 0.0); |
| } |
| |
| numaDestroy(&na); |
| return napeak; |
| } |
| |
| |
| /*! |
| * numaFindExtrema() |
| * |
| * Input: nas (input values) |
| * delta (relative amount to resolve peaks and valleys) |
| * Return: nad (locations of extrema), or null on error |
| * |
| * Notes: |
| * (1) This returns a sequence of extrema (peaks and valleys). |
| * (2) The algorithm is analogous to that for determining |
| * mountain peaks. Suppose we have a local peak, with |
| * bumps on the side. Under what conditions can we consider |
| * those 'bumps' to be actual peaks? The answer: if the |
| * bump is separated from the peak by a saddle that is at |
| * least 500 feet below the bump. |
| * (3) Operationally, suppose we are looking for a peak. |
| * We are keeping the largest value we've seen since the |
| * last valley, and are looking for a value that is delta |
| * BELOW our current peak. When we find such a value, |
| * we label the peak, use the current value to label the |
| * valley, and then do the same operation in reverse (looking |
| * for a valley). |
| */ |
| NUMA * |
| numaFindExtrema(NUMA *nas, |
| l_float32 delta) |
| { |
| l_int32 i, n, found, loc, direction; |
| l_float32 startval, val, maxval, minval; |
| NUMA *nad; |
| |
| PROCNAME("numaFindExtrema"); |
| |
| if (!nas) |
| return (NUMA *)ERROR_PTR("nas not defined", procName, NULL); |
| |
| n = numaGetCount(nas); |
| nad = numaCreate(0); |
| |
| /* We don't know if we'll find a peak or valley first, |
| * but use the first element of nas as the reference point. |
| * Break when we deviate by 'delta' from the first point. */ |
| numaGetFValue(nas, 0, &startval); |
| found = FALSE; |
| for (i = 1; i < n; i++) { |
| numaGetFValue(nas, i, &val); |
| if (L_ABS(val - startval) >= delta) { |
| found = TRUE; |
| break; |
| } |
| } |
| |
| if (!found) |
| return nad; /* it's empty */ |
| |
| /* Are we looking for a peak or a valley? */ |
| if (val > startval) { /* peak */ |
| direction = 1; |
| maxval = val; |
| } |
| else { |
| direction = -1; |
| minval = val; |
| } |
| loc = i; |
| |
| /* Sweep through the rest of the array, recording alternating |
| * peak/valley extrema. */ |
| for (i = i + 1; i < n; i++) { |
| numaGetFValue(nas, i, &val); |
| if (direction == 1 && val > maxval ) { /* new local max */ |
| maxval = val; |
| loc = i; |
| } |
| else if (direction == -1 && val < minval ) { /* new local min */ |
| minval = val; |
| loc = i; |
| } |
| else if (direction == 1 && (maxval - val >= delta)) { |
| numaAddNumber(nad, loc); /* save the current max location */ |
| direction = -1; /* reverse: start looking for a min */ |
| minval = val; |
| loc = i; /* current min location */ |
| } |
| else if (direction == -1 && (val - minval >= delta)) { |
| numaAddNumber(nad, loc); /* save the current min location */ |
| direction = 1; /* reverse: start looking for a max */ |
| maxval = val; |
| loc = i; /* current max location */ |
| } |
| } |
| |
| /* Save the final extremum */ |
| /* numaAddNumber(nad, loc); */ |
| return nad; |
| } |
| |
| |
| /*----------------------------------------------------------------------* |
| * Threshold crossings and frequency analysis * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * numaSelectCrossingThreshold() |
| * |
| * Input: nax (<optional> numa of abscissa values; can be NULL) |
| * nay (signal) |
| * estthresh (estimated pixel threshold for crossing: e.g., for |
| * images, white <--> black; typ. ~120) |
| * &bestthresh (<return> robust estimate of threshold to use) |
| * Return: 0 if OK, 1 on error |
| * |
| * Note: |
| * (1) When a valid threshold is used, the number of crossings is |
| * a maximum, because none are missed. If no threshold intersects |
| * all the crossings, the crossings must be determined with |
| * numaCrossingsByPeaks(). |
| * (2) @estthresh is an input estimate of the threshold that should |
| * be used. We compute the crossings with 41 thresholds |
| * (20 below and 20 above). There is a range in which the |
| * number of crossings is a maximum. Return a threshold |
| * in the center of this stable plateau of crossings. |
| * This can then be used with numaCrossingsByThreshold() |
| * to get a good estimate of crossing locations. |
| */ |
| l_int32 |
| numaSelectCrossingThreshold(NUMA *nax, |
| NUMA *nay, |
| l_float32 estthresh, |
| l_float32 *pbestthresh) |
| { |
| l_int32 i, inrun, istart, iend, maxstart, maxend, runlen, maxrunlen; |
| l_int32 val, maxval, nmax, count; |
| l_float32 thresh, fmaxval, fmodeval; |
| NUMA *nat, *nac; |
| |
| PROCNAME("numaSelectCrossingThreshold"); |
| |
| if (!nay) |
| return ERROR_INT("nay not defined", procName, 1); |
| |
| /* Compute the number of crossings for different thresholds */ |
| nat = numaCreate(41); |
| for (i = 0; i < 41; i++) { |
| thresh = estthresh - 80.0 + 4.0 * i; |
| nac = numaCrossingsByThreshold(nax, nay, thresh); |
| numaAddNumber(nat, numaGetCount(nac)); |
| numaDestroy(&nac); |
| } |
| |
| /* Find the center of the plateau of max crossings, which |
| * extends from thresh[istart] to thresh[iend]. */ |
| numaGetMax(nat, &fmaxval, NULL); |
| maxval = (l_int32)fmaxval; |
| nmax = 0; |
| for (i = 0; i < 41; i++) { |
| numaGetIValue(nat, i, &val); |
| if (val == maxval) |
| nmax++; |
| } |
| if (nmax < 3) { /* likely accidental max; try the mode */ |
| numaGetMode(nat, &fmodeval, &count); |
| if (count > nmax && fmodeval > 0.5 * fmaxval) |
| maxval = (l_int32)fmodeval; /* use the mode */ |
| } |
| |
| inrun = FALSE; |
| iend = 40; |
| maxrunlen = 0; |
| for (i = 0; i < 41; i++) { |
| numaGetIValue(nat, i, &val); |
| if (val == maxval) { |
| if (!inrun) { |
| istart = i; |
| inrun = TRUE; |
| } |
| continue; |
| } |
| if (inrun && (val != maxval)) { |
| iend = i - 1; |
| runlen = iend - istart + 1; |
| inrun = FALSE; |
| if (runlen > maxrunlen) { |
| maxstart = istart; |
| maxend = iend; |
| maxrunlen = runlen; |
| } |
| } |
| } |
| if (inrun) { |
| runlen = i - istart; |
| if (runlen > maxrunlen) { |
| maxstart = istart; |
| maxend = i - 1; |
| maxrunlen = runlen; |
| } |
| } |
| |
| #if 0 |
| foundfirst = FALSE; |
| iend = 40; |
| for (i = 0; i < 41; i++) { |
| numaGetIValue(nat, i, &val); |
| if (val == maxval) { |
| if (!foundfirst) { |
| istart = i; |
| foundfirst = TRUE; |
| } |
| } |
| if ((val != maxval) && foundfirst) { |
| iend = i - 1; |
| break; |
| } |
| } |
| nmax = iend - istart + 1; |
| #endif |
| |
| *pbestthresh = estthresh - 80.0 + 2.0 * (l_float32)(maxstart + maxend); |
| |
| #if DEBUG_CROSSINGS |
| fprintf(stderr, "\nCrossings attain a maximum at %d thresholds, between:\n" |
| " thresh[%d] = %5.1f and thresh[%d] = %5.1f\n", |
| nmax, maxstart, estthresh - 80.0 + 4.0 * maxstart, |
| maxend, estthresh - 80.0 + 4.0 * maxend); |
| fprintf(stderr, "The best choice: %5.1f\n", *pbestthresh); |
| fprintf(stderr, "Number of crossings at the 41 thresholds:"); |
| numaWriteStream(stderr, nat); |
| #endif /* DEBUG_CROSSINGS */ |
| |
| numaDestroy(&nat); |
| return 0; |
| } |
| |
| |
| /*! |
| * numaCrossingsByThreshold() |
| * |
| * Input: nax (<optional> numa of abscissa values; can be NULL) |
| * nay (numa of ordinate values, corresponding to nax) |
| * thresh (threshold value for nay) |
| * Return: nad (abscissa pts at threshold), or null on error |
| * |
| * Notes: |
| * (1) If nax == NULL, we use startx and delx from nay to compute |
| * the crossing values in nad. |
| */ |
| NUMA * |
| numaCrossingsByThreshold(NUMA *nax, |
| NUMA *nay, |
| l_float32 thresh) |
| { |
| l_int32 i, n; |
| l_float32 startx, delx; |
| l_float32 xval1, xval2, yval1, yval2, delta1, delta2, crossval, fract; |
| NUMA *nad; |
| |
| PROCNAME("numaCrossingsByThreshold"); |
| |
| if (!nay) |
| return (NUMA *)ERROR_PTR("nay not defined", procName, NULL); |
| n = numaGetCount(nay); |
| |
| if (nax && (numaGetCount(nax) != n)) |
| return (NUMA *)ERROR_PTR("nax and nay sizes differ", procName, NULL); |
| |
| nad = numaCreate(0); |
| numaGetFValue(nay, 0, &yval1); |
| numaGetXParameters(nay, &startx, &delx); |
| if (nax) |
| numaGetFValue(nax, 0, &xval1); |
| else |
| xval1 = startx; |
| for (i = 1; i < n; i++) { |
| numaGetFValue(nay, i, &yval2); |
| if (nax) |
| numaGetFValue(nax, i, &xval2); |
| else |
| xval2 = startx + i * delx; |
| delta1 = yval1 - thresh; |
| delta2 = yval2 - thresh; |
| if (delta1 == 0.0) |
| numaAddNumber(nad, xval1); |
| else if (delta2 == 0.0) |
| numaAddNumber(nad, xval2); |
| else if (delta1 * delta2 < 0.0) { /* crossing */ |
| fract = L_ABS(delta1) / L_ABS(yval1 - yval2); |
| crossval = xval1 + fract * (xval2 - xval1); |
| numaAddNumber(nad, crossval); |
| } |
| xval1 = xval2; |
| yval1 = yval2; |
| } |
| |
| return nad; |
| } |
| |
| |
| /*! |
| * numaCrossingsByPeaks() |
| * |
| * Input: nax (<optional> numa of abscissa values) |
| * nay (numa of ordinate values, corresponding to nax) |
| * delta (parameter used to identify when a new peak can be found) |
| * Return: nad (abscissa pts at threshold), or null on error |
| * |
| * Notes: |
| * (1) If nax == NULL, we use startx and delx from nay to compute |
| * the crossing values in nad. |
| */ |
| NUMA * |
| numaCrossingsByPeaks(NUMA *nax, |
| NUMA *nay, |
| l_float32 delta) |
| { |
| l_int32 i, j, n, np, previndex, curindex; |
| l_float32 startx, delx; |
| l_float32 xval1, xval2, yval1, yval2, delta1, delta2; |
| l_float32 prevval, curval, thresh, crossval, fract; |
| NUMA *nap, *nad; |
| |
| PROCNAME("numaCrossingsByPeaks"); |
| |
| if (!nax) |
| return (NUMA *)ERROR_PTR("nax not defined", procName, NULL); |
| if (!nay) |
| return (NUMA *)ERROR_PTR("nay not defined", procName, NULL); |
| |
| n = numaGetCount(nax); |
| if (numaGetCount(nay) != n) |
| return (NUMA *)ERROR_PTR("nax and nay sizes differ", procName, NULL); |
| |
| /* Find the extrema. Also add last point in nay to get |
| * the last transition (from the last peak to the end). |
| * The number of crossings is 1 more than the number of extrema. */ |
| nap = numaFindExtrema(nay, delta); |
| numaAddNumber(nap, n - 1); |
| np = numaGetCount(nap); |
| L_INFO_INT("Number of crossings: %d", procName, np); |
| |
| /* Do all computation in index units of nax */ |
| nad = numaCreate(np); /* output crossings, in nax units */ |
| previndex = 0; /* prime the search with 1st point */ |
| numaGetFValue(nay, 0, &prevval); /* prime the search with 1st point */ |
| numaGetXParameters(nay, &startx, &delx); |
| for (i = 0; i < np; i++) { |
| numaGetIValue(nap, i, &curindex); |
| numaGetFValue(nay, curindex, &curval); |
| thresh = (prevval + curval) / 2.0; |
| /* fprintf(stderr, "thresh[%d] = %7.3f\n", i, thresh); */ |
| if (nax) |
| numaGetFValue(nax, previndex, &xval1); |
| else |
| xval1 = startx + previndex * delx; |
| numaGetFValue(nay, previndex, &yval1); |
| for (j = previndex + 1; j <= curindex; j++) { |
| if (nax) |
| numaGetFValue(nax, j, &xval2); |
| else |
| xval2 = startx + j * delx; |
| numaGetFValue(nay, j, &yval2); |
| delta1 = yval1 - thresh; |
| delta2 = yval2 - thresh; |
| if (delta1 == 0.0) { |
| numaAddNumber(nad, xval1); |
| break; |
| } |
| else if (delta2 == 0.0) { |
| numaAddNumber(nad, xval2); |
| break; |
| } |
| else if (delta1 * delta2 < 0.0) { /* crossing */ |
| fract = L_ABS(delta1) / L_ABS(yval1 - yval2); |
| crossval = xval1 + fract * (xval2 - xval1); |
| numaAddNumber(nad, crossval); |
| break; |
| } |
| xval1 = xval2; |
| yval1 = yval2; |
| } |
| previndex = curindex; |
| prevval = curval; |
| } |
| |
| numaDestroy(&nap); |
| return nad; |
| } |
| |
| |
| /*! |
| * numaEvalBestHaarParameters() |
| * |
| * Input: nas (numa of non-negative signal values) |
| * relweight (relative weight of (-1 comb) / (+1 comb) |
| * contributions to the 'convolution'. In effect, |
| * the convolution kernel is a comb consisting of |
| * alternating +1 and -weight.) |
| * nwidth (number of widths to consider) |
| * nshift (number of shifts to consider for each width) |
| * minwidth (smallest width to consider) |
| * maxwidth (largest width to consider) |
| * &bestwidth (<return> width giving largest score) |
| * &bestshift (<return> shift giving largest score) |
| * &bestscore (<optional return> convolution with |
| * "Haar"-like comb) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) This does a linear sweep of widths, evaluating at @nshift |
| * shifts for each width, computing the score from a convolution |
| * with a long comb, and finding the (width, shift) pair that |
| * gives the maximum score. The best width is the "half-wavelength" |
| * of the signal. |
| * (2) The convolving function is a comb of alternating values |
| * +1 and -1 * relweight, separated by the width and phased by |
| * the shift. This is similar to a Haar transform, except |
| * there the convolution is performed with a square wave. |
| * (3) The function is useful for finding the line spacing |
| * and strength of line signal from pixel sum projections. |
| * (4) The score is normalized to the size of nas divided by |
| * the number of half-widths. For image applications, the input is |
| * typically an array of pixel projections, so one should |
| * normalize by dividing the score by the image width in the |
| * pixel projection direction. |
| */ |
| l_int32 |
| numaEvalBestHaarParameters(NUMA *nas, |
| l_float32 relweight, |
| l_int32 nwidth, |
| l_int32 nshift, |
| l_float32 minwidth, |
| l_float32 maxwidth, |
| l_float32 *pbestwidth, |
| l_float32 *pbestshift, |
| l_float32 *pbestscore) |
| { |
| l_int32 i, j; |
| l_float32 delwidth, delshift, width, shift, score; |
| l_float32 bestwidth, bestshift, bestscore; |
| |
| PROCNAME("numaEvalBestHaarParameters"); |
| |
| if (!nas) |
| return ERROR_INT("nas not defined", procName, 1); |
| if (!pbestwidth || !pbestshift) |
| return ERROR_INT("&bestwidth and &bestshift not defined", procName, 1); |
| |
| bestscore = 0.0; |
| delwidth = (maxwidth - minwidth) / (nwidth - 1.0); |
| for (i = 0; i < nwidth; i++) { |
| width = minwidth + delwidth * i; |
| delshift = width / (l_float32)(nshift); |
| for (j = 0; j < nshift; j++) { |
| shift = j * delshift; |
| numaEvalHaarSum(nas, width, shift, relweight, &score); |
| if (score > bestscore) { |
| bestscore = score; |
| bestwidth = width; |
| bestshift = shift; |
| #if DEBUG_FREQUENCY |
| fprintf(stderr, "width = %7.3f, shift = %7.3f, score = %7.3f\n", |
| width, shift, score); |
| #endif /* DEBUG_FREQUENCY */ |
| } |
| } |
| } |
| |
| *pbestwidth = bestwidth; |
| *pbestshift = bestshift; |
| if (pbestscore) |
| *pbestscore = bestscore; |
| return 0; |
| } |
| |
| |
| /*! |
| * numaEvalHaarSum() |
| * |
| * Input: nas (numa of non-negative signal values) |
| * width (distance between +1 and -1 in convolution comb) |
| * shift (phase of the comb: location of first +1) |
| * relweight (relative weight of (-1 comb) / (+1 comb) |
| * contributions to the 'convolution'. In effect, |
| * the convolution kernel is a comb consisting of |
| * alternating +1 and -weight.) |
| * &score (<return> convolution with "Haar"-like comb) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) This does a convolution with a comb of alternating values |
| * +1 and -relweight, separated by the width and phased by the shift. |
| * This is similar to a Haar transform, except that for Haar, |
| * (1) the convolution kernel is symmetric about 0, so the |
| * relweight is 1.0, and |
| * (2) the convolution is performed with a square wave. |
| * (2) The score is normalized to the size of nas divided by |
| * twice the "width". For image applications, the input is |
| * typically an array of pixel projections, so one should |
| * normalize by dividing the score by the image width in the |
| * pixel projection direction. |
| * (3) To get a Haar-like result, use relweight = 1.0. For detecting |
| * signals where you expect every other sample to be close to |
| * zero, as with barcodes or filtered text lines, you can |
| * use relweight > 1.0. |
| */ |
| l_int32 |
| numaEvalHaarSum(NUMA *nas, |
| l_float32 width, |
| l_float32 shift, |
| l_float32 relweight, |
| l_float32 *pscore) |
| { |
| l_int32 i, n, nsamp, index; |
| l_float32 score, weight, val; |
| |
| PROCNAME("numaEvalHaarSum"); |
| |
| if (!pscore) |
| return ERROR_INT("&score not defined", procName, 1); |
| *pscore = 0.0; |
| if (!nas) |
| return ERROR_INT("nas not defined", procName, 1); |
| if ((n = numaGetCount(nas)) < 2 * width) |
| return ERROR_INT("nas size too small", procName, 1); |
| |
| score = 0.0; |
| nsamp = (l_int32)((n - shift) / width); |
| for (i = 0; i < nsamp; i++) { |
| index = (l_int32)(shift + i * width); |
| weight = (i % 2) ? 1.0 : -1.0 * relweight; |
| numaGetFValue(nas, index, &val); |
| score += weight * val; |
| } |
| |
| *pscore = 2.0 * width * score / (l_float32)n; |
| return 0; |
| } |
| |
