# 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;
z = -mxGetInf();
z = 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]];
}
}

```