BufferedImage img = ImageIO.read(new File("model3.jpg"));
        int pixelsize = img.getWidth() * img.getHeight();

        double[][][] YCbCr = new double[img.getWidth()][img.getHeight()][3];
        double Mr = 0, Mb = 0, Ymax = 0;
        for (int i = 0; i < img.getWidth(); i++) {
            for (int j = 0; j < img.getHeight(); j++) {
                YCbCr[i][j] = toYCbCr(img.getRGB(i, j));
                Mr += YCbCr[i][j][2];
                Mb += YCbCr[i][j][1];
                Ymax = Math.max(Ymax, YCbCr[i][j][0]);
            }
        }

        Mr /= pixelsize;
        Mb /= pixelsize;

        double Dr = 0, Db = 0;
        for (int i = 0; i < YCbCr.length; i++) {
            for (int j = 0; j < YCbCr[i].length; j++) {
                Db += Math.abs(YCbCr[i][j][1] - Mb);
                Dr += Math.abs(YCbCr[i][j][2] - Mr);
            }
        }
        Dr /= pixelsize;
        Db /= pixelsize;

        double[][] Y = new double[img.getWidth()][img.getHeight()];
        double[] Yhistogram = new double[256];
        double Ysum = 0;
        for (int i = 0; i < Y.length; i++) {
            for (int j = 0; j < Y[i].length; j++) {
                int value = (Math.abs(YCbCr[i][j][1] - (Mb + Db * Math.signum(Mb))) < 1.5 * Db) & //
                        (Math.abs(YCbCr[i][j][2]) - (1.5 * Mr + Dr * Math.signum(Mr))) < 1.5 * Dr ? 1 : 0;
                if (value <= 0)
                    continue;
                double y = YCbCr[i][j][0];
                Y[i][j] = y;
                Yhistogram[(int) Y[i][j]]++;
                Ysum++;
            }
        }

        double Yhistogramsum = 0;
        double Ymin = 0;
        for (int i = Yhistogram.length - 1; i >= 0; i--) {
            Yhistogramsum += Yhistogram[i];
            if (Yhistogramsum > 0.1 * Ysum) {
                Ymin = i;
                break;
            }
        }

        double Raver = 0, Gaver = 0, Baver = 0;
        double averSum = 0;
        for (int i = 0; i < Y.length; i++) {
            for (int j = 0; j < Y[i].length; j++) {
                if (Y[i][j] > Ymin) {

                    int color = img.getRGB(i, j);
                    int r = (color >> 16) & 0xFF;
                    int g = (color >> 8) & 0xFF;
                    int b = color & 0xFF;
                    Raver += r;
                    Gaver += g;
                    Baver += b;
                    averSum++;
                }
            }
        }
        Raver /= averSum;
        Gaver /= averSum;
        Baver /= averSum;

        double Rgain = Ymax / Raver, Ggain = Ymax / Gaver, Bgain = Ymax / Baver;
        for (int i = 0; i < img.getWidth(); i++) {
            for (int j = 0; j < img.getHeight(); j++) {
                Color color = new Color(img.getRGB(i, j));
                int r = ensureColor((int) Math.floor(color.getRed() * Rgain));
                int g = ensureColor((int) Math.floor(color.getGreen() * Ggain));
                int b = ensureColor((int) Math.floor(color.getBlue() * Bgain));
                img.setRGB(i, j, new Color(r, g, b).getRGB());
            }
        }

        ImageIO.write(img, "jpg", new File("xxx.jpg"));

    // https://mathematica.stackexchange.com/questions/29786/how-to-convert-rgb-to-ycbcr
    private double[] toYCbCr(int color) {

        int r = (color >> 16) & 0xFF;
        int g = (color >> 8) & 0xFF;
        int b = color & 0xFF;

        double Y = 16 + (65.481 * r / 255 + 128.553 * g / 255 + 24.966 * b / 255);
        double Cb = 128 + (-37.797 * r / 255 - 74.203 * g / 255 + 112 * b / 255);
        double Cr = 128 + (112 * r / 255 - 93.786 * g / 255 - 18.214 * b / 255);

        return new double[] { Y, Cb, Cr };
    }

    private int ensureColor(double color) {
        if (color < 0)
            return 0;
        if (color > 255)
            return 255;
        return (int) color;
    }