Mon, 03 Aug 2009 14:09:20 +0100
added P-touch decoder source
1 /*
2 #
3 # File : integral_line.h
4 # ( C++ header file - CImg plug-in )
5 #
6 # Description : This CImg plug-in defines function to track integral lines.
7 # This file is a part of the CImg Library project.
8 # ( http://cimg.sourceforge.net )
9 #
10 # Copyright : David Tschumperle
11 # ( http://www.greyc.ensicaen.fr/~dtschump/ )
12 #
13 # License : CeCILL v2.0
14 # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
15 #
16 # This software is governed by the CeCILL license under French law and
17 # abiding by the rules of distribution of free software. You can use,
18 # modify and/ or redistribute the software under the terms of the CeCILL
19 # license as circulated by CEA, CNRS and INRIA at the following URL
20 # "http://www.cecill.info".
21 #
22 # As a counterpart to the access to the source code and rights to copy,
23 # modify and redistribute granted by the license, users are provided only
24 # with a limited warranty and the software's author, the holder of the
25 # economic rights, and the successive licensors have only limited
26 # liability.
27 #
28 # In this respect, the user's attention is drawn to the risks associated
29 # with loading, using, modifying and/or developing or reproducing the
30 # software by the user in light of its specific status of free software,
31 # that may mean that it is complicated to manipulate, and that also
32 # therefore means that it is reserved for developers and experienced
33 # professionals having in-depth computer knowledge. Users are therefore
34 # encouraged to load and test the software's suitability as regards their
35 # requirements in conditions enabling the security of their systems and/or
36 # data to be ensured and, more generally, to use and operate it in the
37 # same conditions as regards security.
38 #
39 # The fact that you are presently reading this means that you have had
40 # knowledge of the CeCILL license and that you accept its terms.
41 #
42 */
44 #ifndef cimg_plugin_integral_line
45 #define cimg_plugin_integral_line
47 #define pcimg_valign2d(i,j) \
48 { restype &u = W(i,j,0,0), &v = W(i,j,0,1); \
49 if (u*curru + v*currv<0) { u=-u; v=-v; }}
50 #define pcimg_valign3d(i,j,k) \
51 { restype &u = W(i,j,k,0), &v = W(i,j,k,1), &w = W(i,j,k,2); \
52 if (u*curru + v*currv + w*currw<0) { u=-u; v=-v; w=-w; }}
54 CImgList<typename cimg::superset<float,T>::type>
55 get_integral_line(const float x, const float y, const float z=0,
56 const float L=100, const float dl=0.5f, const unsigned int interpolation=3,
57 const bool orientations_only=false) const {
59 typedef typename cimg::superset<float,T>::type restype;
60 CImgList<restype> tracking;
61 CImg<restype> W = (*this)*dl;
63 const unsigned int
64 dx1 = width-1,
65 dy1 = height-1;
66 const float
67 L2 = L/2,
68 cu = (float)(dl*W((int)x,(int)y,(int)z,0)),
69 cv = (float)(dl*W((int)x,(int)y,(int)z,1));
70 float
71 pu = cu,
72 pv = cv,
73 X = x,
74 Y = y;
76 // 3D integral lines
77 //-------------------
78 switch (W.dimv()) {
80 case 3: {
81 const unsigned int
82 dz1 = depth-1;
83 const float
84 cw = (float)(dl*W((int)x,(int)y,(int)z,2));
85 float
86 pw = cw,
87 Z = z;
89 switch (interpolation) {
90 case 0: { // Nearest neighbor
91 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
92 tracking.insert(CImg<restype>::vector(X,Y,Z));
93 const int
94 cx = (int)(X+0.5f),
95 cy = (int)(Y+0.5f),
96 cz = (int)(Z+0.5f);
97 float
98 u = (float)(dl*W(cx,cy,cz,0)),
99 v = (float)(dl*W(cx,cy,cz,1)),
100 w = (float)(dl*W(cx,cy,cz,2));
101 if (orientations_only && (pu*u + pv*v + pw*w)<0) { u=-u; v=-v; w=-w; }
102 X+=(pu=u); Y+=(pv=v); Z+=(pw=w);
103 }
104 pu = cu;
105 pv = cv;
106 pw = cw;
107 X = x;
108 Y = y;
109 Z = z;
110 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
111 const int
112 cx = (int)(X+0.5f),
113 cy = (int)(Y+0.5f),
114 cz = (int)(Z+0.5f);
115 float
116 u = (float)(dl*W(cx,cy,cz,0)),
117 v = (float)(dl*W(cx,cy,cz,1)),
118 w = (float)(dl*W(cx,cy,cz,2));
119 if (orientations_only && (pu*u + pv*v + pw*w)<0) { u=-u; v=-v; w=-w; }
120 X-=(pu=u); Y-=(pv=v); Z-=(pw=w);
121 tracking.insert(CImg<restype>::vector(X,Y,Z),0);
122 }
123 } break;
125 case 1: { // Linear
126 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
127 tracking.insert(CImg<restype>::vector(X,Y,Z));
128 const int
129 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
130 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1,
131 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1;
132 if (orientations_only) {
133 const float
134 curru = (float)W(cx,cy,cz,0),
135 currv = (float)W(cx,cy,cz,1),
136 currw = (float)W(cx,cy,cz,2);
137 pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);
138 pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);
139 pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);
140 pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);
141 pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz);
142 pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);
143 pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);
144 pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);
145 pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);
146 }
147 float
148 u = (float)(dl*W._linear_atXYZ(X,Y,Z,0)),
149 v = (float)(dl*W._linear_atXYZ(X,Y,Z,1)),
150 w = (float)(dl*W._linear_atXYZ(X,Y,Z,2));
151 if (orientations_only && (pu*u + pv*v + pw*w)<0) { u=-u; v=-v; w=-w; }
152 X+=(pu=u); Y+=(pv=v); Z+=(pw=w);
153 }
154 pu = cu;
155 pv = cv;
156 pw = cw;
157 X = x;
158 Y = y;
159 Z = z;
160 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
161 const int
162 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
163 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1,
164 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1;
165 if (orientations_only) {
166 const float
167 curru = (float)W(cx,cy,cz,0),
168 currv = (float)W(cx,cy,cz,1),
169 currw = (float)W(cx,cy,cz,2);
170 pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);
171 pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);
172 pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);
173 pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);
174 pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz);
175 pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);
176 pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);
177 pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);
178 pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);
179 }
180 float
181 u = (float)(dl*W._linear_atXYZ(X,Y,Z,0)),
182 v = (float)(dl*W._linear_atXYZ(X,Y,Z,1)),
183 w = (float)(dl*W._linear_atXYZ(X,Y,Z,2));
184 if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }
185 X-=(pu=u); Y-=(pv=v); Z-=(pw=w);
186 tracking.insert(CImg<restype>::vector(X,Y,Z),0);
187 }
189 } break;
191 case 2: { // 2nd order Runge Kutta
192 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
193 tracking.insert(CImg<restype>::vector(X,Y,Z));
194 const int
195 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
196 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1,
197 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1;
198 if (orientations_only) {
199 const float
200 curru = (float)W(cx,cy,cz,0),
201 currv = (float)W(cx,cy,cz,1),
202 currw = (float)W(cx,cy,cz,2);
203 pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);
204 pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);
205 pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);
206 pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);
207 pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz);
208 pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);
209 pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);
210 pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);
211 pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);
212 }
213 const float
214 u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)),
215 v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)),
216 w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2));
217 float
218 u = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)),
219 v = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)),
220 w = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2));
221 if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }
222 X+=(pu=u); Y+=(pv=v); Z+=(pw=w);
223 }
224 pu = cu;
225 pv = cv;
226 pw = cw;
227 X = x;
228 Y = y;
229 Z = z;
230 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
231 const int
232 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
233 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1,
234 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1;
235 if (orientations_only) {
236 const float
237 curru = (float)W(cx,cy,cz,0),
238 currv = (float)W(cx,cy,cz,1),
239 currw = (float)W(cx,cy,cz,2);
240 pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);
241 pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);
242 pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);
243 pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);
244 pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz);
245 pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);
246 pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);
247 pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);
248 pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);
249 }
250 const float
251 u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)),
252 v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)),
253 w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2));
254 float
255 u = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)),
256 v = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)),
257 w = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2));
258 if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }
259 X-=(pu=u); Y-=(pv=v); Z-=(pw=w);
260 tracking.insert(CImg<restype>::vector(X,Y,Z),0);
261 }
262 } break;
264 case 3: { // 4nd order Runge Kutta
265 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
266 tracking.insert(CImg<restype>::vector(X,Y,Z));
267 const int
268 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
269 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1,
270 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1;
271 if (orientations_only) {
272 const float
273 curru = (float)W(cx,cy,cz,0),
274 currv = (float)W(cx,cy,cz,1),
275 currw = (float)W(cx,cy,cz,2);
276 pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);
277 pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);
278 pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);
279 pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);
280 pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz);
281 pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);
282 pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);
283 pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);
284 pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);
285 }
286 const float
287 u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)),
288 v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)),
289 w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2)),
290 u1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)),
291 v1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)),
292 w1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2)),
293 u2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,0)),
294 v2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,1)),
295 w2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,2)),
296 u3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,0)),
297 v3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,1)),
298 w3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,2));
299 float
300 u = u0/6 + u1/3 + u2/3 + u3/6,
301 v = v0/6 + v1/3 + v2/3 + v3/6,
302 w = w0/6 + w1/3 + w2/3 + w3/6;
303 if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }
304 X+=(pu=u); Y+=(pv=v); Z+=(pw=w);
305 }
306 pu = cu;
307 pv = cv;
308 pw = cw;
309 X = x;
310 Y = y;
311 Z = z;
312 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {
313 const int
314 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
315 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1,
316 cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1;
317 if (orientations_only) {
318 const float
319 curru = (float)W(cx,cy,cz,0),
320 currv = (float)W(cx,cy,cz,1),
321 currw = (float)W(cx,cy,cz,2);
322 pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);
323 pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);
324 pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);
325 pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);
326 pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz);
327 pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);
328 pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);
329 pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);
330 pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);
331 }
332 const float
333 u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)),
334 v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)),
335 w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2)),
336 u1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)),
337 v1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)),
338 w1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2)),
339 u2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,0)),
340 v2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,1)),
341 w2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,2)),
342 u3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,0)),
343 v3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,1)),
344 w3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,2));
345 float
346 u = u0/6 + u1/3 + u2/3 + u3/6,
347 v = v0/6 + v1/3 + v2/3 + v3/6,
348 w = w0/6 + w1/3 + w2/3 + w3/6;
349 if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }
350 X-=(pu=u); Y-=(pv=v); Z-=(pw=w);
351 tracking.insert(CImg<restype>::vector(X,Y,Z),0);
352 }
353 } break;
354 }
356 } break;
358 // 2D integral lines
359 //-------------------
360 case 2: {
362 switch (interpolation) {
363 case 0: { // Nearest neighbor
364 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
365 tracking.insert(CImg<restype>::vector(X,Y));
366 const int
367 cx = (int)(X+0.5f),
368 cy = (int)(Y+0.5f);
369 float
370 u = (float)(dl*W(cx,cy,0,0)),
371 v = (float)(dl*W(cx,cy,0,1));
372 if (orientations_only && (pu*u + pv*v)<0) { u=-u; v=-v; }
373 X+=(pu=u); Y+=(pv=v);
374 }
375 pu = cu;
376 pv = cv;
377 X = x;
378 Y = y;
379 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
380 const int
381 cx = (int)(X+0.5f),
382 cy = (int)(Y+0.5f);
383 float
384 u = (float)(dl*W(cx,cy,0,0)),
385 v = (float)(dl*W(cx,cy,0,1));
386 if (orientations_only && (pu*u + pv*v)<0) { u=-u; v=-v; }
387 X-=(pu=u); Y-=(pv=v);
388 tracking.insert(CImg<restype>::vector(X,Y),0);
389 }
390 } break;
392 case 1: { // Linear
393 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
394 tracking.insert(CImg<restype>::vector(X,Y));
395 const int
396 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
397 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1;
398 if (orientations_only) {
399 const float
400 curru = (float)W(cx,cy,0,0),
401 currv = (float)W(cx,cy,0,1);
402 pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);
403 pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy);
404 pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);
405 }
406 float
407 u = (float)(dl*W._linear_atXY(X,Y,0,0)),
408 v = (float)(dl*W._linear_atXY(X,Y,0,1));
409 if (orientations_only && (pu*u + pv*v)<0) { u=-u; v=-v; }
410 X+=(pu=u); Y+=(pv=v);
411 }
412 pu = cu;
413 pv = cv;
414 X = x;
415 Y = y;
416 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
417 const int
418 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
419 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1;
420 if (orientations_only) {
421 const float
422 curru = (float)W(cx,cy,0,0),
423 currv = (float)W(cx,cy,0,1);
424 pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);
425 pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy);
426 pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);
427 }
428 float
429 u = (float)(dl*W._linear_atXY(X,Y,0,0)),
430 v = (float)(dl*W._linear_atXY(X,Y,0,1));
431 if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; }
432 X-=(pu=u); Y-=(pv=v);
433 tracking.insert(CImg<restype>::vector(X,Y),0);
434 }
435 } break;
437 case 2: { // 2nd order Runge Kutta
438 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
439 tracking.insert(CImg<restype>::vector(X,Y));
440 const int
441 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
442 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1;
443 if (orientations_only) {
444 const float
445 curru = (float)W(cx,cy,0,0),
446 currv = (float)W(cx,cy,0,1);
447 pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);
448 pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy);
449 pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);
450 }
451 const float
452 u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)),
453 v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1));
454 float
455 u = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,0)),
456 v = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,1));
457 if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; }
458 X+=(pu=u); Y+=(pv=v);
459 }
460 pu = cu;
461 pv = cv;
462 X = x;
463 Y = y;
464 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
465 const int
466 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
467 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1;
468 if (orientations_only) {
469 const float
470 curru = (float)W(cx,cy,0,0),
471 currv = (float)W(cx,cy,0,1);
472 pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);
473 pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy);
474 pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);
475 }
476 const float
477 u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)),
478 v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1));
479 float
480 u = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,0)),
481 v = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,1));
482 if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; }
483 X-=(pu=u); Y-=(pv=v);
484 tracking.insert(CImg<restype>::vector(X,Y),0);
485 }
486 } break;
488 case 3: { // 4nd order Runge Kutta
489 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
490 tracking.insert(CImg<restype>::vector(X,Y));
491 const int
492 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
493 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1;
494 if (orientations_only) {
495 const float
496 curru = (float)W(cx,cy,0,0),
497 currv = (float)W(cx,cy,0,1);
498 pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);
499 pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy);
500 pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);
501 }
502 const float
503 u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)),
504 v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1)),
505 u1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,0)),
506 v1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,1)),
507 u2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,0)),
508 v2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,1)),
509 u3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,0)),
510 v3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,1));
511 float
512 u = u0/6 + u1/3 + u2/3 + u3/6,
513 v = v0/6 + v1/3 + v2/3 + v3/6;
514 if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; }
515 X+=(pu=u); Y+=(pv=v);
516 }
517 pu = cu;
518 pv = cv;
519 X = x;
520 Y = y;
521 for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {
522 const int
523 cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1,
524 cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1;
525 if (orientations_only) {
526 const float
527 curru = (float)W(cx,cy,0,0),
528 currv = (float)W(cx,cy,0,1);
529 pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);
530 pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy);
531 pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);
532 }
533 const float
534 u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)),
535 v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1)),
536 u1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,0)),
537 v1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,1)),
538 u2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,0)),
539 v2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,1)),
540 u3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,0)),
541 v3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,1));
542 float
543 u = u0/6 + u1/3 + u2/3 + u3/6,
544 v = v0/6 + v1/3 + v2/3 + v3/6;
545 if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; }
546 X-=(pu=u); Y-=(pv=v);
547 tracking.insert(CImg<restype>::vector(X,Y),0);
548 }
549 } break;
550 }
552 } break;
554 default:
555 throw CImgInstanceException("CImg<%s>::get_integral_line() : Instance image must have dimv()=2 or 3 (current is %u).",
556 pixel_type(),dim);
557 break;
558 }
560 return tracking;
561 }
563 #endif