You have two points,

p1 = (x1, y1, z1), and

p2 = (x2, y2, z2) .

You compute a vector

v = p2 - p1 = (x2-x1, y2-y1, z2-z1) .

from the first point to the second.

You scale this vector so that it is one on its largest (in absolute value) dimension, to get a new vector,

s = [ 1 / max( abs(v) ) ] v

(I.e., you divide each element of v by the largest absolute value found in v). This is how much you will step in each dimension on each iteration. The reason for this choice of scaling is to avoid gaps (if you stepped by more than one pixel in any dimension, you would get gaps).

Finally, you execute the following loop (pseudocode):

p = p1 for i = 1 to max(abs(v)) { image(round(p)) = color p = p + s } // Now p = p2.

That's the basic idea.

DDA does this with integer math. I wouldn't worry about it until I'd implemented the floating-point version, but, I'll write a little about it here and include some MATLAB code at the end in case you're curious.

DDA, conceptually, just adds the vector "s" to the vector "p" over and over, but it splits all the numbers into two parts, which you should think of as being like "digits:"

- DDA stores everything to the left of the decimal point in units of 1 pixel.

- It stores everything to the right of the decimal point in units of 1/V-ths of a pixel, where V = max(abs(v)).

- It implements addition between these "mixed base" numbers, using the standard grade-school algorithm: You add the low parts, carry any overflow, and add the high parts.

Here's some (MATLAB) code [EDIT: fixed code; added comments]:

function img = rasterLineB(img, x1, x2, color) % New, improved; conceptually cleaner, and without (obvious) bugs. n = numel(x1); % Number of dimensions we're working in. v = x2 - x1; % vector from x1 to x2 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 = x1; % High (most significant) part of position xl = zeros(n,1, 'int32'); % Low (least significant) part of position for k=1:V+1 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(abs(xl(j)) >= V) % Overflow 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