NIFTI_20090325/bresenham_line3d.m
87b08fd7
 %  Generate X Y Z coordinates of a 3D Bresenham's line between
 %  two given points.
 %
 %  A very useful application of this algorithm can be found in the
 %  implementation of Fischer's Bresenham interpolation method in my
 %  another program that can rotate three dimensional image volume
 %  with an affine matrix:
 %  http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21080
 %
 %  Usage: [X Y Z] = bresenham_line3d(P1, P2, [precision]);
 %
 %  P1	- vector for Point1, where P1 = [x1 y1 z1]
 %
 %  P2	- vector for Point2, where P2 = [x2 y2 z2]
 %
 %  precision (optional) - Although according to Bresenham's line
 %	algorithm, point coordinates x1 y1 z1 and x2 y2 z2 should
 %	be integer numbers, this program extends its limit to all
 %	real numbers. If any of them are floating numbers, you
 %	should specify how many digits of decimal that you would
 %	like to preserve. Be aware that the length of output X Y
 %	Z coordinates will increase in 10 times for each decimal
 %	digit that you want to preserve. By default, the precision
 %	is 0, which means that they will be rounded to the nearest
 %	integer.
 %
 %  X	- a set of x coordinates on Bresenham's line
 %
 %  Y	- a set of y coordinates on Bresenham's line
 %
 %  Z	- a set of z coordinates on Bresenham's line
 %
 %  Therefore, all points in XYZ set (i.e. P(i) = [X(i) Y(i) Z(i)])
 %  will constitute the Bresenham's line between P1 and P1.
 %
 %  Example:
 %	P1 = [12 37 6];     P2 = [46 3 35];
 %	[X Y Z] = bresenham_line3d(P1, P2);
 %	figure; plot3(X,Y,Z,'s','markerface','b');
 %
 %  This program is ported to MATLAB from:
 %
 %  B.Pendleton.  line3d - 3D Bresenham's (a 3D line drawing algorithm)
 %  ftp://ftp.isc.org/pub/usenet/comp.sources.unix/volume26/line3d, 1992
 %
 %  Which is also referenced by:
 %
 %  Fischer, J., A. del Rio (2004).  A Fast Method for Applying Rigid
 %  Transformations to Volume Data, WSCG2004 Conference.
 %  http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
 %
 %  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
 %
 function [X,Y,Z] = bresenham_line3d(P1, P2, precision)
 
    if ~exist('precision','var') | isempty(precision) | round(precision) == 0
       precision = 0;
       P1 = round(P1);
       P2 = round(P2);
    else
       precision = round(precision);
       P1 = round(P1*(10^precision));
       P2 = round(P2*(10^precision));
    end
 
    d = max(abs(P2-P1)+1);
    X = zeros(1, d);
    Y = zeros(1, d);
    Z = zeros(1, d);
 
    x1 = P1(1);
    y1 = P1(2);
    z1 = P1(3);
 
    x2 = P2(1);
    y2 = P2(2);
    z2 = P2(3);
 
    dx = x2 - x1;
    dy = y2 - y1;
    dz = z2 - z1;
 
    ax = abs(dx)*2;
    ay = abs(dy)*2;
    az = abs(dz)*2;
 
    sx = sign(dx);
    sy = sign(dy);
    sz = sign(dz);
 
    x = x1;
    y = y1;
    z = z1;
    idx = 1;
 
    if(ax>=max(ay,az))			% x dominant
       yd = ay - ax/2;
       zd = az - ax/2;
 
       while(1)
          X(idx) = x;
          Y(idx) = y;
          Z(idx) = z;
          idx = idx + 1;
 
          if(x == x2)		% end
             break;
          end
 
          if(yd >= 0)		% move along y
             y = y + sy;
             yd = yd - ax;
          end
 
          if(zd >= 0)		% move along z
             z = z + sz;
             zd = zd - ax;
          end
 
          x  = x  + sx;		% move along x
          yd = yd + ay;
          zd = zd + az;
       end
    elseif(ay>=max(ax,az))		% y dominant
       xd = ax - ay/2;
       zd = az - ay/2;
 
       while(1)
          X(idx) = x;
          Y(idx) = y;
          Z(idx) = z;
          idx = idx + 1;
 
          if(y == y2)		% end
             break;
          end
 
          if(xd >= 0)		% move along x
             x = x + sx;
             xd = xd - ay;
          end
 
          if(zd >= 0)		% move along z
             z = z + sz;
             zd = zd - ay;
          end
 
          y  = y  + sy;		% move along y
          xd = xd + ax;
          zd = zd + az;
       end
    elseif(az>=max(ax,ay))		% z dominant
       xd = ax - az/2;
       yd = ay - az/2;
 
       while(1)
          X(idx) = x;
          Y(idx) = y;
          Z(idx) = z;
          idx = idx + 1;
 
          if(z == z2)		% end
             break;
          end
 
          if(xd >= 0)		% move along x
             x = x + sx;
             xd = xd - az;
          end
 
          if(yd >= 0)		% move along y
             y = y + sy;
             yd = yd - az;
          end
 
          z  = z  + sz;		% move along z
          xd = xd + ax;
          yd = yd + ay;
       end
    end
 
    if precision ~= 0
       X = X/(10^precision);
       Y = Y/(10^precision);
       Z = Z/(10^precision);
    end
 
    return;					% bresenham_line3d