0001: /*
0002: * $RCSfile: Histogram.java,v $
0003: *
0004: * Copyright (c) 2005 Sun Microsystems, Inc. All rights reserved.
0005: *
0006: * Use is subject to license terms.
0007: *
0008: * $Revision: 1.1 $
0009: * $Date: 2005/02/11 04:57:08 $
0010: * $State: Exp $
0011: */
0012: package javax.media.jai;
0013:
0014: import java.awt.Rectangle;
0015: import java.awt.image.DataBuffer;
0016: import java.awt.image.Raster;
0017: import java.awt.image.SampleModel;
0018: import java.io.Serializable;
0019: import java.util.LinkedList;
0020: import java.util.ListIterator;
0021:
0022: /**
0023: * This class represents a histogram accumulated from a
0024: * <code>RenderedImage</code>.
0025: *
0026: * <p> A "bin" is a container, where each element stores the total number
0027: * of pixel samples of an image whose
0028: * values lie within a given range. A histogram of an image consists of
0029: * a list of such bins whose range does not overlap with each other
0030: * (mutually exclusive). For an image that has multiple samples per
0031: * pixel (multi-banded images), a separate list of bins represents each
0032: * individual band.
0033: *
0034: * <p> A "low-value" specifies the lowest inclusive pixel value to be
0035: * checked, and a "high-value" specifies the highest exclusive pixel
0036: * value to be checked. Therefore, the width of a bin
0037: * (<code>binWidth</code>) is determined by
0038: * <code>(highValue - lowValue) / numberOfBins</code>. The range
0039: * of bin <code>i</code> is defined as from
0040: * <code>lowValue + i * binWidth</code> inclusive to
0041: * <code>lowValue + (i + 1) * binWidth</code> exclusive.
0042: *
0043: * <p> The image may have any data type. Its histogram may be accumulated
0044: * over the entire image, or over a specific region-of-interest (ROI)
0045: * within the image's bounds. Furthermore, the horizontal and vertical
0046: * subsampling factors specify the rate of sampling in the two directions,
0047: * so that only every <i>n</i>th pixel will be counted. This allows the
0048: * accuracy of the histogram to be traded for the speed of the computation.
0049: * Of course a subsampling rate of 1 means every pixel will be counted.
0050: *
0051: * <p> The "Histogram" operator generates the histogram data of an image
0052: * and uses this object to store the final pixel counts. The operator
0053: * returns an instance of this class when a request is made via the
0054: * <code>getProperty</code> method for the "histogram" property. The actual
0055: * bins may be obtained by calling the <code>getBins</code> method.
0056: *
0057: * @see ROI
0058: * @see javax.media.jai.operator.HistogramDescriptor
0059: *
0060: */
0061: public class Histogram implements Serializable {
0062:
0063: /** The number of bins used for each band of the image. */
0064: private int[] numBins;
0065:
0066: /**
0067: * The lowest inclusive pixel value of the image checked for each band.
0068: */
0069: private double[] lowValue;
0070:
0071: /**
0072: * The highest exclusive pixel value of the image checked for each band.
0073: */
0074: private double[] highValue;
0075:
0076: /**
0077: * The number of bands of the image from which the histogram is
0078: * accumulated. This is the same as the number of bands of the
0079: * bins (bins.length).
0080: */
0081: private int numBands;
0082:
0083: /** The width of a bin for each band. */
0084: private double[] binWidth;
0085:
0086: /** The bins for each band, used to hold the pixel counts. */
0087: private int[][] bins = null;
0088:
0089: /** The total bin count over all bins for each band. */
0090: private int[] totals = null;
0091:
0092: /** The mean value over all bins for each band. */
0093: private double[] mean = null;
0094:
0095: /**
0096: * Copy an int array into a new int array of a given length padding
0097: * with zeroth element if needed.
0098: *
0099: * @throws IllegalArgumentException If <code>array</code>
0100: * is <code>null</code> or its length is <code>0</code>.
0101: */
0102: private static final int[] fill(int[] array, int newLength) {
0103: int[] newArray = null;
0104:
0105: if (array == null || array.length == 0) {
0106: throw new IllegalArgumentException(JaiI18N
0107: .getString("Generic0"));
0108: } else if (newLength > 0) {
0109: newArray = new int[newLength];
0110: int oldLength = array.length;
0111: for (int i = 0; i < newLength; i++) {
0112: if (i < oldLength) {
0113: newArray[i] = array[i];
0114: } else {
0115: newArray[i] = array[0];
0116: }
0117: }
0118: }
0119: return newArray;
0120: }
0121:
0122: /**
0123: * Copy an double array into a new double array of a given length padding
0124: * with zeroth element if needed.
0125: *
0126: * @throws IllegalArgumentException If <code>array</code>
0127: * is <code>null</code> or its length is <code>0</code>.
0128: */
0129: private static final double[] fill(double[] array, int newLength) {
0130: double[] newArray = null;
0131:
0132: if (array == null || array.length == 0) {
0133: throw new IllegalArgumentException(JaiI18N
0134: .getString("Generic0"));
0135: } else if (newLength > 0) {
0136: newArray = new double[newLength];
0137: int oldLength = array.length;
0138: for (int i = 0; i < newLength; i++) {
0139: if (i < oldLength) {
0140: newArray[i] = array[i];
0141: } else {
0142: newArray[i] = array[0];
0143: }
0144: }
0145: }
0146: return newArray;
0147: }
0148:
0149: /**
0150: * Constructor.
0151: *
0152: * <p> This constructor should be used when <code>numBins</code>,
0153: * <code>lowValue</code>, and/or <code>highValues</code> are
0154: * different for each band of the image. The length of the arrays
0155: * indicates the number of bands the image has, and the three arrays
0156: * must have the same length.
0157: *
0158: * <p> Since this constructor has no way of knowing the actual number
0159: * of bands of the image, the length of the arrays is not checked
0160: * against anything. Therefore, it is very important that the caller
0161: * supplies the correct arrays of correct lengths, or errors will occur.
0162: *
0163: * @param numBins The number of bins for each band of the image.
0164: * The length of this array indicates the number of bands for
0165: * this histogram.
0166: * @param lowValue The lowest inclusive pixel value checked for
0167: * each band.
0168: * @param highValue The highest exclusive pixel value checked for
0169: * each band.
0170: *
0171: * @throws IllegalArgumentException If <code>numBins</code>,
0172: * <code>lowValue</code>, or <code>highValue</code> is
0173: * <code>null</code>.
0174: * @throws IllegalArgumentException If the array length of the three
0175: * arguments are not the same, or any array length is 0.
0176: * @throws IllegalArgumentException If the number of bins for any band
0177: * is less than or equal to 0.
0178: * @throws IllegalArgumentException If the low-value of any band is
0179: * greater than or equal to its corresponding high-value.
0180: */
0181: public Histogram(int[] numBins, double[] lowValue,
0182: double[] highValue) {
0183:
0184: if (numBins == null || lowValue == null || highValue == null) {
0185: throw new IllegalArgumentException(JaiI18N
0186: .getString("Generic0"));
0187: }
0188:
0189: numBands = numBins.length;
0190:
0191: if (lowValue.length != numBands || highValue.length != numBands) {
0192: throw new IllegalArgumentException(JaiI18N
0193: .getString("Histogram0"));
0194: }
0195:
0196: if (numBands == 0) {
0197: throw new IllegalArgumentException(JaiI18N
0198: .getString("Histogram1"));
0199: }
0200:
0201: for (int i = 0; i < numBands; i++) {
0202: if (numBins[i] <= 0) {
0203: throw new IllegalArgumentException(JaiI18N
0204: .getString("Histogram2"));
0205: }
0206:
0207: if (lowValue[i] >= highValue[i]) {
0208: throw new IllegalArgumentException(JaiI18N
0209: .getString("Histogram3"));
0210: }
0211: }
0212:
0213: this .numBins = (int[]) numBins.clone();
0214: this .lowValue = (double[]) lowValue.clone();
0215: this .highValue = (double[]) highValue.clone();
0216:
0217: binWidth = new double[numBands];
0218:
0219: // Compute binWidth for all bands.
0220: for (int i = 0; i < numBands; i++) {
0221: binWidth[i] = (highValue[i] - lowValue[i]) / numBins[i];
0222: }
0223: }
0224:
0225: /**
0226: * Constructor.
0227: *
0228: * <p> This constructor should be used when <code>numBins</code>,
0229: * <code>lowValue</code>, and/or <code>highValues</code> may be
0230: * different for each band of the image. The number of bands in the
0231: * image is provided explicitly. If any of the arrays provided
0232: * has a length which is less than the number of bands, the first
0233: * element in that array is used to fill out the array to a length
0234: * of <code>numBands</code>.
0235: *
0236: *
0237: * @param numBins The number of bins for each band of the image.
0238: * @param lowValue The lowest inclusive pixel value checked for
0239: * each band.
0240: * @param highValue The highest exclusive pixel value checked for
0241: * each band.
0242: * @param numBands The number of bands in the image.
0243: *
0244: * @throws IllegalArgumentException If <code>numBins</code>,
0245: * <code>lowValue</code>, or <code>highValue</code> is
0246: * <code>null</code>.
0247: * @throws IllegalArgumentException If any array length is 0.
0248: * @throws IllegalArgumentException If the number of bins for any band
0249: * is less than or equal to 0.
0250: * @throws IllegalArgumentException If the low-value of any band is
0251: * greater than or equal to its corresponding high-value.
0252: * @throws IllegalArgumentException If <code>numBands</code> is less
0253: * than or equal to 0.
0254: *
0255: * @since JAI 1.1
0256: */
0257: public Histogram(int[] numBins, double[] lowValue,
0258: double[] highValue, int numBands) {
0259: this (fill(numBins, numBands), fill(lowValue, numBands), fill(
0260: highValue, numBands));
0261: }
0262:
0263: /**
0264: * Constructor.
0265: *
0266: * <p> The same <code>numBins</code>, <code>lowValue</code>, and
0267: * <code>highValue</code> is applied to every band of the image.
0268: *
0269: * @param numBins The number of bins for all bands of the image.
0270: * @param lowValue The lowest inclusive pixel value checked for
0271: * all bands.
0272: * @param highValue The highest exclusive pixel value checked for
0273: * all bands.
0274: * @param numBands The number of bands of the image.
0275: *
0276: * @throws IllegalArgumentException If <code>numBins</code> or
0277: * <code>numBands</code> is less than or equal to 0.
0278: * @throws IllegalArgumentException If <code>lowValue</code>
0279: * is less than or equal to <code>highValue</code>.
0280: *
0281: * @since JAI 1.1
0282: */
0283: public Histogram(int numBins, double lowValue, double highValue,
0284: int numBands) {
0285: if (numBands <= 0) {
0286: throw new IllegalArgumentException(JaiI18N
0287: .getString("Histogram1"));
0288: }
0289:
0290: if (numBins <= 0) {
0291: throw new IllegalArgumentException(JaiI18N
0292: .getString("Histogram2"));
0293: }
0294:
0295: if (lowValue >= highValue) {
0296: throw new IllegalArgumentException(JaiI18N
0297: .getString("Histogram3"));
0298: }
0299:
0300: this .numBands = numBands;
0301: this .numBins = new int[numBands];
0302: this .lowValue = new double[numBands];
0303: this .highValue = new double[numBands];
0304: this .binWidth = new double[numBands];
0305:
0306: double bw = (highValue - lowValue) / numBins; // binWidth
0307:
0308: for (int i = 0; i < numBands; i++) {
0309: this .numBins[i] = numBins;
0310: this .lowValue[i] = lowValue;
0311: this .highValue[i] = highValue;
0312: this .binWidth[i] = bw;
0313: }
0314: }
0315:
0316: /** Returns the number of bins of the histogram for all bands. */
0317: public int[] getNumBins() {
0318: return (int[]) numBins.clone();
0319: }
0320:
0321: /**
0322: * Returns the number of bins of the histogram for a specific band.
0323: *
0324: * @param band The index of the band whose <code>numBins</code>
0325: * is to be returned.
0326: *
0327: * @throws ArrayIndexOutOfBoundsException If an invalid band index
0328: * is specified.
0329: */
0330: public int getNumBins(int band) {
0331: return numBins[band];
0332: }
0333:
0334: /** Returns the lowest inclusive pixel value checked for all bands. */
0335: public double[] getLowValue() {
0336: return (double[]) lowValue.clone();
0337: }
0338:
0339: /**
0340: * Returns the lowest inclusive pixel value checked for a specific band.
0341: *
0342: * @param band The index of the band whose <code>lowValue</code>
0343: * is to be returned.
0344: *
0345: * @throws ArrayIndexOutOfBoundsException If an invalid band index
0346: * is specified.
0347: */
0348: public double getLowValue(int band) {
0349: return lowValue[band];
0350: }
0351:
0352: /** Returns the highest exclusive pixel value checked for all bands. */
0353: public double[] getHighValue() {
0354: return (double[]) highValue.clone();
0355: }
0356:
0357: /**
0358: * Returns the highest exclusive pixel value checked for a specific band.
0359: *
0360: * @param band The index of the band whose <code>highValue</code>
0361: * is to be returned.
0362: *
0363: * @throws ArrayIndexOutOfBoundsException If an invalid band index
0364: * is specified.
0365: */
0366: public double getHighValue(int band) {
0367: return highValue[band];
0368: }
0369:
0370: /**
0371: * Returns the number of bands of the histogram.
0372: * This value is the same as the number of bands of the
0373: * bins, <code>bins.length</code>.
0374: */
0375: public int getNumBands() {
0376: return numBands;
0377: }
0378:
0379: /**
0380: * Returns the array of bands of bins, each bin is the histogram for a
0381: * band, i.e., the format of the returned array is
0382: * <code>int[bands][bins]</code>.
0383: */
0384: public synchronized int[][] getBins() {
0385: if (bins == null) {
0386: bins = new int[numBands][];
0387:
0388: for (int i = 0; i < numBands; i++) {
0389: bins[i] = new int[numBins[i]];
0390: }
0391: }
0392:
0393: return bins;
0394: }
0395:
0396: /**
0397: * Returns the bins of the histogram for a specific band by reference.
0398: *
0399: * @param band The index of the band whose <code>bins</code>
0400: * are to be returned.
0401: *
0402: * @throws ArrayIndexOutOfBoundsException If an invalid band index
0403: * is specified.
0404: */
0405: public int[] getBins(int band) {
0406: getBins(); // make sure bins is initialized
0407: return bins[band];
0408: }
0409:
0410: /**
0411: * Returns the number of pixel samples found in a given bin for a
0412: * specific band.
0413: *
0414: * @param band The index of the band-of-interest.
0415: * @param bin The index of the bin whose value is to be returned.
0416: *
0417: * @throws ArrayIndexOutOfBoundsException If an invalid band or
0418: * bin index is specified.
0419: */
0420: public int getBinSize(int band, int bin) {
0421: getBins(); // make sure bins is initialized
0422: return bins[band][bin];
0423: }
0424:
0425: /**
0426: * Returns the lowest inclusive pixel value of a given bin
0427: * for a specific band.
0428: *
0429: * @param band The index of the band-of-interest.
0430: * @param bin The index of the bin whose <code>lowValue</code>
0431: * is to be returned.
0432: *
0433: * @throws ArrayIndexOutOfBoundsException If an invalid band
0434: * index is specified.
0435: */
0436: public double getBinLowValue(int band, int bin) {
0437: return lowValue[band] + bin * binWidth[band];
0438: }
0439:
0440: /**
0441: * Resets the values of all bins to zero. If <code>bins</code> has
0442: * not been initialized (<code>null</code>), this method does nothing.
0443: */
0444: public void clearHistogram() {
0445: if (bins != null) {
0446: synchronized (bins) {
0447: for (int i = 0; i < numBands; i++) {
0448: int[] b = bins[i];
0449: int length = b.length;
0450:
0451: for (int j = 0; j < length; j++) {
0452: b[j] = 0;
0453: }
0454: }
0455: }
0456: }
0457: }
0458:
0459: /**
0460: * Returns the total bin count over all bins for all bands.
0461: *
0462: * <p> An array which stores the total bin count is kept in this class
0463: * and a reference to this array is returned by this method for
0464: * performance reasons. The elements of the returned array should
0465: * not be modified or undefined errors may occur. The array format is
0466: * <code>int[numBands]</code>.
0467: *
0468: * @since JAI 1.1
0469: */
0470: public int[] getTotals() {
0471: if (totals == null) {
0472: getBins(); // make sure bins is initialized
0473:
0474: synchronized (this ) {
0475: totals = new int[numBands];
0476:
0477: for (int i = 0; i < numBands; i++) {
0478: int[] b = bins[i];
0479: int length = b.length;
0480: int t = 0;
0481:
0482: for (int j = 0; j < length; j++) {
0483: t += b[j];
0484: }
0485:
0486: totals[i] = t;
0487: }
0488: }
0489: }
0490:
0491: return totals;
0492: }
0493:
0494: /**
0495: * Returns the total bin count for the specified sub-range (via "min"
0496: * and "max" bin) of the indicated band. The sub-ragne must fall
0497: * within the actual range of the bins.
0498: *
0499: * @param band The index of the band-of-interest.
0500: * @param minBin The minimum bin index to be counted.
0501: * @param maxBin The maximum bin index to be counted.
0502: *
0503: * @throws ArrayIndexOutOfBoundsException If an invalid band index
0504: * is specified.
0505: * @throws IllegalArgumentException If <code>minBin</code> is greater than
0506: * <code>maxBin</code>.
0507: *
0508: * @since JAI 1.1
0509: */
0510: public int getSubTotal(int band, int minBin, int maxBin) {
0511: if (minBin < 0 || maxBin >= numBins[band]) {
0512: throw new ArrayIndexOutOfBoundsException(JaiI18N
0513: .getString("Histogram5"));
0514: }
0515:
0516: if (minBin > maxBin) {
0517: throw new IllegalArgumentException(JaiI18N
0518: .getString("Histogram10"));
0519: }
0520:
0521: int[] b = getBins(band);
0522: int total = 0;
0523:
0524: for (int i = minBin; i <= maxBin; i++) {
0525: total += b[i];
0526: }
0527:
0528: return total;
0529: }
0530:
0531: /**
0532: * Returns the mean values for all bands of the histogram.
0533: *
0534: * @since JAI 1.1
0535: */
0536: public double[] getMean() {
0537: if (mean == null) {
0538: getTotals(); // make sure totals is computed
0539:
0540: synchronized (this ) {
0541: mean = new double[numBands];
0542:
0543: for (int i = 0; i < numBands; i++) {
0544: int[] counts = getBins(i);
0545: int nBins = numBins[i];
0546: double level = getLowValue(i);
0547: double bw = binWidth[i];
0548:
0549: double mu = 0.0;
0550: double total = totals[i];
0551:
0552: for (int b = 0; b < nBins; b++) {
0553: mu += (counts[b] / total) * level;
0554: level += bw;
0555: }
0556:
0557: mean[i] = mu;
0558: }
0559: }
0560: }
0561:
0562: return mean;
0563: }
0564:
0565: /**
0566: * Accumulates the histogram of the pixels within a specific
0567: * region-of-interest (ROI) by counting them based on the
0568: * indicated horizontal and vertical sampling period. The result
0569: * is stored in the <code>bins</code> array and may be obtained
0570: * by calling the <code>getBins</code> method.
0571: *
0572: * <p> The <code>ROI</code> specifies the region within which
0573: * the pixels are counted. If it is <code>null</code>, the
0574: * entire <code>Raster</code> is counted. If it is not
0575: * <code>null</code> and does not intersect with the
0576: * <code>Raster</code>, then this method returns without changing
0577: * the <code>bins</code>.
0578: *
0579: * <p> The set of pixels to be counted may be obtained by
0580: * intersecting the grid <code>(xStart + i * xPeriod,
0581: * yStart + j * yPeriod); i, j >= 0</code> with the ROI
0582: * and the bounding rectangle of the <code>Raster</code>.
0583: *
0584: * @param raster The Raster that contains the pixels to be counted.
0585: * @param roi The region-of-interest within which the pixels are counted.
0586: * @param xStart The initial X sample coordinate.
0587: * @param yStart The initial Y sample coordinate.
0588: * @param xPeriod The X sampling period.
0589: * @param yPeriod The Y sampling period.
0590: *
0591: * @throws IllegalArgumentException If <code>raster</code> is
0592: * <code>null</code>.
0593: * @throws IllegalArgumentException If the pixels stored in the
0594: * <code>raster</code> do not have the same number of bands
0595: * (samples per pixel) as this histogram's bins.
0596: * @thows RuntimeException if the data type is not supported
0597: * (not in DataBuffer.TYPE_BYTE,..., DataBuff.TYPE_DOUBLE.
0598: */
0599: public void countPixels(Raster raster, ROI roi, int xStart,
0600: int yStart, int xPeriod, int yPeriod) {
0601:
0602: if (raster == null) {
0603: throw new IllegalArgumentException(JaiI18N
0604: .getString("Generic0"));
0605: }
0606:
0607: SampleModel sampleModel = raster.getSampleModel();
0608:
0609: if (sampleModel.getNumBands() != numBands) {
0610: throw new IllegalArgumentException(JaiI18N
0611: .getString("Histogram4"));
0612: }
0613:
0614: Rectangle bounds = raster.getBounds();
0615:
0616: LinkedList rectList;
0617: if (roi == null) { // ROI is the whole Raster
0618: rectList = new LinkedList();
0619: rectList.addLast(bounds);
0620: } else {
0621: rectList = roi.getAsRectangleList(bounds.x, bounds.y,
0622: bounds.width, bounds.height);
0623: if (rectList == null) {
0624: return; // ROI does not intersect with Raster boundary.
0625: }
0626: }
0627:
0628: PixelAccessor accessor = new PixelAccessor(sampleModel, null);
0629:
0630: ListIterator iterator = rectList.listIterator(0);
0631:
0632: while (iterator.hasNext()) {
0633: Rectangle r = (Rectangle) iterator.next();
0634: int tx = r.x;
0635: int ty = r.y;
0636:
0637: // Find the actual ROI based on start and period.
0638: r.x = startPosition(tx, xStart, xPeriod);
0639: r.y = startPosition(ty, yStart, yPeriod);
0640: r.width = tx + r.width - r.x;
0641: r.height = ty + r.height - r.y;
0642:
0643: if (r.width <= 0 || r.height <= 0) {
0644: continue; // no pixel to count in this rectangle
0645: }
0646:
0647: switch (accessor.sampleType) {
0648: case PixelAccessor.TYPE_BIT:
0649: case DataBuffer.TYPE_BYTE:
0650: countPixelsByte(accessor, raster, r, xPeriod, yPeriod);
0651: break;
0652: case DataBuffer.TYPE_USHORT:
0653: countPixelsUShort(accessor, raster, r, xPeriod, yPeriod);
0654: break;
0655: case DataBuffer.TYPE_SHORT:
0656: countPixelsShort(accessor, raster, r, xPeriod, yPeriod);
0657: break;
0658: case DataBuffer.TYPE_INT:
0659: countPixelsInt(accessor, raster, r, xPeriod, yPeriod);
0660: break;
0661: case DataBuffer.TYPE_FLOAT:
0662: countPixelsFloat(accessor, raster, r, xPeriod, yPeriod);
0663: break;
0664: case DataBuffer.TYPE_DOUBLE:
0665: countPixelsDouble(accessor, raster, r, xPeriod, yPeriod);
0666: break;
0667: default:
0668: throw new RuntimeException(JaiI18N
0669: .getString("Histogram11"));
0670: }
0671: }
0672: }
0673:
0674: private void countPixelsByte(PixelAccessor accessor, Raster raster,
0675: Rectangle rect, int xPeriod, int yPeriod) {
0676: UnpackedImageData uid = accessor.getPixels(raster, rect,
0677: DataBuffer.TYPE_BYTE, false);
0678:
0679: byte[][] byteData = uid.getByteData();
0680: int pixelStride = uid.pixelStride * xPeriod;
0681: int lineStride = uid.lineStride * yPeriod;
0682: int[] offsets = uid.bandOffsets;
0683:
0684: for (int b = 0; b < numBands; b++) {
0685: byte[] data = byteData[b];
0686: int lineOffset = offsets[b]; // line offset
0687:
0688: int[] bin = new int[numBins[b]];
0689: double low = lowValue[b];
0690: double high = highValue[b];
0691: double bwidth = binWidth[b];
0692:
0693: for (int h = 0; h < rect.height; h += yPeriod) {
0694: int pixelOffset = lineOffset; // pixel offset
0695: lineOffset += lineStride;
0696:
0697: for (int w = 0; w < rect.width; w += xPeriod) {
0698: int d = data[pixelOffset] & 0xff;
0699: pixelOffset += pixelStride;
0700:
0701: if (d >= low && d < high) {
0702: int i = (int) ((d - low) / bwidth);
0703: bin[i]++;
0704: }
0705: }
0706: }
0707:
0708: mergeBins(b, bin); // merge this band to the whole bins
0709: }
0710: }
0711:
0712: private void countPixelsUShort(PixelAccessor accessor,
0713: Raster raster, Rectangle rect, int xPeriod, int yPeriod) {
0714: UnpackedImageData uid = accessor.getPixels(raster, rect,
0715: DataBuffer.TYPE_USHORT, false);
0716:
0717: short[][] shortData = uid.getShortData();
0718: int pixelStride = uid.pixelStride * xPeriod;
0719: int lineStride = uid.lineStride * yPeriod;
0720: int[] offsets = uid.bandOffsets;
0721:
0722: for (int b = 0; b < numBands; b++) {
0723: short[] data = shortData[b];
0724: int lineOffset = offsets[b]; // line offset
0725:
0726: int[] bin = new int[numBins[b]];
0727: double low = lowValue[b];
0728: double high = highValue[b];
0729: double bwidth = binWidth[b];
0730:
0731: for (int h = 0; h < rect.height; h += yPeriod) {
0732: int pixelOffset = lineOffset; // pixel offset
0733: lineOffset += lineStride;
0734:
0735: for (int w = 0; w < rect.width; w += xPeriod) {
0736: int d = data[pixelOffset] & 0xffff;
0737: pixelOffset += pixelStride;
0738:
0739: if (d >= low && d < high) {
0740: int i = (int) ((d - low) / bwidth);
0741: bin[i]++;
0742: }
0743: }
0744: }
0745:
0746: mergeBins(b, bin); // merge this band to the whole bins
0747: }
0748: }
0749:
0750: private void countPixelsShort(PixelAccessor accessor,
0751: Raster raster, Rectangle rect, int xPeriod, int yPeriod) {
0752: UnpackedImageData uid = accessor.getPixels(raster, rect,
0753: DataBuffer.TYPE_SHORT, false);
0754:
0755: short[][] shortData = uid.getShortData();
0756: int pixelStride = uid.pixelStride * xPeriod;
0757: int lineStride = uid.lineStride * yPeriod;
0758: int[] offsets = uid.bandOffsets;
0759:
0760: for (int b = 0; b < numBands; b++) {
0761: short[] data = shortData[b];
0762: int lineOffset = offsets[b]; // line offset
0763:
0764: int[] bin = new int[numBins[b]];
0765: double low = lowValue[b];
0766: double high = highValue[b];
0767: double bwidth = binWidth[b];
0768:
0769: for (int h = 0; h < rect.height; h += yPeriod) {
0770: int pixelOffset = lineOffset; // pixel offset
0771: lineOffset += lineStride;
0772:
0773: for (int w = 0; w < rect.width; w += xPeriod) {
0774: int d = data[pixelOffset];
0775: pixelOffset += pixelStride;
0776:
0777: if (d >= low && d < high) {
0778: int i = (int) ((d - low) / bwidth);
0779: bin[i]++;
0780: }
0781: }
0782: }
0783:
0784: mergeBins(b, bin); // merge this band to the whole bins
0785: }
0786: }
0787:
0788: private void countPixelsInt(PixelAccessor accessor, Raster raster,
0789: Rectangle rect, int xPeriod, int yPeriod) {
0790: UnpackedImageData uid = accessor.getPixels(raster, rect,
0791: DataBuffer.TYPE_INT, false);
0792:
0793: int[][] intData = uid.getIntData();
0794: int pixelStride = uid.pixelStride * xPeriod;
0795: int lineStride = uid.lineStride * yPeriod;
0796: int[] offsets = uid.bandOffsets;
0797:
0798: for (int b = 0; b < numBands; b++) {
0799: int[] data = intData[b];
0800: int lineOffset = offsets[b]; // line offset
0801:
0802: int[] bin = new int[numBins[b]];
0803: double low = lowValue[b];
0804: double high = highValue[b];
0805: double bwidth = binWidth[b];
0806:
0807: for (int h = 0; h < rect.height; h += yPeriod) {
0808: int pixelOffset = lineOffset; // pixel offset
0809: lineOffset += lineStride;
0810:
0811: for (int w = 0; w < rect.width; w += xPeriod) {
0812: int d = data[pixelOffset];
0813: pixelOffset += pixelStride;
0814:
0815: if (d >= low && d < high) {
0816: int i = (int) ((d - low) / bwidth);
0817: bin[i]++;
0818: }
0819: }
0820: }
0821:
0822: mergeBins(b, bin); // merge this band to the whole bins
0823: }
0824: }
0825:
0826: private void countPixelsFloat(PixelAccessor accessor,
0827: Raster raster, Rectangle rect, int xPeriod, int yPeriod) {
0828: UnpackedImageData uid = accessor.getPixels(raster, rect,
0829: DataBuffer.TYPE_FLOAT, false);
0830:
0831: float[][] floatData = uid.getFloatData();
0832: int pixelStride = uid.pixelStride * xPeriod;
0833: int lineStride = uid.lineStride * yPeriod;
0834: int[] offsets = uid.bandOffsets;
0835:
0836: for (int b = 0; b < numBands; b++) {
0837: float[] data = floatData[b];
0838: int lineOffset = offsets[b]; // line offset
0839:
0840: int[] bin = new int[numBins[b]];
0841: double low = lowValue[b];
0842: double high = highValue[b];
0843: double bwidth = binWidth[b];
0844:
0845: for (int h = 0; h < rect.height; h += yPeriod) {
0846: int pixelOffset = lineOffset; // pixel offset
0847: lineOffset += lineStride;
0848:
0849: for (int w = 0; w < rect.width; w += xPeriod) {
0850: float d = data[pixelOffset];
0851: pixelOffset += pixelStride;
0852:
0853: if (d >= low && d < high) {
0854: int i = (int) ((d - low) / bwidth);
0855: bin[i]++;
0856: }
0857: }
0858: }
0859:
0860: mergeBins(b, bin); // merge this band to the whole bins
0861: }
0862: }
0863:
0864: private void countPixelsDouble(PixelAccessor accessor,
0865: Raster raster, Rectangle rect, int xPeriod, int yPeriod) {
0866: UnpackedImageData uid = accessor.getPixels(raster, rect,
0867: DataBuffer.TYPE_DOUBLE, false);
0868:
0869: double[][] doubleData = uid.getDoubleData();
0870: int pixelStride = uid.pixelStride * xPeriod;
0871: int lineStride = uid.lineStride * yPeriod;
0872: int[] offsets = uid.bandOffsets;
0873:
0874: for (int b = 0; b < numBands; b++) {
0875: double[] data = doubleData[b];
0876: int lineOffset = offsets[b]; // line offset
0877:
0878: int[] bin = new int[numBins[b]];
0879: double low = lowValue[b];
0880: double high = highValue[b];
0881: double bwidth = binWidth[b];
0882:
0883: for (int h = 0; h < rect.height; h += yPeriod) {
0884: int pixelOffset = lineOffset; // pixel offset
0885: lineOffset += lineStride;
0886:
0887: for (int w = 0; w < rect.width; w += xPeriod) {
0888: double d = data[pixelOffset];
0889: pixelOffset += pixelStride;
0890:
0891: if (d >= low && d < high) {
0892: int i = (int) ((d - low) / bwidth);
0893: bin[i]++;
0894: }
0895: }
0896: }
0897:
0898: mergeBins(b, bin); // merge this band to the whole bins
0899: }
0900: }
0901:
0902: /** Finds the first pixel at or after <code>pos</code> to be counted. */
0903: private int startPosition(int pos, int start, int Period) {
0904: int t = (pos - start) % Period;
0905: return t == 0 ? pos : pos + (Period - t);
0906: }
0907:
0908: /** Merges bin count for a band. Synchronized on bins for MT-safe. */
0909: private void mergeBins(int band, int[] bin) {
0910: getBins();
0911:
0912: synchronized (bins) {
0913: int[] b = bins[band];
0914: int length = b.length;
0915:
0916: for (int i = 0; i < length; i++) {
0917: b[i] += bin[i];
0918: }
0919: }
0920: }
0921:
0922: /**
0923: * Returns a specified (absolute) (central) moment of the histogram.
0924: *
0925: * <p> The <i>i</i>th moment in each band is defined to be the mean
0926: * of the image pixel values raised to the <i>i</i>th power in that
0927: * band. For central moments the average of the <i>i</i>th power
0928: * of the deviation from the mean is used. For absolute moments the
0929: * absolute value of the exponentiated term is used.
0930: *
0931: * <p> Note that the mean is the first moment, the average energy the
0932: * second central moment, etc.
0933: *
0934: * @param moment The moment number or index which must be positive or
0935: * an <code>IllegalArgumentException</code> will be thrown.
0936: * @param isAbsolute Whether to calculate the absolute moment.
0937: * @param isCentral Whether to calculate the central moment.
0938: * @return The requested (absolute) (central) moment of the histogram.
0939: *
0940: * @since JAI 1.1
0941: */
0942: public double[] getMoment(int moment, boolean isAbsolute,
0943: boolean isCentral) {
0944: // Check for non-positive moment number.
0945: if (moment < 1) {
0946: throw new IllegalArgumentException(JaiI18N
0947: .getString("Histogram6"));
0948: }
0949:
0950: // If the mean is required but has not yet been calculated
0951: // then calculate it first.
0952: if ((moment == 1 || isCentral) && mean == null) {
0953: getMean();
0954: }
0955:
0956: // If it's the first non-absolute, non-central moment return the mean.
0957: if ((moment == 1) && !isAbsolute && !isCentral) {
0958: return mean;
0959: }
0960:
0961: double[] moments = new double[numBands];
0962:
0963: // For the trivial case of first central moment return zeros.
0964: if (moment == 1 && isCentral) {
0965: for (int band = 0; band < numBands; band++) {
0966: moments[band] = 0.0;
0967: }
0968: } else {
0969: // Get the total counts for all bands.
0970: getTotals();
0971:
0972: for (int band = 0; band < numBands; band++) {
0973: // Cache some band-dependent quantities.
0974: int[] counts = getBins(band);
0975: int nBins = numBins[band];
0976: double level = getLowValue(band);
0977: double bw = binWidth[band];
0978: double total = totals[band];
0979:
0980: // Clear the moment value for this band.
0981: double mmt = 0.0;
0982:
0983: if (isCentral) {
0984: // Cache the mean value for this band.
0985: double mu = mean[band];
0986:
0987: if (isAbsolute && moment % 2 == 0) {
0988: // Even moment so absolute value has no effect.
0989: for (int b = 0; b < nBins; b++) {
0990: mmt += Math.pow(level - mu, moment)
0991: * counts[b] / total;
0992: level += bw;
0993: }
0994: } else {
0995: // Odd moment so need to use absolute value.
0996: for (int b = 0; b < nBins; b++) {
0997: mmt += Math.abs(Math
0998: .pow(level - mu, moment))
0999: * counts[b] / total;
1000: level += bw;
1001: }
1002: }
1003: } else if (isAbsolute && moment % 2 != 0) {
1004: // Odd moment so need to use absolute value.
1005: for (int b = 0; b < nBins; b++) {
1006: mmt += Math.abs(Math.pow(level, moment))
1007: * counts[b] / total;
1008: level += bw;
1009: }
1010: } else {
1011: // Even or non-absolute non-central moment.
1012: for (int b = 0; b < nBins; b++) {
1013: mmt += Math.pow(level, moment) * counts[b]
1014: / total;
1015: level += bw;
1016: }
1017: }
1018:
1019: // Save the result.
1020: moments[band] = mmt;
1021: }
1022: }
1023:
1024: return moments;
1025: }
1026:
1027: /**
1028: * Returns the standard deviation for all bands of the histogram.
1029: * This is a convenience method as the returned values
1030: * could be calculated using the first and second moments which
1031: * are available via the moment generation function.
1032: *
1033: * @return The standard deviation values for all bands.
1034: *
1035: * @since JAI 1.1
1036: */
1037: public double[] getStandardDeviation() {
1038: getMean();
1039:
1040: double[] variance = getMoment(2, false, false);
1041:
1042: double[] stdev = new double[numBands];
1043:
1044: for (int i = 0; i < variance.length; i++) {
1045: stdev[i] = Math.sqrt(variance[i] - mean[i] * mean[i]);
1046: }
1047:
1048: return stdev;
1049: }
1050:
1051: /**
1052: * Returns the entropy of the histogram.
1053: *
1054: * <p> The histogram entropy is defined to be the negation of the sum
1055: * of the products of the probability associated with each bin with
1056: * the base-2 log of the probability.
1057: *
1058: * @return The entropy of the histogram.
1059: *
1060: * @since JAI 1.1
1061: */
1062: public double[] getEntropy() {
1063: // Get the total counts for all bands.
1064: getTotals();
1065:
1066: double log2 = Math.log(2.0);
1067:
1068: double[] entropy = new double[numBands];
1069:
1070: for (int band = 0; band < numBands; band++) {
1071: int[] counts = getBins(band);
1072: int nBins = numBins[band];
1073: double total = totals[band];
1074:
1075: double H = 0.0;
1076:
1077: for (int b = 0; b < nBins; b++) {
1078: double p = counts[b] / total;
1079: if (p != 0.0) {
1080: H -= p * (Math.log(p) / log2);
1081: }
1082: }
1083:
1084: entropy[band] = H;
1085: }
1086:
1087: return entropy;
1088: }
1089:
1090: /**
1091: * Computes a smoothed version of the histogram.
1092: *
1093: * <p> Each band of the histogram is smoothed by averaging over a
1094: * moving window of a size specified by the method parameter: if
1095: * the value of the parameter is <i>k</i> then the width of the window
1096: * is <i>2*k + 1</i>. If the window runs off the end of the histogram
1097: * only those values which intersect the histogram are taken into
1098: * consideration. The smoothing may optionally be weighted to favor
1099: * the central value using a "triangular" weighting. For example,
1100: * for a value of <i>k</i> equal to 2 the central bin would have weight
1101: * 1/3, the adjacent bins 2/9, and the next adjacent bins 1/9.
1102: *
1103: * @param isWeighted Whether bins will be weighted using a triangular
1104: * weighting scheme favoring bins near the central bin.
1105: * @param k The smoothing parameter which must be non-negative or an
1106: * <code>IllegalArgumentException</code> will be thrown. If zero, the
1107: * histogram object will be returned with no smoothing applied.
1108: * @return A smoothed version of the histogram.
1109: *
1110: * @since JAI 1.1
1111: */
1112: public Histogram getSmoothed(boolean isWeighted, int k) {
1113: if (k < 0) {
1114: throw new IllegalArgumentException(JaiI18N
1115: .getString("Histogram7"));
1116: } else if (k == 0) {
1117: return this ;
1118: }
1119:
1120: // Create a new, identical but empty Histogram.
1121: Histogram smoothedHistogram = new Histogram(getNumBins(),
1122: getLowValue(), getHighValue());
1123:
1124: // Get a reference to the bins of the new Histogram.
1125: int[][] smoothedBins = smoothedHistogram.getBins();
1126:
1127: // Get the total counts for all bands.
1128: getTotals();
1129:
1130: // Initialize the smoothing weights if needed.
1131: double[] weights = null;
1132: if (isWeighted) {
1133: int numWeights = 2 * k + 1;
1134: double denom = numWeights * numWeights;
1135: weights = new double[numWeights];
1136: for (int i = 0; i <= k; i++) {
1137: weights[i] = (i + 1) / denom;
1138: }
1139: for (int i = k + 1; i < numWeights; i++) {
1140: weights[i] = weights[numWeights - 1 - i];
1141: }
1142: }
1143:
1144: for (int band = 0; band < numBands; band++) {
1145: // Cache bin-dependent values and references.
1146: int[] counts = getBins(band);
1147: int[] smoothedCounts = smoothedBins[band];
1148: int nBins = smoothedHistogram.getNumBins(band);
1149:
1150: // Clear the band total count for the smoothed histogram.
1151: int sum = 0;
1152:
1153: if (isWeighted) {
1154: for (int b = 0; b < nBins; b++) {
1155: // Determine clipped range.
1156: int min = Math.max(b - k, 0);
1157: int max = Math.min(b + k, nBins);
1158:
1159: // Calculate the offset into the weight array.
1160: int offset = k > b ? k - b : 0;
1161:
1162: // Accumulate the total for the range.
1163: double acc = 0;
1164: double weightTotal = 0;
1165: for (int i = min; i < max; i++) {
1166: double w = weights[offset++];
1167: acc += counts[i] * w;
1168: weightTotal += w;
1169: }
1170:
1171: // Round the accumulated value.
1172: smoothedCounts[b] = (int) (acc / weightTotal + 0.5);
1173:
1174: // Accumulate total for band.
1175: sum += smoothedCounts[b];
1176: }
1177: } else {
1178: for (int b = 0; b < nBins; b++) {
1179: // Determine clipped range.
1180: int min = Math.max(b - k, 0);
1181: int max = Math.min(b + k, nBins);
1182:
1183: // Accumulate the total for the range.
1184: int acc = 0;
1185: for (int i = min; i < max; i++) {
1186: acc += counts[i];
1187: }
1188:
1189: // Calculate the average for the range.
1190: smoothedCounts[b] = (int) (acc
1191: / (double) (max - min + 1) + 0.5);
1192:
1193: // Accumulate total for band.
1194: sum += smoothedCounts[b];
1195: }
1196: }
1197:
1198: // Rescale the counts such that the band total is approximately
1199: // the same as for the same band of the original histogram.
1200: double factor = (double) totals[band] / (double) sum;
1201: for (int b = 0; b < nBins; b++) {
1202: smoothedCounts[b] = (int) (smoothedCounts[b] * factor + 0.5);
1203: }
1204: }
1205:
1206: return smoothedHistogram;
1207: }
1208:
1209: /**
1210: * Computes a Gaussian smoothed version of the histogram.
1211: *
1212: * <p> Each band of the histogram is smoothed by discrete convolution
1213: * with a kernel approximating a Gaussian impulse response with the
1214: * specified standard deviation.
1215: *
1216: * @param standardDeviation The standard deviation of the Gaussian
1217: * smoothing kernel which must be non-negative or an
1218: * <code>IllegalArgumentException</code> will be thrown. If zero, the
1219: * histogram object will be returned with no smoothing applied.
1220: * @return A Gaussian smoothed version of the histogram.
1221: *
1222: * @since JAI 1.1
1223: */
1224: public Histogram getGaussianSmoothed(double standardDeviation) {
1225: if (standardDeviation < 0.0) {
1226: throw new IllegalArgumentException(JaiI18N
1227: .getString("Histogram8"));
1228: } else if (standardDeviation == 0.0) {
1229: return this ;
1230: }
1231:
1232: // Create a new, identical but empty Histogram.
1233: Histogram smoothedHistogram = new Histogram(getNumBins(),
1234: getLowValue(), getHighValue());
1235:
1236: // Get a reference to the bins of the new Histogram.
1237: int[][] smoothedBins = smoothedHistogram.getBins();
1238:
1239: // Get the total counts for all bands.
1240: getTotals();
1241:
1242: // Determine the number of weights (must be odd).
1243: int numWeights = (int) (2 * 2.58 * standardDeviation + 0.5);
1244: if (numWeights % 2 == 0) {
1245: numWeights++;
1246: }
1247:
1248: // Initialize the smoothing weights.
1249: double[] weights = new double[numWeights];
1250: int m = numWeights / 2;
1251: double var = standardDeviation * standardDeviation;
1252: double gain = 1.0 / Math.sqrt(2.0 * Math.PI * var);
1253: double exp = -1.0 / (2.0 * var);
1254: for (int i = m; i < numWeights; i++) {
1255: double del = i - m;
1256: weights[i] = weights[numWeights - 1 - i] = gain
1257: * Math.exp(exp * del * del);
1258: }
1259:
1260: for (int band = 0; band < numBands; band++) {
1261: // Cache bin-dependent values and references.
1262: int[] counts = getBins(band);
1263: int[] smoothedCounts = smoothedBins[band];
1264: int nBins = smoothedHistogram.getNumBins(band);
1265:
1266: // Clear the band total count for the smoothed histogram.
1267: int sum = 0;
1268:
1269: for (int b = 0; b < nBins; b++) {
1270: // Determine clipped range.
1271: int min = Math.max(b - m, 0);
1272: int max = Math.min(b + m, nBins);
1273:
1274: // Calculate the offset into the weight array.
1275: int offset = m > b ? m - b : 0;
1276:
1277: // Accumulate the total for the range.
1278: double acc = 0;
1279: double weightTotal = 0;
1280: for (int i = min; i < max; i++) {
1281: double w = weights[offset++];
1282: acc += counts[i] * w;
1283: weightTotal += w;
1284: }
1285:
1286: // Round the accumulated value.
1287: smoothedCounts[b] = (int) (acc / weightTotal + 0.5);
1288:
1289: // Accumulate total for band.
1290: sum += smoothedCounts[b];
1291: }
1292:
1293: // Rescale the counts such that the band total is approximately
1294: // the same as for the same band of the original histogram.
1295: double factor = (double) totals[band] / (double) sum;
1296: for (int b = 0; b < nBins; b++) {
1297: smoothedCounts[b] = (int) (smoothedCounts[b] * factor + 0.5);
1298: }
1299: }
1300:
1301: return smoothedHistogram;
1302: }
1303:
1304: /**
1305: * Calculates the <i>p-tile</i> threshold.
1306: *
1307: * <p> Computes thresholds such that a specified proportion of the sample
1308: * values in each band are below the threshold.
1309: *
1310: * @param p The proportion of samples in each band which should be below
1311: * the threshold in the band. If <code>p</code> is not in the range
1312: * (0.0, 1.0) an <code>IllegalArgumentException</code> will be thrown.
1313: * @return The requested <i>p-tile</i> thresholds.
1314: *
1315: * @since JAI 1.1
1316: */
1317: public double[] getPTileThreshold(double p) {
1318: if (p <= 0.0 || p >= 1.0) {
1319: throw new IllegalArgumentException(JaiI18N
1320: .getString("Histogram9"));
1321: }
1322:
1323: double[] thresholds = new double[numBands];
1324: getTotals();
1325:
1326: for (int band = 0; band < numBands; band++) {
1327: // Cache some band-dependent values.
1328: int nBins = numBins[band];
1329: int[] counts = getBins(band);
1330:
1331: // Calculate the total count for this band.
1332: int totalCount = totals[band];
1333:
1334: // Determine the number of binWidths to add to the lowValue
1335: // to get the desired threshold.
1336: int numBinWidths = 0;
1337: int count = counts[0];
1338: int idx = 0;
1339: while ((double) count / (double) totalCount < p) {
1340: numBinWidths++;
1341: count += counts[++idx];
1342: }
1343:
1344: // Calculate the threshold.
1345: thresholds[band] = getLowValue(band) + numBinWidths
1346: * binWidth[band];
1347: }
1348:
1349: return thresholds;
1350: }
1351:
1352: /**
1353: * Calculates the threshold using the mode method.
1354: *
1355: * <p> The threshold is defined to be the minimum between two peaks.
1356: * The first peak is the highest peak in the histogram. The second
1357: * peak is the highest peak in the histogram weighted by a specified
1358: * power of the distance from the first peak.
1359: *
1360: * @param power The exponent of the distance weighting from the
1361: * first peak.
1362: * @return The requested thresholds.
1363: *
1364: * @since JAI 1.1
1365: */
1366: public double[] getModeThreshold(double power) {
1367: double[] thresholds = new double[numBands];
1368: getTotals();
1369:
1370: for (int band = 0; band < numBands; band++) {
1371: // Cache some band-dependent values.
1372: int nBins = numBins[band];
1373: int[] counts = getBins(band);
1374:
1375: // Find the primary mode (highest peak).
1376: int mode1 = 0;
1377: int mode1Count = counts[0];
1378: for (int b = 1; b < nBins; b++) {
1379: if (counts[b] > mode1Count) {
1380: mode1 = b;
1381: mode1Count = counts[b];
1382: }
1383: }
1384:
1385: // Find the secondary mode (highest weighted peak).
1386: int mode2 = -1;
1387: double mode2count = 0.0;
1388: for (int b = 0; b < nBins; b++) {
1389: double d = counts[b]
1390: * Math.pow(Math.abs(b - mode1), power);
1391: if (d > mode2count) {
1392: mode2 = b;
1393: mode2count = d;
1394: }
1395: }
1396:
1397: // Find the minimum value between the two peaks.
1398: int min = mode1;
1399: int minCount = counts[mode1];
1400: for (int b = mode1 + 1; b <= mode2; b++) {
1401: if (counts[b] < minCount) {
1402: min = b;
1403: minCount = counts[b];
1404: }
1405: }
1406:
1407: thresholds[band] = (int) ((mode1 + mode2) / 2.0 + 0.5);
1408: }
1409:
1410: return thresholds;
1411: }
1412:
1413: /**
1414: * Calculates the threshold using iterative bisection.
1415: *
1416: * <p> For each band an initial threshold is defined to be the midpoint
1417: * of the range of data represented by the histogram. The mean value is
1418: * calculated for each sub-histogram and a new threshold is defined as the
1419: * arithmetic mean of the two sub-histogram means. This process is
1420: * repeated until the threshold value no longer changes.
1421: *
1422: * @return The requested thresholds.
1423: *
1424: * @since JAI 1.1
1425: */
1426: public double[] getIterativeThreshold() {
1427: double[] thresholds = new double[numBands];
1428: getTotals();
1429:
1430: for (int band = 0; band < numBands; band++) {
1431: // Cache some band-dependent values.
1432: int nBins = numBins[band];
1433: int[] counts = getBins(band);
1434: double bw = binWidth[band];
1435:
1436: // Set intial threshold to midpoint of data range for this band.
1437: double threshold = 0.5 * (getLowValue(band) + getHighValue(band));
1438: double mid1 = 0.5 * (getLowValue(band) + threshold);
1439: double mid2 = 0.5 * (threshold + getHighValue(band));
1440:
1441: // Iterate only if total is nonzero.
1442: if (totals[band] != 0) {
1443:
1444: // Loop until threshold estimate no longer changes.
1445: int countDown = 1000;
1446: do {
1447: // Save band threshold for this iteration.
1448: thresholds[band] = threshold;
1449:
1450: // Cache the total count for this band.
1451: double total = totals[band];
1452:
1453: // Initialize the level corresponding to a bin.
1454: double level = getLowValue(band);
1455:
1456: // Clear mean values for sub-ranges.
1457: double mean1 = 0.0;
1458: double mean2 = 0.0;
1459:
1460: // Clear sub-range 1 count.
1461: int count1 = 0;
1462:
1463: // Calculate the mean values for the two sub-ranges.
1464: for (int b = 0; b < nBins; b++) {
1465: // Update the mean value for the appropriate sub-range.
1466: if (level <= threshold) {
1467: int c = counts[b];
1468: mean1 += c * level;
1469: count1 += c;
1470: } else {
1471: mean2 += counts[b] * level;
1472: }
1473:
1474: // Augment the level for the current bin by the bin width.
1475: level += bw;
1476: }
1477:
1478: // Rescale values using sub-range totals.
1479: if (count1 != 0) {
1480: mean1 /= count1;
1481: } else {
1482: mean1 = mid1;
1483: }
1484: if (total != count1) {
1485: mean2 /= (total - count1);
1486: } else {
1487: mean2 = mid2;
1488: }
1489:
1490: // Update the threshold estimate.
1491: threshold = 0.5 * (mean1 + mean2);
1492: } while (Math.abs(threshold - thresholds[band]) > 1e-6
1493: && --countDown > 0);
1494: } else {
1495: thresholds[band] = threshold;
1496: }
1497: }
1498:
1499: return thresholds;
1500: }
1501:
1502: /**
1503: * Calculates the threshold which maximizes the ratio of the between-class
1504: * variance to the within-class variance for each band.
1505: *
1506: * @return The requested thresholds.
1507: *
1508: * @since JAI 1.1
1509: */
1510: public double[] getMaxVarianceThreshold() {
1511: double[] thresholds = new double[numBands];
1512: getTotals();
1513: getMean();
1514: double[] variance = getMoment(2, false, false);
1515:
1516: for (int band = 0; band < numBands; band++) {
1517: // Cache some band-dependent values.
1518: int nBins = numBins[band];
1519: int[] counts = getBins(band);
1520: double total = totals[band];
1521: double mBand = mean[band];
1522: double bw = binWidth[band];
1523:
1524: double prob0 = 0.0;
1525: double mean0 = 0.0;
1526: double lv = getLowValue(band);
1527: double level = lv;
1528: double maxRatio = -Double.MAX_VALUE;
1529: int maxIndex = 0;
1530: int runLength = 0;
1531:
1532: for (int t = 0; t < nBins; t++, level += bw) {
1533: double p = counts[t] / total;
1534: prob0 += p;
1535: if (prob0 == 0.0) {
1536: continue;
1537: }
1538:
1539: double m0 = (mean0 += p * level) / prob0;
1540:
1541: double prob1 = 1.0 - prob0;
1542:
1543: if (prob1 == 0.0) {
1544: continue;
1545: }
1546:
1547: double m1 = (mBand - mean0) / prob1;
1548:
1549: double var0 = 0.0;
1550: double g = lv;
1551: for (int b = 0; b <= t; b++, g += bw) {
1552: double del = g - m0;
1553: var0 += del * del * counts[b];
1554: }
1555: var0 /= total;
1556:
1557: double var1 = 0.0;
1558: for (int b = t + 1; b < nBins; b++, g += bw) {
1559: double del = g - m1;
1560: var1 += del * del * counts[b];
1561: }
1562: var1 /= total;
1563:
1564: if (var0 == 0.0 && var1 == 0.0 && m1 != 0.0) {
1565: maxIndex = (int) (((m0 + m1) / 2.0 - getLowValue(band))
1566: / bw + 0.5);
1567: runLength = 0;
1568: break;
1569: }
1570:
1571: if (var0 / prob0 < 0.5 || var1 / prob1 < 0.5) {
1572: continue;
1573: }
1574:
1575: double mdel = m0 - m1;
1576: double ratio = prob0 * prob1 * mdel * mdel
1577: / (var0 + var1);
1578:
1579: if (ratio > maxRatio) {
1580: maxRatio = ratio;
1581: maxIndex = t;
1582: runLength = 0;
1583: } else if (ratio == maxRatio) {
1584: runLength++;
1585: }
1586: }
1587:
1588: thresholds[band] = getLowValue(band)
1589: + (maxIndex + runLength / 2.0 + 0.5) * bw;
1590: }
1591:
1592: return thresholds;
1593: }
1594:
1595: /**
1596: * Calculates the threshold which maximizes the entropy.
1597: *
1598: * <p> The <i>entropy</i> of a range of gray levels is defined to
1599: * be the negation of the sum of products of the probability and
1600: * the logarithm thereof over all gray levels in the range. The
1601: * maximum entropy threshold is defined to be that value which maximizes
1602: * the sum of the entropy of the two ranges which are above and below the
1603: * threshold, respectively. This computation is effected for each band.
1604: *
1605: * @return The requested thresholds.
1606: *
1607: * @since JAI 1.1
1608: */
1609: public double[] getMaxEntropyThreshold() {
1610: double[] thresholds = new double[numBands];
1611: getTotals();
1612:
1613: double[] entropy = getEntropy();
1614:
1615: double log2 = Math.log(2.0);
1616:
1617: for (int band = 0; band < numBands; band++) {
1618: // Cache some band-dependent values.
1619: int nBins = numBins[band];
1620: int[] counts = getBins(band);
1621: double total = totals[band];
1622: double H = entropy[band];
1623:
1624: double P1 = 0.0;
1625: double H1 = 0.0;
1626:
1627: double maxCriterion = -Double.MAX_VALUE;
1628: int maxIndex = 0;
1629: int runLength = 0;
1630:
1631: for (int t = 0; t < nBins; t++) {
1632: double p = counts[t] / total;
1633:
1634: if (p == 0.0) {
1635: continue;
1636: }
1637:
1638: P1 += p;
1639:
1640: H1 -= p * Math.log(p) / log2;
1641:
1642: double max1 = 0.0;
1643: for (int b = 0; b <= t; b++) {
1644: if (counts[b] > max1) {
1645: max1 = counts[b];
1646: }
1647: }
1648:
1649: if (max1 == 0.0) {
1650: continue;
1651: }
1652:
1653: double max2 = 0.0;
1654: for (int b = t + 1; b < nBins; b++) {
1655: if (counts[b] > max2) {
1656: max2 = counts[b];
1657: }
1658: }
1659:
1660: if (max2 == 0.0) {
1661: continue;
1662: }
1663:
1664: double ratio = H1 / H;
1665:
1666: double criterion = ratio * Math.log(P1)
1667: / Math.log(max1 / total) + (1.0 - ratio)
1668: * Math.log(1.0 - P1) / Math.log(max2 / total);
1669:
1670: if (criterion > maxCriterion) {
1671: maxCriterion = criterion;
1672: maxIndex = t;
1673: runLength = 0;
1674: } else if (criterion == maxCriterion) {
1675: runLength++;
1676: }
1677: }
1678:
1679: thresholds[band] = getLowValue(band)
1680: + (maxIndex + runLength / 2.0 + 0.5)
1681: * binWidth[band];
1682: }
1683:
1684: return thresholds;
1685: }
1686:
1687: /**
1688: * Calculates the threshold which minimizes the probability of error.
1689: *
1690: * <p> For each band the histogram is modeled as the sum of two Gaussian
1691: * distributions and the threshold which minimizes the misclassification
1692: * error is computed. If the underlying histogram is unimodal the mean
1693: * value of each band will be returned as the threshold.
1694: * The bimodality of the histogram for that band will be identically zero.
1695: *
1696: * @return The requested thresholds.
1697: *
1698: * @since JAI 1.1
1699: */
1700: public double[] getMinErrorThreshold() {
1701: double[] thresholds = new double[numBands];
1702: getTotals();
1703: getMean();
1704:
1705: for (int band = 0; band < numBands; band++) {
1706: // Cache some band-dependent values.
1707: int nBins = numBins[band];
1708: int[] counts = getBins(band);
1709: double total = totals[band];
1710: double lv = getLowValue(band);
1711: double bw = binWidth[band];
1712:
1713: int total1 = 0;
1714: int total2 = totals[band];
1715: double sum1 = 0.0;
1716: double sum2 = mean[band] * total;
1717:
1718: double level = lv;
1719:
1720: double minCriterion = Double.MAX_VALUE;
1721: int minIndex = 0;
1722: int runLength = 0;
1723:
1724: double J0 = Double.MAX_VALUE;
1725: double J1 = Double.MAX_VALUE;
1726: double J2 = Double.MAX_VALUE;
1727: int Jcount = 0;
1728:
1729: for (int t = 0; t < nBins; t++, level += bw) {
1730: int c = counts[t];
1731:
1732: total1 += c;
1733: total2 -= c;
1734:
1735: double incr = level * c;
1736: sum1 += incr;
1737: sum2 -= incr;
1738:
1739: if (total1 == 0 || sum1 == 0.0) {
1740: continue;
1741: } else if (total2 == 0 || sum2 == 0.0) {
1742: break;
1743: }
1744:
1745: double m1 = sum1 / total1;
1746: double m2 = sum2 / total2;
1747:
1748: double s1 = 0.0;
1749: double g = lv;
1750: for (int b = 0; b <= t; b++, g += bw) {
1751: double v = g - m1;
1752: s1 += counts[b] * v * v;
1753: }
1754: s1 /= total1;
1755:
1756: if (s1 < 0.5) {
1757: continue;
1758: }
1759:
1760: double s2 = 0.0;
1761: //g = lv;
1762: for (int b = t + 1; b < nBins; b++, g += bw) {
1763: double v = g - m2;
1764: s2 += counts[b] * v * v;
1765: }
1766: s2 /= total2;
1767:
1768: if (s2 < 0.5) {
1769: continue;
1770: }
1771:
1772: double P1 = total1 / total;
1773: double P2 = total2 / total;
1774: double J = 1.0 + P1 * Math.log(s1) + P2 * Math.log(s2)
1775: - 2.0 * (P1 * Math.log(P1) + P2 * Math.log(P2));
1776:
1777: Jcount++;
1778:
1779: J0 = J1;
1780: J1 = J2;
1781: J2 = J;
1782:
1783: if (Jcount >= 3) {
1784: if (J1 <= J0 && J1 <= J2) {
1785: if (J1 < minCriterion) {
1786: minCriterion = J1;
1787: minIndex = t - 1;
1788: runLength = 0;
1789: } else if (J1 == minCriterion) {
1790: runLength++;
1791: }
1792: }
1793: }
1794: }
1795:
1796: thresholds[band] = minIndex == 0 ? mean[band]
1797: : getLowValue(band)
1798: + (minIndex + runLength / 2.0 + 0.5) * bw;
1799: }
1800:
1801: return thresholds;
1802: }
1803:
1804: /**
1805: * Calculates the threshold which minimizes the <i>fuzziness</i>.
1806: *
1807: * @since JAI 1.1
1808: */
1809: public double[] getMinFuzzinessThreshold() {
1810: double[] thresholds = new double[numBands];
1811: getTotals();
1812: getMean();
1813:
1814: for (int band = 0; band < numBands; band++) {
1815: // Cache some band-dependent values.
1816: int nBins = numBins[band];
1817: int[] counts = getBins(band);
1818: double total = totals[band];
1819:
1820: double bw = binWidth[band];
1821:
1822: int total1 = 0;
1823: int total2 = totals[band];
1824: double sum1 = 0.0;
1825: double sum2 = mean[band] * total;
1826:
1827: double lv = getLowValue(band);
1828: double level = lv;
1829: double C = getHighValue(band) - lv;
1830:
1831: double minCriterion = Double.MAX_VALUE;
1832: int minIndex = 0;
1833: int runLength = 0;
1834:
1835: for (int t = 0; t < nBins; t++, level += bw) {
1836: int c = counts[t];
1837:
1838: total1 += c;
1839: total2 -= c;
1840:
1841: double incr = level * c;
1842: sum1 += incr;
1843: sum2 -= incr;
1844:
1845: if (total1 == 0 || total2 == 0) {
1846: continue;
1847: }
1848:
1849: double m1 = sum1 / total1;
1850: double m2 = sum2 / total2;
1851:
1852: double g = lv;
1853: double E = 0.0;
1854: for (int b = 0; b < nBins; b++, g += bw) {
1855: double u = b <= t ? 1.0 / (1.0 + Math.abs(g - m1)
1856: / C) : 1.0 / (1.0 + Math.abs(g - m2) / C);
1857: double v = 1 - u;
1858: E += (-u * Math.log(u) - v * Math.log(v))
1859: * (counts[b] / total);
1860: }
1861:
1862: if (E < minCriterion) {
1863: minCriterion = E;
1864: minIndex = t;
1865: runLength = 0;
1866: } else if (E == minCriterion) {
1867: runLength++;
1868: }
1869: }
1870:
1871: thresholds[band] = lv + (minIndex + runLength / 2.0 + 0.5)
1872: * bw;
1873: }
1874:
1875: return thresholds;
1876: }
1877: }
|