Code covered by the BSD License  

Highlights from
Tools for NIfTI and ANALYZE image

image thumbnail

Tools for NIfTI and ANALYZE image

by

Jimmy Shen (view profile)

 

23 Oct 2005 (Updated )

Load, save, make, reslice, view (and edit) both NIfTI and ANALYZE data on any platform

bresenham_line3d(P1, P2, precision)
%  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

Contact us