distance transform

dt.cpp

/*-------------------------------------------------------------------------------------
 * 									dt 1-d
 * -----------------------------------------------------------------------------------*/
template <typename Dtype>
void MessagePassingLayer<Dtype>::dt1d(const Dtype *src, Dtype *dst, int *ptr, int step, int len,
        Dtype a_c, Dtype b_c, Dtype a_p, Dtype b_p, Dtype dshift_c, Dtype dshift_p, int dlen) {
  int   *v = new int[len];
  float *z = new float[len+1];
  int k = 0;
  int q = 0;
  v[0] = 0;
  z[0] = -INF;
  z[1] = +INF;

  for (q = 1; q <= len-1; q++) {
    float s = ( (src[q*step] - src[v[k]*step])
    - b_c * -(q - v[k]) + a_c * (square(q) - square(v[k]))
    - b_p * (q - v[k]) + a_p * (square(q) - square(v[k]))
    + 2*a_c * (q-v[k])*(-dshift_c) + 2*a_p * (q-v[k])*(dshift_p) )
    / ( 2*a_c*(q-v[k]) + 2*a_p*(q-v[k]) );
    while (s <= z[k]) {
      k--;
      s = ( (src[q*step] - src[v[k]*step])
      - b_c * -(q - v[k]) + a_c * (square(q) - square(v[k]))
      - b_p * (q - v[k]) + a_p * (square(q) - square(v[k]))
      + 2*a_c * (q-v[k])*(-dshift_c) + 2*a_p * (q-v[k])*(dshift_p) )
      / ( 2*a_c*(q-v[k]) + 2*a_p*(q-v[k]) );
    }
    k++;
    v[k]   = q;
    z[k]   = s;
    z[k+1] = +INF;
  }

  k = 0;
  for (q = 0; q <= dlen-1; q++) {
    while (z[k+1] < q)
      k++;
    dst[q*step] = src[v[k]*step] + a_c * square(q + dshift_c - v[k]) + b_c * -(q + dshift_c - v[k])
    + a_p * square(q - dshift_p - v[k]) + b_p * (q - dshift_p - v[k]);
    ptr[q*step] = v[k];
  }

  delete [] v;
  delete [] z;
}


/**
 * Distance Transform
 * --------------------------------------
 * Compute distance transform of map vals
 * Input params
 * - vals:  input map
 * - sizx, sizy: size of the inputmap vals
 * - defw_c, defw_p: deformation weight child->parent and parent->child
 * - mean_c, mean_p: child->parent mean / parent->child mean
 * - var_c, var_p: child->parent var / parent->child var (could set as [1, 1])
 * - lenx, leny: ????
 *
 * Output Params
 * - M
 * - Ix, Iy
 */
template <typename Dtype>
void MessagePassingLayer<Dtype>::distance_transform(const Dtype* vals, int sizx, int sizy,
                        const Dtype* defw_c, const Dtype* defw_p,
                        Dtype* mean_c, Dtype* var_c,
                        Dtype* mean_p, Dtype* var_p,
                        int32_t lenx, int32_t leny,
                        Dtype *M, int32_t *Ix, int32_t *Iy
                        ) {
  // ---- deformation weight child->parent
  Dtype ax_c = -defw_c[0] ;            // 2nd order
  Dtype bx_c = -defw_c[1] ;            // 1st order
  Dtype ay_c = -defw_c[2];
  Dtype by_c = -defw_c[3];

  // ---- deformation weight parent->child
  Dtype ax_p = -defw_p[0];
  Dtype bx_p = -defw_p[1];
  Dtype ay_p = -defw_p[2];
  Dtype by_p = -defw_p[3];

  Dtype   *tmpM =  new Dtype[leny*sizx];
  int32_t *tmpIy = new int32_t[leny*sizx];

  for (int x = 0; x < sizx; x++)
    dt1d(vals+x*sizy, tmpM+x*leny, tmpIy+x*leny, 1, sizy,
            ay_c/square(var_c[1]), by_c/var_c[1],
            ay_p/square(var_p[1]), by_p/var_p[1],
            mean_c[1], mean_p[1], leny);

  for (int y = 0; y < leny; y++)
    dt1d(tmpM+y, M+y, Ix+y, leny, sizx,
            ax_c/square(var_c[0]), bx_c/var_c[0],
            ax_p/square(var_p[0]), bx_p/var_p[0],
            mean_c[0], mean_p[0], lenx);

  // get argmins and adjust for matlab indexing from 1
  for (int x = 0; x < lenx; x++) {
    for (int y = 0; y < leny; y++) {
      int p = x*leny+y;
      Iy[p] = tmpIy[Ix[p]*leny+y] + 1;
      Ix[p] = Ix[p] + 1;
    }
  }
}