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.