MultiSource/Benchmarks/Rodinia/srad/sradKernel.c

#include "srad.h"
void srad_kernel(float dN[ROWS][COLS], float dS[ROWS][COLS],
                 float dW[ROWS][COLS], float dE[ROWS][COLS],
                 float I[ROWS][COLS], float J[ROWS][COLS],
                 float c[ROWS][COLS]) {
  int size_R, iter;
  float q0sqr, sum, sum2, tmp, meanROI, varROI;
  float Jc, G2, L, num, den, qsqr;
  float cN, cS, cW, cE;
  float D;
  int i, j;

  size_R = (Y2 - Y1 + 1) * (X2 - X1 + 1);

  for (iter = 0; iter < ITER; iter++) {
    sum = 0;
    sum2 = 0;
    for (i = Y1; i <= Y2; i++) {
      for (j = X1; j <= X2; j++) {
        tmp = J[i][j];
        sum += tmp;
        sum2 += tmp * tmp;
      }
    }
    meanROI = sum / size_R;
    varROI = (sum2 / size_R) - meanROI * meanROI;
    q0sqr = varROI / (meanROI * meanROI);

    {
      int i = 0;
      int j = 0;
      Jc = J[i][j];

      // directional derivates
      dN[i][j] = J[0][j] - Jc;
      dS[i][j] = J[1][j] - Jc;
      dW[i][j] = J[i][0] - Jc;
      dE[i][j] = J[i][1] - Jc;

      G2 = (dN[i][j] * dN[i][j] + dS[i][j] * dS[i][j] + dW[i][j] * dW[i][j] +
            dE[i][j] * dE[i][j]) /
           (Jc * Jc);

      L = (dN[i][j] + dS[i][j] + dW[i][j] + dE[i][j]) / Jc;

      num = (0.5 * G2) - ((1.0 / 16.0) * (L * L));
      den = 1 + (.25 * L);
      qsqr = num / (den * den);

      // diffusion coefficent (equ 33)
      den = (qsqr - q0sqr) / (q0sqr * (1 + q0sqr));
      c[i][j] = 1.0 / (1.0 + den);

      // saturate diffusion coefficent
      if (c[i][j] < 0) {
        c[i][j] = 0;
      } else if (c[i][j] > 1) {
        c[i][j] = 1;
      }

      j = COLS - 1;

      Jc = J[i][j];
      // directional derivates
      dN[i][j] = J[0][j] - Jc;
      dS[i][j] = J[1][j] - Jc;
      dW[i][j] = J[i][COLS - 2] - Jc;
      dE[i][j] = J[i][COLS - 1] - Jc;

      G2 = (dN[i][j] * dN[i][j] + dS[i][j] * dS[i][j] + dW[i][j] * dW[i][j] +
            dE[i][j] * dE[i][j]) /
           (Jc * Jc);

      L = (dN[i][j] + dS[i][j] + dW[i][j] + dE[i][j]) / Jc;

      num = (0.5 * G2) - ((1.0 / 16.0) * (L * L));
      den = 1 + (.25 * L);
      qsqr = num / (den * den);

      // diffusion coefficent (equ 33)
      den = (qsqr - q0sqr) / (q0sqr * (1 + q0sqr));
      c[i][j] = 1.0 / (1.0 + den);

      // saturate diffusion coefficent
      if (c[i][j] < 0) {
        c[i][j] = 0;
      } else if (c[i][j] > 1) {
        c[i][j] = 1;
      }
    }
    {
      int i = ROWS - 1;
      int j = 0;
      Jc = J[i][j];

      // directional derivates
      dN[i][j] = J[ROWS - 2][j] - Jc;
      dS[i][j] = J[ROWS - 1][j] - Jc;
      dW[i][j] = J[i][0] - Jc;
      dE[i][j] = J[i][1] - Jc;

      G2 = (dN[i][j] * dN[i][j] + dS[i][j] * dS[i][j] + dW[i][j] * dW[i][j] +
            dE[i][j] * dE[i][j]) /
           (Jc * Jc);

      L = (dN[i][j] + dS[i][j] + dW[i][j] + dE[i][j]) / Jc;

      num = (0.5 * G2) - ((1.0 / 16.0) * (L * L));
      den = 1 + (.25 * L);
      qsqr = num / (den * den);

      // diffusion coefficent (equ 33)
      den = (qsqr - q0sqr) / (q0sqr * (1 + q0sqr));
      c[i][j] = 1.0 / (1.0 + den);

      // saturate diffusion coefficent
      if (c[i][j] < 0) {
        c[i][j] = 0;
      } else if (c[i][j] > 1) {
        c[i][j] = 1;
      }

      j = COLS - 1;

      Jc = J[i][j];
      // directional derivates
      dN[i][j] = J[ROWS - 2][j] - Jc;
      dS[i][j] = J[ROWS - 1][j] - Jc;
      dW[i][j] = J[i][COLS - 2] - Jc;
      dE[i][j] = J[i][COLS - 1] - Jc;

      G2 = (dN[i][j] * dN[i][j] + dS[i][j] * dS[i][j] + dW[i][j] * dW[i][j] +
            dE[i][j] * dE[i][j]) /
           (Jc * Jc);

      L = (dN[i][j] + dS[i][j] + dW[i][j] + dE[i][j]) / Jc;

      num = (0.5 * G2) - ((1.0 / 16.0) * (L * L));
      den = 1 + (.25 * L);
      qsqr = num / (den * den);

      // diffusion coefficent (equ 33)
      den = (qsqr - q0sqr) / (q0sqr * (1 + q0sqr));
      c[i][j] = 1.0 / (1.0 + den);

      // saturate diffusion coefficent
      if (c[i][j] < 0) {
        c[i][j] = 0;
      } else if (c[i][j] > 1) {
        c[i][j] = 1;
      }
    }

    for (int i = 1; i < ROWS - 1; i++) {
      for (int j = 1; j < COLS - 1; j++) {

        Jc = J[i][j];

        // directional derivates
        dN[i][j] = J[i - 1][j] - Jc;
        dS[i][j] = J[i + 1][j] - Jc;
        dW[i][j] = J[i][j - 1] - Jc;
        dE[i][j] = J[i][j + 1] - Jc;

        G2 = (dN[i][j] * dN[i][j] + dS[i][j] * dS[i][j] + dW[i][j] * dW[i][j] +
              dE[i][j] * dE[i][j]) /
             (Jc * Jc);

        L = (dN[i][j] + dS[i][j] + dW[i][j] + dE[i][j]) / Jc;

        num = (0.5 * G2) - ((1.0 / 16.0) * (L * L));
        den = 1 + (.25 * L);
        qsqr = num / (den * den);

        // diffusion coefficent (equ 33)
        den = (qsqr - q0sqr) / (q0sqr * (1 + q0sqr));
        c[i][j] = 1.0 / (1.0 + den);

        // saturate diffusion coefficent
        if (c[i][j] < 0) {
          c[i][j] = 0;
        } else if (c[i][j] > 1) {
          c[i][j] = 1;
        }
      }
    }

    {
      int i = ROWS - 1;
      for (int j = 0; j < COLS - 1; j++) {
        // diffusion coefficent
        cN = c[i][j];
        cS = c[i][j];
        cW = c[i][j];
        cE = c[i][j + 1];

        // divergence (equ 58)
        D = cN * dN[i][j] + cS * dS[i][j] + cW * dW[i][j] + cE * dE[i][j];

        // image update (equ 61)
        J[i][j] = J[i][j] + 0.25 * LAMDBA * D;
      }
      i = ROWS - 1;
      int j = COLS - 1;

      // diffusion coefficent
      cN = c[i][j];
      cS = c[i][j];
      cW = c[i][j];
      cE = c[i][j];

      // divergence (equ 58)
      D = cN * dN[i][j] + cS * dS[i][j] + cW * dW[i][j] + cE * dE[i][j];

      // image update (equ 61)
      J[i][j] = J[i][j] + 0.25 * LAMDBA * D;
    }

    for (int i = 0; i < ROWS - 1; i++) {
      for (int j = 0; j < COLS - 1; j++) {
        // diffusion coefficent
        cN = c[i][j];
        cS = c[i + 1][j];
        cW = c[i][j];
        cE = c[i][j + 1];

        // divergence (equ 58)
        D = cN * dN[i][j] + cS * dS[i][j] + cW * dW[i][j] + cE * dE[i][j];

        // image update (equ 61)
        J[i][j] = J[i][j] + 0.25 * LAMDBA * D;
      }
    }
  }
}