 Christoph Budziszewski authored on30/03/2009 17:54:25
 1 1 `new file mode 100644` ... ... `@@ -0,0 +1,189 @@` 1 ```+% Generate X Y Z coordinates of a 3D Bresenham's line between ``` 2 ```+% two given points. ``` 3 ```+% ``` 4 ```+% A very useful application of this algorithm can be found in the ``` 5 ```+% implementation of Fischer's Bresenham interpolation method in my ``` 6 ```+% another program that can rotate three dimensional image volume ``` 7 ```+% with an affine matrix: ``` 8 ```+% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21080 ``` 9 ```+% ``` 10 ```+% Usage: [X Y Z] = bresenham_line3d(P1, P2, [precision]); ``` 11 ```+% ``` 12 ```+% P1 - vector for Point1, where P1 = [x1 y1 z1] ``` 13 ```+% ``` 14 ```+% P2 - vector for Point2, where P2 = [x2 y2 z2] ``` 15 ```+% ``` 16 ```+% precision (optional) - Although according to Bresenham's line ``` 17 ```+% algorithm, point coordinates x1 y1 z1 and x2 y2 z2 should ``` 18 ```+% be integer numbers, this program extends its limit to all ``` 19 ```+% real numbers. If any of them are floating numbers, you ``` 20 ```+% should specify how many digits of decimal that you would ``` 21 ```+% like to preserve. Be aware that the length of output X Y ``` 22 ```+% Z coordinates will increase in 10 times for each decimal ``` 23 ```+% digit that you want to preserve. By default, the precision ``` 24 ```+% is 0, which means that they will be rounded to the nearest ``` 25 ```+% integer. ``` 26 ```+% ``` 27 ```+% X - a set of x coordinates on Bresenham's line ``` 28 ```+% ``` 29 ```+% Y - a set of y coordinates on Bresenham's line ``` 30 ```+% ``` 31 ```+% Z - a set of z coordinates on Bresenham's line ``` 32 ```+% ``` 33 ```+% Therefore, all points in XYZ set (i.e. P(i) = [X(i) Y(i) Z(i)]) ``` 34 ```+% will constitute the Bresenham's line between P1 and P1. ``` 35 ```+% ``` 36 ```+% Example: ``` 37 ```+% P1 = [12 37 6]; P2 = [46 3 35]; ``` 38 ```+% [X Y Z] = bresenham_line3d(P1, P2); ``` 39 ```+% figure; plot3(X,Y,Z,'s','markerface','b'); ``` 40 ```+% ``` 41 ```+% This program is ported to MATLAB from: ``` 42 ```+% ``` 43 ```+% B.Pendleton. line3d - 3D Bresenham's (a 3D line drawing algorithm) ``` 44 ```+% ftp://ftp.isc.org/pub/usenet/comp.sources.unix/volume26/line3d, 1992 ``` 45 ```+% ``` 46 ```+% Which is also referenced by: ``` 47 ```+% ``` 48 ```+% Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid ``` 49 ```+% Transformations to Volume Data, WSCG2004 Conference. ``` 50 ```+% http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf ``` 51 ```+% ``` 52 ```+% - Jimmy Shen (jimmy@rotman-baycrest.on.ca) ``` 53 ```+% ``` 54 ```+function [X,Y,Z] = bresenham_line3d(P1, P2, precision) ``` 55 ```+ ``` 56 ```+ if ~exist('precision','var') | isempty(precision) | round(precision) == 0 ``` 57 ```+ precision = 0; ``` 58 ```+ P1 = round(P1); ``` 59 ```+ P2 = round(P2); ``` 60 ```+ else ``` 61 ```+ precision = round(precision); ``` 62 ```+ P1 = round(P1*(10^precision)); ``` 63 ```+ P2 = round(P2*(10^precision)); ``` 64 ```+ end ``` 65 ```+ ``` 66 ```+ d = max(abs(P2-P1)+1); ``` 67 ```+ X = zeros(1, d); ``` 68 ```+ Y = zeros(1, d); ``` 69 ```+ Z = zeros(1, d); ``` 70 ```+ ``` 71 ```+ x1 = P1(1); ``` 72 ```+ y1 = P1(2); ``` 73 ```+ z1 = P1(3); ``` 74 ```+ ``` 75 ```+ x2 = P2(1); ``` 76 ```+ y2 = P2(2); ``` 77 ```+ z2 = P2(3); ``` 78 ```+ ``` 79 ```+ dx = x2 - x1; ``` 80 ```+ dy = y2 - y1; ``` 81 ```+ dz = z2 - z1; ``` 82 ```+ ``` 83 ```+ ax = abs(dx)*2; ``` 84 ```+ ay = abs(dy)*2; ``` 85 ```+ az = abs(dz)*2; ``` 86 ```+ ``` 87 ```+ sx = sign(dx); ``` 88 ```+ sy = sign(dy); ``` 89 ```+ sz = sign(dz); ``` 90 ```+ ``` 91 ```+ x = x1; ``` 92 ```+ y = y1; ``` 93 ```+ z = z1; ``` 94 ```+ idx = 1; ``` 95 ```+ ``` 96 ```+ if(ax>=max(ay,az)) % x dominant ``` 97 ```+ yd = ay - ax/2; ``` 98 ```+ zd = az - ax/2; ``` 99 ```+ ``` 100 ```+ while(1) ``` 101 ```+ X(idx) = x; ``` 102 ```+ Y(idx) = y; ``` 103 ```+ Z(idx) = z; ``` 104 ```+ idx = idx + 1; ``` 105 ```+ ``` 106 ```+ if(x == x2) % end ``` 107 ```+ break; ``` 108 ```+ end ``` 109 ```+ ``` 110 ```+ if(yd >= 0) % move along y ``` 111 ```+ y = y + sy; ``` 112 ```+ yd = yd - ax; ``` 113 ```+ end ``` 114 ```+ ``` 115 ```+ if(zd >= 0) % move along z ``` 116 ```+ z = z + sz; ``` 117 ```+ zd = zd - ax; ``` 118 ```+ end ``` 119 ```+ ``` 120 ```+ x = x + sx; % move along x ``` 121 ```+ yd = yd + ay; ``` 122 ```+ zd = zd + az; ``` 123 ```+ end ``` 124 ```+ elseif(ay>=max(ax,az)) % y dominant ``` 125 ```+ xd = ax - ay/2; ``` 126 ```+ zd = az - ay/2; ``` 127 ```+ ``` 128 ```+ while(1) ``` 129 ```+ X(idx) = x; ``` 130 ```+ Y(idx) = y; ``` 131 ```+ Z(idx) = z; ``` 132 ```+ idx = idx + 1; ``` 133 ```+ ``` 134 ```+ if(y == y2) % end ``` 135 ```+ break; ``` 136 ```+ end ``` 137 ```+ ``` 138 ```+ if(xd >= 0) % move along x ``` 139 ```+ x = x + sx; ``` 140 ```+ xd = xd - ay; ``` 141 ```+ end ``` 142 ```+ ``` 143 ```+ if(zd >= 0) % move along z ``` 144 ```+ z = z + sz; ``` 145 ```+ zd = zd - ay; ``` 146 ```+ end ``` 147 ```+ ``` 148 ```+ y = y + sy; % move along y ``` 149 ```+ xd = xd + ax; ``` 150 ```+ zd = zd + az; ``` 151 ```+ end ``` 152 ```+ elseif(az>=max(ax,ay)) % z dominant ``` 153 ```+ xd = ax - az/2; ``` 154 ```+ yd = ay - az/2; ``` 155 ```+ ``` 156 ```+ while(1) ``` 157 ```+ X(idx) = x; ``` 158 ```+ Y(idx) = y; ``` 159 ```+ Z(idx) = z; ``` 160 ```+ idx = idx + 1; ``` 161 ```+ ``` 162 ```+ if(z == z2) % end ``` 163 ```+ break; ``` 164 ```+ end ``` 165 ```+ ``` 166 ```+ if(xd >= 0) % move along x ``` 167 ```+ x = x + sx; ``` 168 ```+ xd = xd - az; ``` 169 ```+ end ``` 170 ```+ ``` 171 ```+ if(yd >= 0) % move along y ``` 172 ```+ y = y + sy; ``` 173 ```+ yd = yd - az; ``` 174 ```+ end ``` 175 ```+ ``` 176 ```+ z = z + sz; % move along z ``` 177 ```+ xd = xd + ax; ``` 178 ```+ yd = yd + ay; ``` 179 ```+ end ``` 180 ```+ end ``` 181 ```+ ``` 182 ```+ if precision ~= 0 ``` 183 ```+ X = X/(10^precision); ``` 184 ```+ Y = Y/(10^precision); ``` 185 ```+ Z = Z/(10^precision); ``` 186 ```+ end ``` 187 ```+ ``` 188 ```+ return; % bresenham_line3d ``` 189 ```+ ```