Examples with FloatArray2D mpi.fruitfly.math.datastructures.FloatArray2D used on opensource projects

Example 1 with FloatArray2D

use of mpi.fruitfly.math.datastructures.FloatArray2D in project TrakEM2 by trakem2.

the class StitchingTEM method correlate.

/**
 * @param scale For optimizing the speed of phase- and cross-correlation.
 * @param percent_overlap The minimum chunk of adjacent images to compare with, will automatically and gradually increase to 100% if no good matches are found.
 * @return a double[4] array containing:<ul>
 * <li>x2: relative X position of the second Patch</li>
 * <li>y2: relative Y position of the second Patch</li>
 * <li>flag: ERROR or SUCCESS</li>
 * <li>R: cross-correlation coefficient</li>
 * </ul>
 */
public static double[] correlate(final Patch base, final Patch moving, final float percent_overlap, final double scale, final int direction, final double default_dx, final double default_dy, final double min_R) {
    // PhaseCorrelation2D pc = null;
    final double R = -2;
    // final int limit = 5; // number of peaks to check in the PhaseCorrelation results
    // final float min_R = 0.40f; // minimum R for phase-correlation to be considered good
    // half this min_R will be considered good for cross-correlation
    // Iterate until PhaseCorrelation correlation coefficient R is over 0.5, or there's no more
    // image overlap to feed
    // Utils.log2("min_R: " + min_R);
    ImageProcessor ip1, ip2;
    final Rectangle b1 = base.getBoundingBox(null);
    final Rectangle b2 = moving.getBoundingBox(null);
    final int w1 = b1.width, h1 = b1.height, w2 = b2.width, h2 = b2.height;
    Roi roi1 = null, roi2 = null;
    float overlap = percent_overlap;
    double dx = default_dx, dy = default_dy;
    do {
        // create rois for the stripes
        switch(direction) {
            case TOP_BOTTOM:
                // bottom
                roi1 = new Roi(0, h1 - (int) (h1 * overlap), w1, (int) (h1 * overlap));
                // top
                roi2 = new Roi(0, 0, w2, (int) (h2 * overlap));
                break;
            case LEFT_RIGHT:
                // right
                roi1 = new Roi(w1 - (int) (w1 * overlap), 0, (int) (w1 * overlap), h1);
                // left
                roi2 = new Roi(0, 0, (int) (w2 * overlap), h2);
                break;
        }
        // Utils.log2("roi1: " + roi1);
        // Utils.log2("roi2: " + roi2);
        // will apply the transform if necessary
        ip1 = makeStripe(base, roi1, scale);
        ip2 = makeStripe(moving, roi2, scale);
        // new ImagePlus("roi1", ip1).show();
        // new ImagePlus("roi2", ip2).show();
        ip1.setPixels(ImageFilter.computeGaussianFastMirror(new FloatArray2D((float[]) ip1.getPixels(), ip1.getWidth(), ip1.getHeight()), 1.0).data);
        ip2.setPixels(ImageFilter.computeGaussianFastMirror(new FloatArray2D((float[]) ip2.getPixels(), ip2.getWidth(), ip2.getHeight()), 1.0).data);
        // 
        final ImagePlus imp1 = new ImagePlus("", ip1);
        final ImagePlus imp2 = new ImagePlus("", ip2);
        final PhaseCorrelationCalculator t = new PhaseCorrelationCalculator(imp1, imp2);
        final PhaseCorrelationPeak peak = t.getPeak();
        final double resultR = peak.getCrossCorrelationPeak();
        final int[] peackPostion = peak.getPosition();
        final java.awt.Point shift = new java.awt.Point(peackPostion[0], peackPostion[1]);
        // Utils.log2("overlap: " + overlap + " R: " + resultR + " shift: " + shift + " dx,dy: " + dx + ", " + dy);
        if (resultR >= min_R) {
            // success
            final int success = SUCCESS;
            switch(direction) {
                case TOP_BOTTOM:
                    // boundary checks:
                    // if (shift.y/scale > default_dy) success = ERROR;
                    dx = shift.x / scale;
                    dy = roi1.getBounds().y + shift.y / scale;
                    break;
                case LEFT_RIGHT:
                    // boundary checks:
                    // if (shift.x/scale > default_dx) success = ERROR;
                    dx = roi1.getBounds().x + shift.x / scale;
                    dy = shift.y / scale;
                    break;
            }
            // Utils.log2("R: " + resultR + " shift: " + shift + " dx,dy: " + dx + ", " + dy);
            return new double[] { dx, dy, success, resultR };
        }
        // new ImagePlus("roi1", ip1.duplicate()).show();
        // new ImagePlus("roi2", ip2.duplicate()).show();
        // try { Thread.sleep(1000000000); } catch (Exception e) {}
        // increase for next iteration
        // increments of 10%
        overlap += 0.10;
    } while (R < min_R && Math.abs(overlap - 1.0f) < 0.001f);
    // Phase-correlation failed, fall back to cross-correlation with a safe overlap
    overlap = percent_overlap * 2;
    if (overlap > 1.0f)
        overlap = 1.0f;
    switch(direction) {
        case TOP_BOTTOM:
            // bottom
            roi1 = new Roi(0, h1 - (int) (h1 * overlap), w1, (int) (h1 * overlap));
            // top
            roi2 = new Roi(0, 0, w2, (int) (h2 * overlap));
            break;
        case LEFT_RIGHT:
            // right
            roi1 = new Roi(w1 - (int) (w1 * overlap), 0, (int) (w1 * overlap), h1);
            // left
            roi2 = new Roi(0, 0, (int) (w2 * overlap), h2);
            break;
    }
    // use one third of the size used for phase-correlation though! Otherwise, it may take FOREVER
    final double scale_cc = scale / 3.0f;
    ip1 = makeStripe(base, roi1, scale_cc);
    ip2 = makeStripe(moving, roi2, scale_cc);
    // gaussian blur them before cross-correlation
    ip1.setPixels(ImageFilter.computeGaussianFastMirror(new FloatArray2D((float[]) ip1.getPixels(), ip1.getWidth(), ip1.getHeight()), 1f).data);
    ip2.setPixels(ImageFilter.computeGaussianFastMirror(new FloatArray2D((float[]) ip2.getPixels(), ip2.getWidth(), ip2.getHeight()), 1f).data);
    // new ImagePlus("CC roi1", ip1).show();
    // new ImagePlus("CC roi2", ip2).show();
    final CrossCorrelation2D cc = new CrossCorrelation2D(ip1, ip2, false);
    double[] cc_result = null;
    switch(direction) {
        case TOP_BOTTOM:
            cc_result = cc.computeCrossCorrelationMT(0.9, 0.3, false);
            break;
        case LEFT_RIGHT:
            cc_result = cc.computeCrossCorrelationMT(0.3, 0.9, false);
            break;
    }
    if (cc_result[2] > min_R / 2) {
        // accepting if R is above half the R accepted for Phase Correlation
        // success
        final int success = SUCCESS;
        switch(direction) {
            case TOP_BOTTOM:
                // boundary checks:
                // if (cc_result[1]/scale_cc > default_dy) success = ERROR;
                dx = cc_result[0] / scale_cc;
                dy = roi1.getBounds().y + cc_result[1] / scale_cc;
                break;
            case LEFT_RIGHT:
                // boundary checks:
                // if (cc_result[0]/scale_cc > default_dx) success = ERROR;
                dx = roi1.getBounds().x + cc_result[0] / scale_cc;
                dy = cc_result[1] / scale_cc;
                break;
        }
        // Utils.log2("\trois: \t" + roi1 + "\n\t\t" + roi2);
        return new double[] { dx, dy, success, cc_result[2] };
    }
    // Utils.log2("Using default");
    return new double[] { default_dx, default_dy, ERROR, 0 };
// / ABOVE: boundary checks don't work if default_dx,dy are zero! And may actually be harmful in anycase
}

Example 2 with FloatArray2D

use of mpi.fruitfly.math.datastructures.FloatArray2D in project TrakEM2 by trakem2.

the class CrossCorrelation2D method computeCrossCorrelation.

/**
 * Computes a translational registration with the help of the cross correlation measure. <br>
 * Limits the overlap to 30% and restricts the vertical shift furthermore by a factor of 16.
 *
 * (NOTE: this method is only single threaded, use computeCrossCorrelationMT instead)
 *
 * @deprecated This method is only single threaded, use computeCrossCorrelationMT instead
 *
 * @param relMinOverlapX double - if you want to scan for less possible translations seen from a direct overlay,
 * give the relative factor here (e.g. 0.3 means DONOT scan the outer 30%)
 * NOTE: Below 0.05 does not really make sense as you then compare only very few pixels (even one) on the edges
 * which gives then an R of 1 (perfect match)
 * @param relMinOverlapY double - if you want to scan for less possible translations seen from a direct overlay,
 * give the relative factor here (e.g. 0.3 means DONOT scan the outer 30%)
 * NOTE: Below 0.05 does not really make sense as you then compare only very few pixels (even one) on the edges
 * which gives then an R of 1 (perfect match)
 * @param showImages boolean - Show the result of the cross correlation translation
 * @return double[] return a double array containing {displaceX, displaceY, R}
 */
@Deprecated
public double[] computeCrossCorrelation(final double relMinOverlapX, final double relMinOverlapY, boolean showImages) {
    double maxR = -2;
    int displaceX = 0, displaceY = 0;
    int w1 = img1.width;
    int w2 = img2.width;
    int h1 = img1.height;
    int h2 = img2.height;
    // int factorY = 16;
    final int min_border_w = (int) (w1 < w2 ? w1 * relMinOverlapX + 0.5 : w2 * relMinOverlapX + 0.5);
    final int min_border_h = (int) (h1 < h2 ? h1 * relMinOverlapY + 0.5 : h2 * relMinOverlapY + 0.5);
    for (int moveX = (-w1 + min_border_w); moveX < (w2 - min_border_w); moveX++) {
        for (int moveY = (-h1 + min_border_h); moveY < (h2 - min_border_h); moveY++) {
            // compute average
            double avg1 = 0, avg2 = 0;
            int count = 0;
            double value1, value2;
            // for (int x1 = 0; x1 < w1; x1++)
            for (int x1 = -min(0, moveX); x1 < min(w1, w2 - moveX); x1++) {
                int x2 = x1 + moveX;
                // for (int y1 = 0; y1 < h1; y1++)
                for (int y1 = -min(0, moveY); y1 < min(h1, h2 - moveY); y1++) {
                    int y2 = y1 + moveY;
                    /*if (y2 < 0 || y2 > h2 - 1)
                            continue;*/
                    value1 = img1.get(x1, y1);
                    if (value1 == -1)
                        continue;
                    value2 = img2.get(x2, y2);
                    if (value2 == -1)
                        continue;
                    avg1 += value1;
                    avg2 += value2;
                    count++;
                }
            }
            if (0 == count)
                continue;
            avg1 /= (double) count;
            avg2 /= (double) count;
            double var1 = 0, var2 = 0;
            double coVar = 0;
            double dist1, dist2;
            // for (int x1 = 0; x1 < w1; x1++)
            for (int x1 = -min(0, moveX); x1 < min(w1, w2 - moveX); x1++) {
                int x2 = x1 + moveX;
                // for (int y1 = 0; y1 < h1; y1++)
                for (int y1 = -min(0, moveY); y1 < min(h1, h2 - moveY); y1++) {
                    int y2 = y1 + moveY;
                    /*if (y2 < 0 || y2 > h2 - 1)
                            continue;*/
                    value1 = img1.get(x1, y1);
                    if (value1 == -1)
                        continue;
                    value2 = img2.get(x2, y2);
                    if (value2 == -1)
                        continue;
                    dist1 = value1 - avg1;
                    dist2 = value2 - avg2;
                    coVar += dist1 * dist2;
                    var1 += dist1 * dist1;
                    var2 += dist2 * dist2;
                }
            }
            var1 /= (double) count;
            var2 /= (double) count;
            coVar /= (double) count;
            double stDev1 = Math.sqrt(var1);
            double stDev2 = Math.sqrt(var2);
            // compute correlation coeffienct
            double R = coVar / (stDev1 * stDev2);
            if (R > maxR) {
                maxR = R;
                displaceX = moveX;
                displaceY = moveY;
            }
            if (R < -1 || R > 1) {
                System.out.println("BIG ERROR! R =" + R);
            }
        }
    // System.out.println(moveX + " [" + (-w2 + min_border_w) + ", " + (w1 - min_border_w) + "] + best R: " + maxR);
    }
    if (showImages) {
        System.out.println(-displaceX + " " + -displaceY);
        FloatArray2D result = drawTranslatedImages(img1, img2, new Point(-displaceX, -displaceY), DRAWTYPE_OVERLAP);
        FloatArrayToImagePlus(result, "result", 0, 0).show();
    }
    return new double[] { -displaceX, -displaceY, maxR };
}

Example 3 with FloatArray2D

use of mpi.fruitfly.math.datastructures.FloatArray2D in project TrakEM2 by trakem2.

the class ImageFilter method createGaussianKernel2D.

public static FloatArray2D createGaussianKernel2D(final float sigma, final boolean normalize) {
    int size = 3;
    FloatArray2D gaussianKernel;
    if (sigma <= 0) {
        gaussianKernel = new FloatArray2D(3, 3);
        gaussianKernel.data[4] = 1;
    } else {
        size = max(3, (int) (2 * (int) (3 * sigma + 0.5) + 1));
        final float two_sq_sigma = 2 * sigma * sigma;
        gaussianKernel = new FloatArray2D(size, size);
        for (int y = size / 2; y >= 0; --y) {
            for (int x = size / 2; x >= 0; --x) {
                final float val = (float) Math.exp(-(float) (y * y + x * x) / two_sq_sigma);
                gaussianKernel.set(val, size / 2 - x, size / 2 - y);
                gaussianKernel.set(val, size / 2 - x, size / 2 + y);
                gaussianKernel.set(val, size / 2 + x, size / 2 - y);
                gaussianKernel.set(val, size / 2 + x, size / 2 + y);
            }
        }
    }
    if (normalize) {
        float sum = 0;
        for (final float value : gaussianKernel.data) sum += value;
        for (int i = 0; i < gaussianKernel.data.length; i++) gaussianKernel.data[i] /= sum;
    }
    return gaussianKernel;
}

Example 4 with FloatArray2D

use of mpi.fruitfly.math.datastructures.FloatArray2D in project TrakEM2 by trakem2.

the class ImageFilter method computeIncreasingGaussianX.

public static FloatArray2D computeIncreasingGaussianX(final FloatArray2D input, final float stDevStart, final float stDevEnd) {
    final FloatArray2D output = new FloatArray2D(input.width, input.height);
    final int width = input.width;
    final float changeFilterSize = (float) (stDevEnd - stDevStart) / (float) width;
    float sigma;
    int filterSize;
    float avg;
    for (int x = 0; x < input.width; x++) {
        sigma = stDevStart + changeFilterSize * (float) x;
        final FloatArray2D kernel = createGaussianKernel2D(sigma, true);
        filterSize = kernel.width;
        for (int y = 0; y < input.height; y++) {
            avg = 0;
            for (int fx = -filterSize / 2; fx <= filterSize / 2; fx++) for (int fy = -filterSize / 2; fy <= filterSize / 2; fy++) {
                try {
                    avg += input.get(x + fx, y + fy) * kernel.get(fx + filterSize / 2, fy + filterSize / 2);
                } catch (final Exception e) {
                }
                ;
            }
            output.set(avg, x, y);
        }
    }
    return output;
}

Example 5 with FloatArray2D

use of mpi.fruitfly.math.datastructures.FloatArray2D in project TrakEM2 by trakem2.

the class ImageFilter method computeLaPlaceFilter5.

/*public static void computeLaPlaceFilterInPlace3(FloatArray2D input)
    {

        float buffer = max(input.height, input.width);

        float derivX, derivY;
        float x1, x2, x3;
        float y1, y2, y3;

        for (int diag = 1; diag < input.width + input.height - 1; diag++)
        {

        }

        for (int y = 1; y < input.height -1 ; y++)
            for (int x = 1; x < input.width -1; x++)
            {
                x1 = input.get(x-1,y);
                x2 = input.get(x,y);
                x3 = input.get(x+1,y);

                derivX = x1 - 2*x2 + x3;

                y1 = input.get(x,y-1);
                y2 = input.get(x,y);
                y3 = input.get(x,y+1);

                derivY = y1 - 2*y2 + y3;

                //output.set((float)Math.sqrt(Math.pow(derivX,2) + Math.pow(derivY,2)), x, y);
            }

        return;
    }*/
public static FloatArray2D computeLaPlaceFilter5(final FloatArray2D input) {
    final FloatArray2D output = new FloatArray2D(input.width, input.height);
    float derivX, derivY;
    float x1, x3, x5;
    float y1, y3, y5;
    for (int y = 2; y < input.height - 2; y++) for (int x = 2; x < input.width - 2; x++) {
        x1 = input.get(x - 2, y);
        x3 = input.get(x, y);
        x5 = input.get(x + 2, y);
        derivX = x1 - 2 * x3 + x5;
        y1 = input.get(x, y - 2);
        y3 = input.get(x, y);
        y5 = input.get(x, y + 2);
        derivY = y1 - 2 * y3 + y5;
        output.set((float) Math.sqrt(Math.pow(derivX, 2) + Math.pow(derivY, 2)), x, y);
    }
    return output;
}