Matlab - Propagate points orthogonally on to the edge of shape boundaries

I have a set of points which I want to propagate on to the edge of shape boundary defined by a binary image. The shape boundary is defined by a 1px wide white edge.

I have the coordinates of these points stored in a 2 row by n column matrix. The shape forms a concave boundary with no holes within itself made of around 2500 points. I have approximately 80 to 150 points that I wish to propagate on the shape boundary.

I want to cast a ray from each point from the set of points in an orthogonal direction and detect at which point it intersects the shape boundary at. The orthogonal direction has already been determined. For the required purposes it is calculated taking the normal of the contour calculated for point, using point-1 and point+1.

What would be the best method to do this? Are there some sort of ray tracing algorithms that could be used?

Thank you very much in advance for any help!

EDIT: I have tried to make the question much clearer and added a image describing the problem. In the image the grey line represents the shape contour, the red dots the points I want to propagate and the green line an imaginary orthongally cast ray.

alt text http://img504.imageshack.us/img504/3107/orth.png

ANOTHER EDIT: For clarification I have posted the code used to calculate the normals for each point. Where the xt and yt are vectors storing the coordinates for each point. After calculating the normal value it can be propagated by using the linspace function and the requested length of the orthogonal line.

%#derivaties of contour
dx=[xt(2)-xt(1) (xt(3:end)-xt(1:end-2))/2 xt(end)-xt(end-1)];
dy=[yt(2)-yt(1) (yt(3:end)-yt(1:end-2))/2 yt(end)-yt(end-1)];

%#normal开发者_如何转开发s of contourpoints
l=sqrt(dx.^2+dy.^2);
nx = -dy./l; 
ny =  dx./l;

normals = [nx,ny];

It depends on how many unit vectors you want to test against one shape. If you have one shape and many tests, the easiest thing to do is probably to convert your shape coordinates to polar coordinates which implicitly represent your solution already. This may not be a very effective solution however if you have different shapes and only a few tests for every shape.

Update based on the edited question:

If the rays can start from arbitrary points, not only from the origin, you have to test against all the points. This can be done easily by transforming your shape boundary such that your ray to test starts in the origin in either coordinate direction (positive x in my example code)

% vector of shape boundary points (assumed to be image coordinates, i.e. integers)
shapeBoundary = [xs, ys];

% define the start point and direction you want to test
startPoint = [xsp, ysp];
testVector = unit([xv, yv]);

% now transform the shape boundary
shapeBoundaryTrans(:,1) = shapeBoundary(:,1)-startPoint(1);
shapeBoundaryTrans(:,2) = shapeBoundary(:,2)-startPoint(2);
rotMatrix = [testVector(2), testVector(1); ...
             testVector(-1), testVector(2)];
% somewhat strange transformation to keep it vectorized
shapeBoundaryTrans = shapeBoundaryTrans * rotMatrix';

% now the test is easy: find the points close to the positive x-axis
selector = (abs(shapeBoundaryTrans(:,2)) < 0.5) & (shapeBoundaryTrans(:,1) > 0);
shapeBoundaryTrans(:,2) = 1:size(shapeBoundaryTrans, 1)';
shapeBoundaryReduced = shapeBoundaryTrans(selector, :);
if (isempty(shapeBoundaryReduced))
    [dummy, idx] = min(shapeBoundaryReduced(:, 1));
    collIdx = shapeBoundaryReduced(idx, 2);
    % you have a collision with point collIdx of your shapeBoundary
else
    % no collision
end

This could be done in a nicer way probably, but you get the idea...

If I understand your problem correctly (project each point onto the closest point of the shape boundary), you can

1. use sub2ind to convert the "2 row by n column matrix" description to a BW image with white pixels, something like

   myimage=zeros(imagesize); myimage(imagesize, x_coords, y_coords) = 1

2. use imfill to fill the outside of the boundary

3. run [D,L] = bwdist(BW) on the resulting image, and just read the answers from L.

Should be fairly straightforward.