Jump to content

  • Log In with Google      Sign In   
  • 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)
			%if(abs(xl(j)) >= V)  % Overflow [DDA (grade school addition)]
				% 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)
		    %if(abs(xl(j)) >= V)  % Overflow [DDA (grade school addition)]
			    % 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