package Algorithms;

import Utils.*;

public class Lacunarity2
{
    public static class Statistics
    {
        public IntVector boxTable = new IntVector();
        public IntVector el3 = new IntVector();
        public IntVector fl3 = new IntVector();
        public DoubleVector eLambda3 = new DoubleVector();
        public DoubleVector fLambda3 = new DoubleVector();
        public DoubleVector logBoxSizes = new DoubleVector();
        public DoubleVector logElambda3p1 = new DoubleVector();
        public DoubleVector logFlambda3p1 = new DoubleVector();

        public double elCurve;
        public double flCurve;
        public double elMedial;
        public double flMedial;
        public double elMean;
        public double flMean;
    }

    public static void computeLacunarity(Statistics stats, byte[] image, int width, int height, int numberOfBins, int minSize, int inc)
    {
        int minX = width;
        int minY = height;
        int maxX = 0;
        int maxY = 0;

        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                int isWhite = image[x + width * y] >> 31; // no AND with 0xff
                minX = Math.min(minX, (x & isWhite) | (width & ~isWhite));
                minY = Math.min(minY, (y & isWhite) | (height & ~isWhite));
                maxX = Math.max(maxX, x & isWhite);
                maxY = Math.max(maxY, y & isWhite);
            }
        }

        final int selWidth = maxX - minX + 1;
        final int selHeight = maxY - minY + 1;
        final int smallDimension = Math.min(selWidth, selHeight);
        final int factor = (smallDimension - minSize) / numberOfBins;

        int minMedialDiff = -1;
        int medialBox = 0;

        double eLacTotal = 0.0;
        double fLacTotal = 0.0;

        stats.boxTable.size = 0;
        stats.el3.size = 0;
        stats.fl3.size = 0;
        stats.eLambda3.size = 0;
        stats.fLambda3.size = 0;
        stats.logBoxSizes.size = 0;
        stats.logElambda3p1.size = 0;
        stats.logFlambda3p1.size = 0;

        for (int i = 0; i < numberOfBins; i++) {
            final int boxSize = (i < numberOfBins - 1) ? (minSize + i*factor) : smallDimension;

            int diff = Math.abs(boxSize - smallDimension / 2);
            if (minMedialDiff < 0 || diff < minMedialDiff) {
                minMedialDiff = diff;
                medialBox = i;
            }

            int numXBox = (int)Math.ceil((double)selWidth / (double)inc - (double)boxSize / (double)inc + 1.0);
            int numYBox = (int)Math.ceil((double)selHeight / (double)inc - (double)boxSize / (double)inc + 1.0);
            stats.boxTable.resize(numXBox * numYBox * 2);

            // 00 case
            int topLeftSlice  = countWhitePixels(image, width, height, minX, minY, inc, inc);
            int topMainSlice  = countWhitePixels(image, width, height, minX + inc, minY, boxSize - inc, inc);
            int mainLeftSlice = countWhitePixels(image, width, height, minX, minY + inc, inc, boxSize - inc);
            int mainMainSlice = countWhitePixels(image, width, height, minX + inc, minY + inc, boxSize - inc, boxSize - inc);

            stats.boxTable.buf[0] = topLeftSlice + topMainSlice + mainLeftSlice + mainMainSlice;
            stats.boxTable.buf[1] = topLeftSlice + topMainSlice;

            for (int col = 1; col < numXBox; col++) {
                // 0j case
                int x = minX + col * inc;
                int topRightSlice  = countWhitePixels(image, width, height, x + boxSize - inc, minY, inc, inc);
                int mainRightSlice = countWhitePixels(image, width, height, x + boxSize - inc, minY + inc, inc, boxSize - inc);

                stats.boxTable.buf[2*col]     = stats.boxTable.buf[2*(col-1)] + (topRightSlice + mainRightSlice) - (topLeftSlice + mainLeftSlice);
                stats.boxTable.buf[2*col + 1] = stats.boxTable.buf[2*(col-1) + 1] + topRightSlice - topLeftSlice;

                if (col < numXBox - 1) {
                    topLeftSlice  = countWhitePixels(image, width, height, x, minY, inc, inc);
                    mainLeftSlice = countWhitePixels(image, width, height, x, minY + inc, inc, boxSize - inc);
                }
            }

            for (int row = 1; row < numYBox; row++) {
                // i0 case
                int y = minY + row * inc;
                int bottomLeftSlice = countWhitePixels(image, width, height, minX, y + boxSize - inc, inc, inc);
                int bottomMainSlice = countWhitePixels(image, width, height, minX + inc, y + boxSize - inc, boxSize - inc, inc);

                topLeftSlice = countWhitePixels(image, width, height, minX, y, inc, inc);
                topMainSlice = countWhitePixels(image, width, height, minX + inc, y, boxSize - inc, inc);

                int bottomSlice = bottomMainSlice + bottomLeftSlice;
                int topSlice = topLeftSlice + topMainSlice;

                stats.boxTable.buf[2*row*numXBox] =
                    stats.boxTable.buf[2*(row-1)*numXBox] +
                    bottomSlice -
                    stats.boxTable.buf[2*(row-1)*numXBox + 1];

                stats.boxTable.buf[2*row*numXBox + 1] = topSlice;

                for (int col = 1; col < numXBox; col++) {
                    // ij case
                    int x = minX + col * inc;
                    int bottomRightSlice = countWhitePixels(image, width, height, x + boxSize - inc, y + boxSize - inc, inc, inc);
                    int topRightSlice = countWhitePixels(image, width, height, x + boxSize - inc, y, inc, inc);

                    bottomSlice += bottomRightSlice - bottomLeftSlice;
                    topSlice += topRightSlice - topLeftSlice;

                    int pos = row * numXBox + col;

                    stats.boxTable.buf[2*pos] =
                        stats.boxTable.buf[2*(pos-numXBox)] +
                        bottomSlice -
                        stats.boxTable.buf[2*(pos-numXBox) + 1];

                    stats.boxTable.buf[2*pos + 1] = topSlice;

                    if (col < numXBox - 1) {
                        bottomLeftSlice = countWhitePixels(image, width, height, x, y + boxSize - inc, inc, inc);
                        topLeftSlice  = countWhitePixels(image, width, height, x, y, inc, inc);
                    }
                }
            }

            stats.el3.size = 0;
            stats.fl3.size = 0;

            int totalBoxes = numXBox * numYBox;
            for (int j = 0; j < totalBoxes; j++) {
                int count = stats.boxTable.buf[j*2];
                stats.el3.add(count);
                if (count > 0)
                    stats.fl3.add(count);
            }

            double eLac = computeCoefficientOfVariationSquared(stats.el3.buf, stats.el3.size);
            double fLac = computeCoefficientOfVariationSquared(stats.fl3.buf, stats.fl3.size);
            eLacTotal += eLac;
            fLacTotal += fLac;
            stats.eLambda3.add(eLac);
            stats.fLambda3.add(fLac);

            stats.logBoxSizes.add(Math.log((double)boxSize));
            stats.logElambda3p1.add(Math.log(eLac + 1.0));
            stats.logFlambda3p1.add(Math.log(fLac + 1.0));
        }

        stats.elMean = stats.eLambda3.size > 0 ? eLacTotal / (double)stats.eLambda3.size : 0.0;
        stats.flMean = stats.fLambda3.size > 0 ? fLacTotal / (double)stats.fLambda3.size : 0.0;
        stats.elMedial = stats.eLambda3.buf[medialBox];
        stats.flMedial = stats.fLambda3.buf[medialBox];
        stats.elCurve = findLinearRegressionFactor(stats.logBoxSizes.buf, stats.logElambda3p1.buf, numberOfBins, null);
        stats.flCurve = findLinearRegressionFactor(stats.logBoxSizes.buf, stats.logFlambda3p1.buf, numberOfBins, null);
    }

    public static int countWhitePixels(byte[] image, int width, int height, int xStart, int yStart, int dx, int dy)
    {
        int count = 0;
        int xLast = Math.min(xStart + dx, width) - 1;
        int yLast = Math.min(yStart + dy, height) - 1;

        for (int y = yStart; y <= yLast; y++) {
            for (int x = xStart; x <= xLast; x++) {
                int pixel = image[x + width * y]; // no AND with 0xff
                count -= pixel >> 31; // -1 if the pixel is white, 0 if not
            }
        }

        return count;
    }

    public static double computeCoefficientOfVariationSquared(int[] buf, int n)
    {
        if (n <= 0)
            return 0.0;

        double total = 0;
        for (int i = 0; i < n; i++)
            total += buf[i];

        double avg = total / (double)n;
        double variation = 0.0;
        for (int i = 0; i < n; i++) {
            double dAvg = (double)buf[i] - avg;
            variation += dAvg * dAvg;
        }

        double stdDev = Math.sqrt(variation / (double)n);
        double cv = stdDev / avg;
        return cv * cv;
    }

    // Modified snippet, originally from ij/measure/CurveFitter.java:
    //   Determine sum of squared residuals with linear regression.
    //   The sum of squared residuals is written to the array element with index 'numParams',
    //   the offset and factor params (if any) are written to their proper positions in the
    //   params array
    public static double findLinearRegressionFactor(double[] xData, double[] yData, int numPoints, double[] params)
    {
        double sumX=0, sumX2=0, sumXY=0; // sums for regression; here 'x' are function values
        double sumY=0, sumY2=0;          // only calculated for 'slope', otherwise we use the values calculated already

        for (int i=0; i<numPoints; i++) {
            double x = xData[i];
            double y = yData[i];
            sumX += x;
            sumX2 += x*x;
            sumXY += x*y;
            sumY2 += y*y;
            sumY += y;
        }

        final double sumWeights = (double)numPoints;

        // full linear regression or offset only. Slope is named 'factor' here
        double factor = (sumXY-sumX*sumY/sumWeights)/(sumX2-sumX*sumX/sumWeights);
        if (Double.isNaN(factor) || Double.isInfinite(factor))
            factor = 0; // all 'x' values are equal, any factor (slope) will fit

        if (params != null) {
            double offset = (sumY-factor*sumX)/sumWeights;

            double factorSqrSumX2 = factor*factor*sumX2;
            double offsetSqrSumWeights = sumWeights*offset*offset;
            double sumResidualsSqr = factorSqrSumX2 + offsetSqrSumWeights + sumY2 + 2*factor*offset*sumX - 2*factor*sumXY - 2*offset*sumY;

            // check for accuracy problem: large difference of small numbers?
            // Don't report unrealistic or even negative values, otherwise minimization could lead
            // into parameters where we have a numeric problem
            sumResidualsSqr = Math.max(sumResidualsSqr, 2e-15*(factorSqrSumX2 + offsetSqrSumWeights + sumY2));

            params[0] = offset;
            params[1] = factor;
            params[2] = sumResidualsSqr;
        }

        return factor;
    }
}