# Psdt1d pseudocode

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]]; } }