harris - OptiTrust Trace

⌵Trace for harris

›Variable.inline ~simpl [ctx; multi cVarDef ["h1"; "w1"; "h2"; "w2"]]; Stencil.fusion_targets ~nest_of:2 ~outputs ~overlaps [ctx; any cFun ops] in
›List.iter fuse [ ["grayscale"], [], ["gray"]; ["sobelX"; "sobelY"], [], ["ix"; "iy"]; ["mul"; "sum3x3"; "coarsity"], overlaps_2x2 ["ixx"; "ixy"; "iyy"], ["out"]; ];
›Stencil.fusion_targets_tile [int 32] ~outputs:["out"] ~overlaps:["gray", [int 4]; "ix", [int 2]] [ctx; nbMulti; cFor "y"];
›simpl_mins [ctx];
›Matrix.storage_folding ~dim:0 ~size:(int 4) [ctx; multi cVarDef ["gray"; "ix"; "iy"]];
›Matrix.elim [ctx; multi cVarDef ["ixx"; "ixy"; "iyy"]];
›List.iter inline ["weights_sobelX"; "weights_sobelY"; "weights_sum3x3"]; Variable.bind_syntactic ~dest:[ctx; tBefore; cVarDef "acc_sxx"] ~fresh_name:(name ^ "${occ}") [ctx; cArrayRead name] in
›List.iter bind_gradient ["ix"; "iy"];
›Matrix.elim_mops [ctx];
›Omp.parallel_for [ctx; cFor "y"];
›Omp.simd ~clause:[Simdlen 8] [ctx; nbMulti; cFor "x"]; )

tmp_before.cpp → tmp_after.cpp RENAMED Viewed

@@ -67,43 +67,107 @@ void coarsity(float* out, int h, int w, const float* sxx, const float* sxy,
  for (int y = 0; y < h; y++) {
   for (int x = 0; x < w; x++) {
    float det = sxx[y][x] * syy[y][x] - sxy[y][x] * sxy[y][x];
    float trace = sxx[y][x] + syy[y][x];
    out[y][x] = det - kappa * trace * trace;
   }
  }
 }
 void harris(float* out, int h, int w, const float* in) {
- const int h1 = h - 2;
- const int w1 = w - 2;
- const int h2 = h1 - 2;
- const int w2 = w1 - 2;
- float* const gray = (float* const)malloc(sizeof(float[h][w]));
- float* const ix = (float* const)malloc(sizeof(float[h1][w1]));
- float* const iy = (float* const)malloc(sizeof(float[h1][w1]));
- float* const ixx = (float* const)malloc(sizeof(float[h1][w1]));
- float* const ixy = (float* const)malloc(sizeof(float[h1][w1]));
- float* const iyy = (float* const)malloc(sizeof(float[h1][w1]));
- float* const sxx = (float* const)malloc(sizeof(float[h2][w2]));
- float* const sxy = (float* const)malloc(sizeof(float[h2][w2]));
- float* const syy = (float* const)malloc(sizeof(float[h2][w2]));
- grayscale(gray, h, w, in);
- sobelX(ix, h, w, gray);
- sobelY(iy, h, w, gray);
- mul(ixx, h1, w1, ix, ix);
- mul(ixy, h1, w1, ix, iy);
- mul(iyy, h1, w1, iy, iy);
- sum3x3(sxx, h1, w1, ixx);
- sum3x3(sxy, h1, w1, ixy);
- sum3x3(syy, h1, w1, iyy);
- coarsity(out, h2, w2, sxx, sxy, syy, 0.0399999991059f);
   free(gray);
   free(ix);
   free(iy);
- free(ixx);
- free(ixy);
- free(iyy);
- free(sxx);
- free(sxy);
- free(syy);
 }

  for (int y = 0; y < h; y++) {
   for (int x = 0; x < w; x++) {
    float det = sxx[y][x] * syy[y][x] - sxy[y][x] * sxy[y][x];
    float trace = sxx[y][x] + syy[y][x];
    out[y][x] = det - kappa * trace * trace;
   }
  }
 }
 void harris(float* out, int h, int w, const float* in) {
+ #pragma omp parallel for
+ for (int y = 0; y < -4 + h; y += 32) {
+  float* const iy = (float* const)malloc(sizeof(float[4][-2 + w]));
+  float* const ix = (float* const)malloc(sizeof(float[4][-2 + w]));
+  float* const gray = (float* const)malloc(sizeof(float[4][w]));
+  for (int y_out = 0; y_out < min(h, 36 + y) - y; y_out++) {
+   #pragma omp simd simdlen( 8 )
+   for (int x = 0; x < w; x++) {
+    gray[y_out % 4 * w + x] =
+      0.298999994993f * in[(y_out + y) * w + x] +
+      0.587000012398f * in[h * w + (y_out + y) * w + x] +
+      0.11400000006f * in[2 * h * w + (y_out + y) * w + x];
+   }
+   if (2 <= y_out) {
+    #pragma omp simd simdlen( 8 )
+    for (int x = 0; x < -2 + w; x++) {
+     float acc_ix = 0.f;
+     acc_ix += -0.0833333333333f * gray[(-2 + y_out) % 4 * w + x];
+     acc_ix += 0.0833333333333f * gray[2 + (-2 + y_out) % 4 * w + x];
+     acc_ix += -0.166666666667f * gray[(-1 + y_out) % 4 * w + x];
+     acc_ix += 0.166666666667f * gray[2 + (-1 + y_out) % 4 * w + x];
+     acc_ix += -0.0833333333333f * gray[y_out % 4 * w + x];
+     acc_ix += 0.0833333333333f * gray[2 + y_out % 4 * w + x];
+     ix[(-2 + y_out) % 4 * (-2 + w) + x] = acc_ix;
+     float acc_iy = 0.f;
+     acc_iy += -0.0833333333333f * gray[(-2 + y_out) % 4 * w + x];
+     acc_iy += -0.166666666667f * gray[1 + (-2 + y_out) % 4 * w + x];
+     acc_iy += -0.0833333333333f * gray[2 + (-2 + y_out) % 4 * w + x];
+     acc_iy += 0.0833333333333f * gray[y_out % 4 * w + x];
+     acc_iy += 0.166666666667f * gray[1 + y_out % 4 * w + x];
+     acc_iy += 0.0833333333333f * gray[2 + y_out % 4 * w + x];
+     iy[(-2 + y_out) % 4 * (-2 + w) + x] = acc_iy;
+    }
+   }
+   if (4 <= y_out) {
+    #pragma omp simd simdlen( 8 )
+    for (int x = 0; x < -4 + w; x++) {
+     float ix0 = ix[(-4 + y_out) % 4 * (-2 + w) + x];
+     float ix1 = ix[1 + (-4 + y_out) % 4 * (-2 + w) + x];
+     float ix2 = ix[2 + (-4 + y_out) % 4 * (-2 + w) + x];
+     float ix3 = ix[(-3 + y_out) % 4 * (-2 + w) + x];
+     float ix4 = ix[1 + (-3 + y_out) % 4 * (-2 + w) + x];
+     float ix5 = ix[2 + (-3 + y_out) % 4 * (-2 + w) + x];
+     float ix6 = ix[(-2 + y_out) % 4 * (-2 + w) + x];
+     float ix7 = ix[1 + (-2 + y_out) % 4 * (-2 + w) + x];
+     float ix8 = ix[2 + (-2 + y_out) % 4 * (-2 + w) + x];
+     float iy0 = iy[(-4 + y_out) % 4 * (-2 + w) + x];
+     float iy1 = iy[1 + (-4 + y_out) % 4 * (-2 + w) + x];
+     float iy2 = iy[2 + (-4 + y_out) % 4 * (-2 + w) + x];
+     float iy3 = iy[(-3 + y_out) % 4 * (-2 + w) + x];
+     float iy4 = iy[1 + (-3 + y_out) % 4 * (-2 + w) + x];
+     float iy5 = iy[2 + (-3 + y_out) % 4 * (-2 + w) + x];
+     float iy6 = iy[(-2 + y_out) % 4 * (-2 + w) + x];
+     float iy7 = iy[1 + (-2 + y_out) % 4 * (-2 + w) + x];
+     float iy8 = iy[2 + (-2 + y_out) % 4 * (-2 + w) + x];
+     float acc_sxx = 0.f;
+     acc_sxx += ix0 * ix0;
+     acc_sxx += ix1 * ix1;
+     acc_sxx += ix2 * ix2;
+     acc_sxx += ix3 * ix3;
+     acc_sxx += ix4 * ix4;
+     acc_sxx += ix5 * ix5;
+     acc_sxx += ix6 * ix6;
+     acc_sxx += ix7 * ix7;
+     acc_sxx += ix8 * ix8;
+     float acc_sxy = 0.f;
+     acc_sxy += ix0 * iy0;
+     acc_sxy += ix1 * iy1;
+     acc_sxy += ix2 * iy2;
+     acc_sxy += ix3 * iy3;
+     acc_sxy += ix4 * iy4;
+     acc_sxy += ix5 * iy5;
+     acc_sxy += ix6 * iy6;
+     acc_sxy += ix7 * iy7;
+     acc_sxy += ix8 * iy8;
+     float acc_syy = 0.f;
+     acc_syy += iy0 * iy0;
+     acc_syy += iy1 * iy1;
+     acc_syy += iy2 * iy2;
+     acc_syy += iy3 * iy3;
+     acc_syy += iy4 * iy4;
+     acc_syy += iy5 * iy5;
+     acc_syy += iy6 * iy6;
+     acc_syy += iy7 * iy7;
+     acc_syy += iy8 * iy8;
+     float det_out = acc_sxx * acc_syy - acc_sxy * acc_sxy;
+     float trace_out = acc_sxx + acc_syy;
+     out[(-4 + y_out + y) * (-4 + w) + x] =
+       det_out - 0.0399999991059f * trace_out * trace_out;
+    }
+   }
+  }
   free(gray);
   free(ix);
   free(iy);
+ }
 }