• Create Account

### #ActualEmergent

Posted 30 September 2012 - 12:15 PM

I just realized you're probably looking for a 2d ray shooting algorithm (whereas in my last post I'd assumed you were trying to shoot rays into a 3d grid). Here; let me just give you some MATLAB code; then, between reading about line drawing algorithms and producing your own C/C++ implementation, you'll figure out how it works.

This function is nearly identical to the one I posted previously; the differences are:
1. "v=precision*[sin(theta); cos(theta)]" instead of "v=x2-x1."
2. Instead of iterating for a fixed number of steps (the length of the line), we iterate until we go out of bounds.

function img = rasterRay2dBresenham(img, x0, theta, color)
% New, improved; conceptually cleaner, and without (obvious) bugs.
n = numel(x0);  % Number of dimensions we're working in.
if(n ~= 2)
error('This function is for 2d rays.\n');
end

precision = 256; % An adjustable parameter.
v = int32(round(precision*[cos(theta); sin(theta)]));
V = max(abs(v)); % Inf-norm of v

% In principle, we are computing the unit (in inf-norm) vector
%	u = (1/V)*v .
% We are splitting the integers in u into a "most significant part," or "high" part, uh
% and a "least significant" or "low" part, ul, so that
%	u = uh + (1/V)*ul .
uh = idivide(v, V, 'fix'); % Round towards zero ('floor' rounds towards negative infinity).
ul = rem(v, V); % Use this, not MATLAB's "mod!"  This gives least significant part, preserving sign.

s = size(img);				 % Used in multidimensional array indexing
cp = int32(cumprod([1 s(1:end-1)]))'; %  to compute a linear index

% Preallocate variables
idx = int32(0);
xh = x0;  % High (most significant) part of position
xl = zeros(n,1, 'int32');  % Low (least significant) part of position
while(true)

if(xh(1) < 1 || xh(1) > size(img,1) || ...
xh(2) < 1 || xh(2) > size(img,2))
break;
end

idx = 1+sum((xh-1).*cp); % Compute linear index into array.
img(idx) = color;

% In principle, we have a point x, and a vector u, and we are just
% doing "x = x + u" over and over in a loop.  Because we have split
% "x" and "u" into "low" and "high" parts, we do this as a "ripple carry"
% sum.
xl = xl + ul;  % Add low parts
for j=1:n	  % Carry any overflow
if(2*abs(xl(j)) >= V) % Overflow [Bresenham (rounds rather than floors/fixes)
% The reason we have some flexibility in when we choose to
% carry (e.g., Bresenham vs. DDA) is that, either way, we
% maintain the loop invariant
%	x = xh + xl/V,
% so we have a sensible representation of our current
% position either way.
xh(j) = xh(j) + sign(xl(j));  % Do the carry
xl(j) = xl(j) - sign(xl(j))*V;  % "Mod" the low part
end
xh(j) = xh(j) + uh(j);  % Add the high parts
end

end
end


Finally, let me give you one important caveat: Bresenham (or DDA) lines do not hit all the cells that the ideal line passes through. So it's not appropriate for use as an intersection routine.

### #1Emergent

Posted 30 September 2012 - 12:09 PM

I just realized you're probably looking for a 2d line drawing algorithm (whereas in my last post I'd assumed you were trying to shoot rays into a 3d grid). Here; let me just give you some MATLAB code; then, between reading about line drawing algorithms and producing your own C/C++ implementation, you'll figure out how it works.

This function is nearly identical to the one I posted previously; the differences are:
1. "v=precision*[sin(theta); cos(theta)]" instead of "v=x2-x1."
2. Instead of iterating for a fixed number of steps (the length of the line), we iterate until we go out of bounds.

function img = rasterRay2dBresenham(img, x0, theta, color)
% New, improved; conceptually cleaner, and without (obvious) bugs.
n = numel(x0);  % Number of dimensions we're working in.
if(n ~= 2)
error('This function is for 2d rays.\n');
end

precision = 256; % An adjustable parameter.
v = int32(round(precision*[cos(theta); sin(theta)]));
V = max(abs(v)); % Inf-norm of v

% In principle, we are computing the unit (in inf-norm) vector
%    u = (1/V)*v .
% We are splitting the integers in u into a "most significant part," or "high" part, uh
% and a "least significant" or "low" part, ul, so that
%    u = uh + (1/V)*ul .
uh = idivide(v, V, 'fix'); % Round towards zero ('floor' rounds towards negative infinity).
ul = rem(v, V); % Use this, not MATLAB's "mod!"  This gives least significant part, preserving sign.

s = size(img);				 % Used in multidimensional array indexing
cp = int32(cumprod([1 s(1:end-1)]))'; %  to compute a linear index

% Preallocate variables
idx = int32(0);
xh = x0;  % High (most significant) part of position
xl = zeros(n,1, 'int32');  % Low (least significant) part of position
while(true)

if(xh(1) < 1 || xh(1) > size(img,1) || ...
xh(2) < 1 || xh(2) > size(img,2))
break;
end

idx = 1+sum((xh-1).*cp); % Compute linear index into array.
img(idx) = color;

% In principle, we have a point x, and a vector u, and we are just
% doing "x = x + u" over and over in a loop.  Because we have split
% "x" and "u" into "low" and "high" parts, we do this as a "ripple carry"
% sum.
xl = xl + ul;  % Add low parts
for j=1:n	  % Carry any overflow
if(2*abs(xl(j)) >= V) % Overflow [Bresenham (rounds rather than floors/fixes)
% The reason we have some flexibility in when we choose to
% carry (e.g., Bresenham vs. DDA) is that, either way, we
% maintain the loop invariant
%    x = xh + xl/V,
% so we have a sensible representation of our current
% position either way.
xh(j) = xh(j) + sign(xl(j));  % Do the carry
xl(j) = xl(j) - sign(xl(j))*V;  % "Mod" the low part
end
xh(j) = xh(j) + uh(j);  % Add the high parts
end

end
end


PARTNERS