Browse code

searchlight ready. missing nifti-image-write support. added timing to FBS

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@167 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on30/03/2009 17:54:25
Showing1 changed files
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
+