このforループの計算速度を最適化する方法は？

私は（少なくとも私にとっては）難しい問題に悩まされています。コードをプロファイリングしたときに、（シングルコアの）計算時間のほぼすべてが、以下の単一のネストされたループ（画像の二重積分）によって消費されていることに気付きました。コンピューティングを加速するための最良の方法は何だと思いますか？

ネストされたストリームにマップしようとしましたが、複数のifブロックをマップする方法がわかりません... OpenCLを使用してGPUでマップしようとすると、問題により適していますか？

ipはImageJImageProcessorであり、その方法.getPixelValue(x,y)もかなりのリソースを消費します。ただし、確立されたlibに属しているため、可能であれば変更は避けたいと思います。

変数宣言：

private ImageProcessor ip = null; //This type comes from ImageJ
private double area;
private double a11, a22;
private double u1, u2;
private double v1, v2;
private double y1, y2;
private static final double HALF_SQRT2 = sqrt(2.0) / 2.0;
private static final double SQRT_TINY =
    sqrt((double)Float.intBitsToFloat((int)0x33FFFFFF));

関数：

private double contrast (
) {
    if (area < 1.0) {
        return(1.0 / SQRT_TINY);
    }
    double c = 0.0;
    final int xmin = max((int)floor(u1), 0);
    final int xmax = min((int)ceil(v1), width - 1);
    final int ymin = max((int)floor(u2), 0);
    final int ymax = min((int)ceil(v2), height - 1);
    if ((u1 < xmin) || (xmax < v1) || (u2 < ymin) || (ymax < v2)){
        return(1.0 / SQRT_TINY);
    }
    if ((xmax <= xmin) || (ymax <= ymin)) {
        return(1.0 / SQRT_TINY);
    }
    for (int y = ymin; (y <= ymax); y++) {
        final double dy = y2 - (double)y;
        final double dy2 = dy * dy;
        for (int x = xmin; (x <= xmax); x++) {
            final double dx = y1 - (double)x;
            final double dx2 = dx * dx;
            final double d = sqrt(dx2 + dy2);
            double z = a11 * dx2 + a12 * dx * dy + a22 * dy2;
            if (z < SQRT_TINY) {
                c -= ip.getPixelValue(x, y);
                continue;
            }
            z = a3 / sqrt(z);
            double d0 = (1.0 - z / SQRT2) * d;
            if (d0 < -HALF_SQRT2) {
                c -= ip.getPixelValue(x, y);
                continue;
            }
            if (d0 < HALF_SQRT2) {
                c += SQRT2 * d0 * ip.getPixelValue(x, y);
                continue;
            }
            d0 = (1.0 - z) * d;
            if (d0 < -1.0) {
                c += ip.getPixelValue(x, y);
                continue;
            }
            if (d0 < 1.0) {
                c += (1.0 - d0) * ip.getPixelValue(x, y) / 2.0;
                continue;
            }
        }
    }
    return(c / area);