Examples with BlockStatistics mpicbg.ij.integral.BlockStatistics used on opensource projects

Example 1 with BlockStatistics

use of mpicbg.ij.integral.BlockStatistics in project TrakEM2 by trakem2.

the class RobustNormalizeLocalContrast method run.

public static final void run(final ImageProcessor ip, final int scaleLevel, final int brx1, final int bry1, final float stds1, final int brx2, final int bry2, final float stds2) {
    final PixelSetter setter;
    if (ByteProcessor.class.isInstance(ip))
        setter = new ByteSetter((ByteProcessor) ip);
    else if (ShortProcessor.class.isInstance(ip))
        setter = new ShortSetter((ShortProcessor) ip);
    else
        setter = new FloatSetter(ip);
    final int scale = (int) Util.pow(2, scaleLevel);
    final int scale2 = scale / 2;
    final double[] f = new double[scale];
    for (int i = 0; i < scale; ++i) f[i] = (i + 0.5) / scale;
    final int sbrx1 = brx1 / scale;
    final int sbry1 = bry1 / scale;
    final int sbrx2 = brx2 / scale;
    final int sbry2 = bry2 / scale;
    final int ipWidth = ip.getWidth();
    final int ipHeight = ip.getHeight();
    final double ipMin = ip.getMin();
    final double ipLength = ip.getMax() - ipMin;
    FloatProcessor ipScaled = (FloatProcessor) ip.convertToFloat();
    for (int i = 0; i < scaleLevel; ++i) ipScaled = Downsampler.downsampleFloatProcessor(ipScaled);
    final FloatProcessor std = ipScaled;
    // new ImagePlus("downsampled", std.duplicate() ).show();
    RemoveOutliers.run(std, sbrx1, sbry1, stds1);
    // new ImagePlus("outlier-filtered", std.duplicate() ).show();
    final BlockStatistics bs = new BlockStatistics(std);
    bs.mean(sbrx2, sbry2);
    // new ImagePlus("mean", std.duplicate() ).show();
    final FloatProcessor mean = (FloatProcessor) std.duplicate();
    bs.std(sbrx2, sbry2);
    // new ImagePlus("std", std.duplicate() ).show();
    final int w = mean.getWidth();
    final int h = mean.getHeight();
    final ExecutorService exec = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors());
    final ArrayList<Future<?>> tasks = new ArrayList<Future<?>>();
    /* the big inside */
    for (int y = 1; y < h; ++y) {
        final int ya = y - 1;
        final int yb = y;
        tasks.add(exec.submit(new Runnable() {

            public final void run() {
                for (int x = 1; x < w; ++x) {
                    final int xa = x - 1;
                    final float meanA = mean.getf(xa, ya);
                    final float meanB = mean.getf(x, ya);
                    final float meanC = mean.getf(xa, yb);
                    final float meanD = mean.getf(x, yb);
                    final float stdA = std.getf(xa, ya);
                    final float stdB = std.getf(x, ya);
                    final float stdC = std.getf(xa, yb);
                    final float stdD = std.getf(x, yb);
                    int ys = yb * scale - scale2;
                    final int xss = x * scale - scale2;
                    for (int yi = 0; yi < scale; ++ys, ++yi) for (int xs = xss, xi = 0; xi < scale; ++xs, ++xi) {
                        final double meanAB = meanA * f[xi] + meanB * (1.0 - f[xi]);
                        final double meanCD = meanC * f[xi] + meanD * (1.0 - f[xi]);
                        final double meanABCD = meanAB * f[yi] + meanCD * (1.0 - f[yi]);
                        final double stdAB = stdA * f[xi] + stdB * (1.0 - f[xi]);
                        final double stdCD = stdC * f[xi] + stdD * (1.0 - f[xi]);
                        final double stdABCD = stdAB * f[yi] + stdCD * (1.0 - f[yi]);
                        final double d = stds2 * stdABCD;
                        final double min = meanABCD - d;
                        final int i = ys * ipWidth + xs;
                        setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
                    }
                }
            }
        }));
    }
    for (Future<?> task : tasks) {
        try {
            task.get();
        } catch (InterruptedException e) {
            exec.shutdownNow();
            return;
        } catch (ExecutionException e) {
            exec.shutdownNow();
            return;
        }
    }
    tasks.clear();
    exec.shutdown();
    /* top and bottom */
    for (int x = 1; x < w; ++x) {
        final int xa = x - 1;
        final float meanA = mean.getf(xa, 0);
        final float meanB = mean.getf(x, 0);
        final float meanC = mean.getf(xa, h - 1);
        final float meanD = mean.getf(x, h - 1);
        final float stdA = std.getf(xa, 0);
        final float stdB = std.getf(x, 0);
        final float stdC = std.getf(xa, h - 1);
        final float stdD = std.getf(x, h - 1);
        /* top */
        final int xss = x * scale - scale2;
        for (int ys = 0; ys < scale2; ++ys) for (int xs = xss, xi = 0; xi < scale; ++xs, ++xi) {
            final double meanAB = meanA * f[xi] + meanB * (1.0 - f[xi]);
            final double stdAB = stdA * f[xi] + stdB * (1.0 - f[xi]);
            final double d = stds2 * stdAB;
            final double min = meanAB - d;
            final int i = ys * ipWidth + xs;
            setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
        }
        /* bottom */
        for (int ys = h * scale - scale2; ys < ipHeight; ++ys) for (int xs = xss, xi = 0; xi < scale; ++xs, ++xi) {
            final double meanCD = meanC * f[xi] + meanD * (1.0 - f[xi]);
            final double stdCD = stdC * f[xi] + stdD * (1.0 - f[xi]);
            final double d = stds2 * stdCD;
            final double min = meanCD - d;
            final int i = ys * ipWidth + xs;
            setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
        }
    }
    final int xss = w * scale - scale2;
    /* left and right */
    for (int y = 1; y < h; ++y) {
        final int ya = y - 1;
        final float meanA = mean.getf(0, ya);
        final float meanB = mean.getf(0, y);
        final float meanC = mean.getf(w - 1, ya);
        final float meanD = mean.getf(w - 1, y);
        final float stdA = std.getf(0, ya);
        final float stdB = std.getf(0, y);
        final float stdC = std.getf(w - 1, ya);
        final float stdD = std.getf(w - 1, y);
        int ys = y * scale - scale2;
        for (int yi = 0; yi < scale; ++ys, ++yi) {
            /* left */
            for (int xs = 0; xs < scale2; ++xs) {
                final double meanAB = meanA * f[yi] + meanB * (1.0 - f[yi]);
                final double stdAB = stdA * f[yi] + stdB * (1.0 - f[yi]);
                final double d = stds2 * stdAB;
                final double min = meanAB - d;
                final int i = ys * ipWidth + xs;
                setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
            }
            /* right */
            for (int xs = xss; xs < ipWidth; ++xs) {
                final double meanCD = meanC * f[yi] + meanD * (1.0 - f[yi]);
                final double stdCD = stdC * f[yi] + stdD * (1.0 - f[yi]);
                final double d = stds2 * stdCD;
                final double min = meanCD - d;
                final int i = ys * ipWidth + xs;
                setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
            }
        }
    }
    /* corners */
    final float meanA = mean.getf(0, 0);
    final float meanB = mean.getf(w - 1, 0);
    final float meanC = mean.getf(0, h - 1);
    final float meanD = mean.getf(w - 1, h - 1);
    final float stdA = std.getf(0, 0);
    final float stdB = std.getf(w - 1, 0);
    final float stdC = std.getf(0, h - 1);
    final float stdD = std.getf(w - 1, h - 1);
    for (int ys = 0; ys < scale2; ++ys) {
        for (int xs = 0; xs < scale2; ++xs) {
            final double d = stds2 * stdA;
            final double min = meanA - d;
            final int i = ys * ipWidth + xs;
            setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
        }
        for (int xs = xss; xs < ipWidth; ++xs) {
            final double d = stds2 * stdB;
            final double min = meanB - d;
            final int i = ys * ipWidth + xs;
            setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
        }
    }
    for (int ys = h * scale - scale2; ys < ipHeight; ++ys) {
        for (int xs = 0; xs < scale2; ++xs) {
            final double d = stds2 * stdC;
            final double min = meanC - d;
            final int i = ys * ipWidth + xs;
            setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
        }
        for (int xs = xss; xs < ipWidth; ++xs) {
            final double d = stds2 * stdD;
            final double min = meanD - d;
            final int i = ys * ipWidth + xs;
            setter.setf(i, (float) ((ip.getf(i) - min) / 2 / d * ipLength + ipMin));
        }
    }
// new ImagePlus( "processed", ip.duplicate() ).show();
}