Psdt1d pseudocode

From CSclasswiki
Jump to: navigation, search

Here is the code for the 1D transform, which will become the kernel. I'll add comments when I get the chance. However, to turn this into a kernel all we need to do is figure out the memory references. In other words, instead of the the five array arguments our kernel will get a pointer to the shared 2D structure, and its own reference number, and have to compute from that which part of the 2D array to look at.


 // psdt2d pseudocode:
 // input A(i,j) is an MxN array, containing the initial point costs
 // outputs C(i,j), L1(i,j) and L2(i,j) are also MxN arrays
 // C contains the squared distance transformed costs
 // L1 and L2 hold the indices of the closest visible point for each element of C
 // By "visible point" I mean a point where A(i,j) == C(i,j).  
 // Perhaps we should call these something else?  Maybe "roots"?

 // transform columns:
 run kernels k from 1 to M:
   [B(k,1..N),L0(k,1..N)] = dt1d(A(k,1..N))
 end kernels

 // transform rows:
 run kernels k from 1 to N:
   [C(1..M,k),L2(1..M,k)] = dt1d(A(1..M,k))
 end kernels

 // fix L1 according to choice of L2:
 for all (i,j):
   L1(i,j) = L0(L2(i,j),j)  // need to check this is right
 end for
 // can probably do above with kernels too

 // pseudocode for psdt1d:
 // input is 1D array (may have been row or column in original 2D array)
 // I'll call its members V(1..N) here, but in practice we will probably have to 
 // compute the actual indices into the 2D array
 // output is two 1D arrays, I'll call D(1..N) and L(1..N).  Again, these are probably
 // actually indices into a 2D array -- but they are identified by a 1D index.
 // The remaining arrays V and Z are temporary scratch space.
 // Depending on what is most efficient, they may be allocated beforehand or not.

 // Quote from the Felzenszwalb paper, page 6:
 // The horizontal grid location of the i-th parabola in the lower envelope is stored in v[i].
 // The range in which the i-th parabola of the lower envelope is below the others is given by z[i] and z[i +1]. 
 // The variable k keeps track of the number of parabolas in the lower envelope.

 // actual code for psdt1t:
 void psdt1d(double *cost, double *out, double *loc, int *v, double *z, int size_x) {
   int k, q;
   double s; 

   // set up
   k = 0;
   // NRH:  before we scan, there are no parabolas we have identified in the lower envelope
   v[0] = 0;
   z[0] = -mxGetInf();
   z[1] = mxGetInf();
 
   // compute
   for (q=1; q<size_x; q++) {
     // NRH:  here we are looping over all points and seeing if we should include their parabolas in the lower envelope
     s = ((cost[q]+q*q)-(cost[v[k]]+v[k]*v[k]))/(2*(q-v[k]));  // intercept
     // NRH:  the above computes s, the intersection of the parabola at q with the last parabola in the envelope.
     // NRH:  the loop below will now look at all the parabolas currently in the lower envelope.  
     // NRH:  If the intersection point s is before the point where the kth parabola fell on the envelope, then it is entirely
     // NRH:  undercut and should be removed.  
     // NRH:  We need a while loop because we may undercut more than one envelope parabola.
     while (s<=z[k]) {
       k = k-1;  // NRH:  remove the last parabola on the envelope
       s = ((cost[q]+q*q)-(cost[v[k]]+v[k]*v[k]))/(2*(q-v[k]));
       // NRH:  Above computes the intersection s with the new last parabola on the envelope.
     }
     // NRH:  Now we add the current parabola to the envelope, updating k, v[k], and z[k..k+1].
     k = k+1;
     v[k] = q;
     z[k] = s;
     z[k+1] = mxGetInf();
   }

  // NRH:  at this point we have figured out which parabolas fall on the lower envelope (stored in v)
  // NRH:  and the intersection points between them (stored in z)
  // NRH:  now we have to go through and use this to compute the output

   k = 0;
   // NRH:  this loop will go through all the points.  As we go, we will move through the parabolas on the lower envelope
   for (q=0; q<size_x; q++) {
     // NRH:  this advances k until it indexes the envelope parabola corresponding to the current point q
     while (z[k+1]<q) {
       k++;
     }
     // NRH:  This stores the coordinate of the center of the parabola that is on the lower envelope at q -- 
     // NRH:  this is what I was calling the "visible point" above.
     if (loc) {
       loc[q] = v[k]+1;
     }
     // NRH:  this computes the actual lower envelope value at q -- 
     // NRH:  its value is the value at the center plus the distance from the center, squared.
     out[q] = (q-v[k])*(q-v[k])+cost[v[k]];
   }
 }