1.1 diff -r 5edfbd3e7a46 -r 1204ebf9340d PTdecode/CImg-1.3.0/plugins/integral_line.h 1.2 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.3 +++ b/PTdecode/CImg-1.3.0/plugins/integral_line.h Mon Aug 03 14:09:20 2009 +0100 1.4 @@ -0,0 +1,563 @@ 1.5 +/* 1.6 + # 1.7 + # File : integral_line.h 1.8 + # ( C++ header file - CImg plug-in ) 1.9 + # 1.10 + # Description : This CImg plug-in defines function to track integral lines. 1.11 + # This file is a part of the CImg Library project. 1.12 + # ( http://cimg.sourceforge.net ) 1.13 + # 1.14 + # Copyright : David Tschumperle 1.15 + # ( http://www.greyc.ensicaen.fr/~dtschump/ ) 1.16 + # 1.17 + # License : CeCILL v2.0 1.18 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.19 + # 1.20 + # This software is governed by the CeCILL license under French law and 1.21 + # abiding by the rules of distribution of free software. You can use, 1.22 + # modify and/ or redistribute the software under the terms of the CeCILL 1.23 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.24 + # "http://www.cecill.info". 1.25 + # 1.26 + # As a counterpart to the access to the source code and rights to copy, 1.27 + # modify and redistribute granted by the license, users are provided only 1.28 + # with a limited warranty and the software's author, the holder of the 1.29 + # economic rights, and the successive licensors have only limited 1.30 + # liability. 1.31 + # 1.32 + # In this respect, the user's attention is drawn to the risks associated 1.33 + # with loading, using, modifying and/or developing or reproducing the 1.34 + # software by the user in light of its specific status of free software, 1.35 + # that may mean that it is complicated to manipulate, and that also 1.36 + # therefore means that it is reserved for developers and experienced 1.37 + # professionals having in-depth computer knowledge. Users are therefore 1.38 + # encouraged to load and test the software's suitability as regards their 1.39 + # requirements in conditions enabling the security of their systems and/or 1.40 + # data to be ensured and, more generally, to use and operate it in the 1.41 + # same conditions as regards security. 1.42 + # 1.43 + # The fact that you are presently reading this means that you have had 1.44 + # knowledge of the CeCILL license and that you accept its terms. 1.45 + # 1.46 +*/ 1.47 + 1.48 +#ifndef cimg_plugin_integral_line 1.49 +#define cimg_plugin_integral_line 1.50 + 1.51 +#define pcimg_valign2d(i,j) \ 1.52 + { restype &u = W(i,j,0,0), &v = W(i,j,0,1); \ 1.53 + if (u*curru + v*currv<0) { u=-u; v=-v; }} 1.54 +#define pcimg_valign3d(i,j,k) \ 1.55 + { restype &u = W(i,j,k,0), &v = W(i,j,k,1), &w = W(i,j,k,2); \ 1.56 + if (u*curru + v*currv + w*currw<0) { u=-u; v=-v; w=-w; }} 1.57 + 1.58 +CImgList<typename cimg::superset<float,T>::type> 1.59 +get_integral_line(const float x, const float y, const float z=0, 1.60 + const float L=100, const float dl=0.5f, const unsigned int interpolation=3, 1.61 + const bool orientations_only=false) const { 1.62 + 1.63 + typedef typename cimg::superset<float,T>::type restype; 1.64 + CImgList<restype> tracking; 1.65 + CImg<restype> W = (*this)*dl; 1.66 + 1.67 + const unsigned int 1.68 + dx1 = width-1, 1.69 + dy1 = height-1; 1.70 + const float 1.71 + L2 = L/2, 1.72 + cu = (float)(dl*W((int)x,(int)y,(int)z,0)), 1.73 + cv = (float)(dl*W((int)x,(int)y,(int)z,1)); 1.74 + float 1.75 + pu = cu, 1.76 + pv = cv, 1.77 + X = x, 1.78 + Y = y; 1.79 + 1.80 + // 3D integral lines 1.81 + //------------------- 1.82 + switch (W.dimv()) { 1.83 + 1.84 + case 3: { 1.85 + const unsigned int 1.86 + dz1 = depth-1; 1.87 + const float 1.88 + cw = (float)(dl*W((int)x,(int)y,(int)z,2)); 1.89 + float 1.90 + pw = cw, 1.91 + Z = z; 1.92 + 1.93 + switch (interpolation) { 1.94 + case 0: { // Nearest neighbor 1.95 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.96 + tracking.insert(CImg<restype>::vector(X,Y,Z)); 1.97 + const int 1.98 + cx = (int)(X+0.5f), 1.99 + cy = (int)(Y+0.5f), 1.100 + cz = (int)(Z+0.5f); 1.101 + float 1.102 + u = (float)(dl*W(cx,cy,cz,0)), 1.103 + v = (float)(dl*W(cx,cy,cz,1)), 1.104 + w = (float)(dl*W(cx,cy,cz,2)); 1.105 + if (orientations_only && (pu*u + pv*v + pw*w)<0) { u=-u; v=-v; w=-w; } 1.106 + X+=(pu=u); Y+=(pv=v); Z+=(pw=w); 1.107 + } 1.108 + pu = cu; 1.109 + pv = cv; 1.110 + pw = cw; 1.111 + X = x; 1.112 + Y = y; 1.113 + Z = z; 1.114 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.115 + const int 1.116 + cx = (int)(X+0.5f), 1.117 + cy = (int)(Y+0.5f), 1.118 + cz = (int)(Z+0.5f); 1.119 + float 1.120 + u = (float)(dl*W(cx,cy,cz,0)), 1.121 + v = (float)(dl*W(cx,cy,cz,1)), 1.122 + w = (float)(dl*W(cx,cy,cz,2)); 1.123 + if (orientations_only && (pu*u + pv*v + pw*w)<0) { u=-u; v=-v; w=-w; } 1.124 + X-=(pu=u); Y-=(pv=v); Z-=(pw=w); 1.125 + tracking.insert(CImg<restype>::vector(X,Y,Z),0); 1.126 + } 1.127 + } break; 1.128 + 1.129 + case 1: { // Linear 1.130 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.131 + tracking.insert(CImg<restype>::vector(X,Y,Z)); 1.132 + const int 1.133 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.134 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1, 1.135 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1; 1.136 + if (orientations_only) { 1.137 + const float 1.138 + curru = (float)W(cx,cy,cz,0), 1.139 + currv = (float)W(cx,cy,cz,1), 1.140 + currw = (float)W(cx,cy,cz,2); 1.141 + pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz); 1.142 + pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz); 1.143 + pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz); 1.144 + pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz); 1.145 + pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz); 1.146 + pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz); 1.147 + pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz); 1.148 + pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz); 1.149 + pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz); 1.150 + } 1.151 + float 1.152 + u = (float)(dl*W._linear_atXYZ(X,Y,Z,0)), 1.153 + v = (float)(dl*W._linear_atXYZ(X,Y,Z,1)), 1.154 + w = (float)(dl*W._linear_atXYZ(X,Y,Z,2)); 1.155 + if (orientations_only && (pu*u + pv*v + pw*w)<0) { u=-u; v=-v; w=-w; } 1.156 + X+=(pu=u); Y+=(pv=v); Z+=(pw=w); 1.157 + } 1.158 + pu = cu; 1.159 + pv = cv; 1.160 + pw = cw; 1.161 + X = x; 1.162 + Y = y; 1.163 + Z = z; 1.164 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.165 + const int 1.166 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.167 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1, 1.168 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1; 1.169 + if (orientations_only) { 1.170 + const float 1.171 + curru = (float)W(cx,cy,cz,0), 1.172 + currv = (float)W(cx,cy,cz,1), 1.173 + currw = (float)W(cx,cy,cz,2); 1.174 + pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz); 1.175 + pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz); 1.176 + pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz); 1.177 + pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz); 1.178 + pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz); 1.179 + pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz); 1.180 + pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz); 1.181 + pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz); 1.182 + pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz); 1.183 + } 1.184 + float 1.185 + u = (float)(dl*W._linear_atXYZ(X,Y,Z,0)), 1.186 + v = (float)(dl*W._linear_atXYZ(X,Y,Z,1)), 1.187 + w = (float)(dl*W._linear_atXYZ(X,Y,Z,2)); 1.188 + if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; } 1.189 + X-=(pu=u); Y-=(pv=v); Z-=(pw=w); 1.190 + tracking.insert(CImg<restype>::vector(X,Y,Z),0); 1.191 + } 1.192 + 1.193 + } break; 1.194 + 1.195 + case 2: { // 2nd order Runge Kutta 1.196 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.197 + tracking.insert(CImg<restype>::vector(X,Y,Z)); 1.198 + const int 1.199 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.200 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1, 1.201 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1; 1.202 + if (orientations_only) { 1.203 + const float 1.204 + curru = (float)W(cx,cy,cz,0), 1.205 + currv = (float)W(cx,cy,cz,1), 1.206 + currw = (float)W(cx,cy,cz,2); 1.207 + pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz); 1.208 + pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz); 1.209 + pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz); 1.210 + pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz); 1.211 + pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz); 1.212 + pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz); 1.213 + pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz); 1.214 + pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz); 1.215 + pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz); 1.216 + } 1.217 + const float 1.218 + u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)), 1.219 + v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)), 1.220 + w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2)); 1.221 + float 1.222 + u = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)), 1.223 + v = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)), 1.224 + w = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2)); 1.225 + if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; } 1.226 + X+=(pu=u); Y+=(pv=v); Z+=(pw=w); 1.227 + } 1.228 + pu = cu; 1.229 + pv = cv; 1.230 + pw = cw; 1.231 + X = x; 1.232 + Y = y; 1.233 + Z = z; 1.234 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.235 + const int 1.236 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.237 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1, 1.238 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1; 1.239 + if (orientations_only) { 1.240 + const float 1.241 + curru = (float)W(cx,cy,cz,0), 1.242 + currv = (float)W(cx,cy,cz,1), 1.243 + currw = (float)W(cx,cy,cz,2); 1.244 + pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz); 1.245 + pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz); 1.246 + pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz); 1.247 + pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz); 1.248 + pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz); 1.249 + pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz); 1.250 + pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz); 1.251 + pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz); 1.252 + pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz); 1.253 + } 1.254 + const float 1.255 + u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)), 1.256 + v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)), 1.257 + w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2)); 1.258 + float 1.259 + u = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)), 1.260 + v = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)), 1.261 + w = (float)(dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2)); 1.262 + if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; } 1.263 + X-=(pu=u); Y-=(pv=v); Z-=(pw=w); 1.264 + tracking.insert(CImg<restype>::vector(X,Y,Z),0); 1.265 + } 1.266 + } break; 1.267 + 1.268 + case 3: { // 4nd order Runge Kutta 1.269 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.270 + tracking.insert(CImg<restype>::vector(X,Y,Z)); 1.271 + const int 1.272 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.273 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1, 1.274 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1; 1.275 + if (orientations_only) { 1.276 + const float 1.277 + curru = (float)W(cx,cy,cz,0), 1.278 + currv = (float)W(cx,cy,cz,1), 1.279 + currw = (float)W(cx,cy,cz,2); 1.280 + pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz); 1.281 + pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz); 1.282 + pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz); 1.283 + pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz); 1.284 + pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz); 1.285 + pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz); 1.286 + pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz); 1.287 + pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz); 1.288 + pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz); 1.289 + } 1.290 + const float 1.291 + u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)), 1.292 + v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)), 1.293 + w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2)), 1.294 + u1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)), 1.295 + v1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)), 1.296 + w1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2)), 1.297 + u2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,0)), 1.298 + v2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,1)), 1.299 + w2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,2)), 1.300 + u3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,0)), 1.301 + v3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,1)), 1.302 + w3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,2)); 1.303 + float 1.304 + u = u0/6 + u1/3 + u2/3 + u3/6, 1.305 + v = v0/6 + v1/3 + v2/3 + v3/6, 1.306 + w = w0/6 + w1/3 + w2/3 + w3/6; 1.307 + if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; } 1.308 + X+=(pu=u); Y+=(pv=v); Z+=(pw=w); 1.309 + } 1.310 + pu = cu; 1.311 + pv = cv; 1.312 + pw = cw; 1.313 + X = x; 1.314 + Y = y; 1.315 + Z = z; 1.316 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) { 1.317 + const int 1.318 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.319 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1, 1.320 + cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>(int)dz1)?(int)dz1:cz+1; 1.321 + if (orientations_only) { 1.322 + const float 1.323 + curru = (float)W(cx,cy,cz,0), 1.324 + currv = (float)W(cx,cy,cz,1), 1.325 + currw = (float)W(cx,cy,cz,2); 1.326 + pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz); 1.327 + pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz); 1.328 + pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz); 1.329 + pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz); 1.330 + pcimg_valign3d(px,cy,cz); pcimg_valign3d(nx,cy,cz); 1.331 + pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz); 1.332 + pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz); 1.333 + pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz); 1.334 + pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz); 1.335 + } 1.336 + const float 1.337 + u0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,0)), 1.338 + v0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,1)), 1.339 + w0 = (float)(0.5f*dl*W._linear_atXYZ(X,Y,Z,2)), 1.340 + u1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,0)), 1.341 + v1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,1)), 1.342 + w1 = (float)(0.5f*dl*W._linear_atXYZ(X+u0,Y+v0,Z+w0,2)), 1.343 + u2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,0)), 1.344 + v2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,1)), 1.345 + w2 = (float)(0.5f*dl*W._linear_atXYZ(X+u1,Y+v1,Z+w1,2)), 1.346 + u3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,0)), 1.347 + v3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,1)), 1.348 + w3 = (float)(0.5f*dl*W._linear_atXYZ(X+u2,Y+v2,Z+w2,2)); 1.349 + float 1.350 + u = u0/6 + u1/3 + u2/3 + u3/6, 1.351 + v = v0/6 + v1/3 + v2/3 + v3/6, 1.352 + w = w0/6 + w1/3 + w2/3 + w3/6; 1.353 + if (orientations_only && (pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; } 1.354 + X-=(pu=u); Y-=(pv=v); Z-=(pw=w); 1.355 + tracking.insert(CImg<restype>::vector(X,Y,Z),0); 1.356 + } 1.357 + } break; 1.358 + } 1.359 + 1.360 + } break; 1.361 + 1.362 + // 2D integral lines 1.363 + //------------------- 1.364 + case 2: { 1.365 + 1.366 + switch (interpolation) { 1.367 + case 0: { // Nearest neighbor 1.368 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.369 + tracking.insert(CImg<restype>::vector(X,Y)); 1.370 + const int 1.371 + cx = (int)(X+0.5f), 1.372 + cy = (int)(Y+0.5f); 1.373 + float 1.374 + u = (float)(dl*W(cx,cy,0,0)), 1.375 + v = (float)(dl*W(cx,cy,0,1)); 1.376 + if (orientations_only && (pu*u + pv*v)<0) { u=-u; v=-v; } 1.377 + X+=(pu=u); Y+=(pv=v); 1.378 + } 1.379 + pu = cu; 1.380 + pv = cv; 1.381 + X = x; 1.382 + Y = y; 1.383 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.384 + const int 1.385 + cx = (int)(X+0.5f), 1.386 + cy = (int)(Y+0.5f); 1.387 + float 1.388 + u = (float)(dl*W(cx,cy,0,0)), 1.389 + v = (float)(dl*W(cx,cy,0,1)); 1.390 + if (orientations_only && (pu*u + pv*v)<0) { u=-u; v=-v; } 1.391 + X-=(pu=u); Y-=(pv=v); 1.392 + tracking.insert(CImg<restype>::vector(X,Y),0); 1.393 + } 1.394 + } break; 1.395 + 1.396 + case 1: { // Linear 1.397 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.398 + tracking.insert(CImg<restype>::vector(X,Y)); 1.399 + const int 1.400 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.401 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1; 1.402 + if (orientations_only) { 1.403 + const float 1.404 + curru = (float)W(cx,cy,0,0), 1.405 + currv = (float)W(cx,cy,0,1); 1.406 + pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py); 1.407 + pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy); 1.408 + pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny); 1.409 + } 1.410 + float 1.411 + u = (float)(dl*W._linear_atXY(X,Y,0,0)), 1.412 + v = (float)(dl*W._linear_atXY(X,Y,0,1)); 1.413 + if (orientations_only && (pu*u + pv*v)<0) { u=-u; v=-v; } 1.414 + X+=(pu=u); Y+=(pv=v); 1.415 + } 1.416 + pu = cu; 1.417 + pv = cv; 1.418 + X = x; 1.419 + Y = y; 1.420 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.421 + const int 1.422 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.423 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1; 1.424 + if (orientations_only) { 1.425 + const float 1.426 + curru = (float)W(cx,cy,0,0), 1.427 + currv = (float)W(cx,cy,0,1); 1.428 + pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py); 1.429 + pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy); 1.430 + pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny); 1.431 + } 1.432 + float 1.433 + u = (float)(dl*W._linear_atXY(X,Y,0,0)), 1.434 + v = (float)(dl*W._linear_atXY(X,Y,0,1)); 1.435 + if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; } 1.436 + X-=(pu=u); Y-=(pv=v); 1.437 + tracking.insert(CImg<restype>::vector(X,Y),0); 1.438 + } 1.439 + } break; 1.440 + 1.441 + case 2: { // 2nd order Runge Kutta 1.442 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.443 + tracking.insert(CImg<restype>::vector(X,Y)); 1.444 + const int 1.445 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.446 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1; 1.447 + if (orientations_only) { 1.448 + const float 1.449 + curru = (float)W(cx,cy,0,0), 1.450 + currv = (float)W(cx,cy,0,1); 1.451 + pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py); 1.452 + pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy); 1.453 + pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny); 1.454 + } 1.455 + const float 1.456 + u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)), 1.457 + v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1)); 1.458 + float 1.459 + u = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,0)), 1.460 + v = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,1)); 1.461 + if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; } 1.462 + X+=(pu=u); Y+=(pv=v); 1.463 + } 1.464 + pu = cu; 1.465 + pv = cv; 1.466 + X = x; 1.467 + Y = y; 1.468 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.469 + const int 1.470 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.471 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1; 1.472 + if (orientations_only) { 1.473 + const float 1.474 + curru = (float)W(cx,cy,0,0), 1.475 + currv = (float)W(cx,cy,0,1); 1.476 + pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py); 1.477 + pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy); 1.478 + pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny); 1.479 + } 1.480 + const float 1.481 + u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)), 1.482 + v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1)); 1.483 + float 1.484 + u = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,0)), 1.485 + v = (float)(dl*W._linear_atXY(X+u0,Y+v0,0,1)); 1.486 + if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; } 1.487 + X-=(pu=u); Y-=(pv=v); 1.488 + tracking.insert(CImg<restype>::vector(X,Y),0); 1.489 + } 1.490 + } break; 1.491 + 1.492 + case 3: { // 4nd order Runge Kutta 1.493 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.494 + tracking.insert(CImg<restype>::vector(X,Y)); 1.495 + const int 1.496 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.497 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1; 1.498 + if (orientations_only) { 1.499 + const float 1.500 + curru = (float)W(cx,cy,0,0), 1.501 + currv = (float)W(cx,cy,0,1); 1.502 + pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py); 1.503 + pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy); 1.504 + pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny); 1.505 + } 1.506 + const float 1.507 + u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)), 1.508 + v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1)), 1.509 + u1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,0)), 1.510 + v1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,1)), 1.511 + u2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,0)), 1.512 + v2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,1)), 1.513 + u3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,0)), 1.514 + v3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,1)); 1.515 + float 1.516 + u = u0/6 + u1/3 + u2/3 + u3/6, 1.517 + v = v0/6 + v1/3 + v2/3 + v3/6; 1.518 + if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; } 1.519 + X+=(pu=u); Y+=(pv=v); 1.520 + } 1.521 + pu = cu; 1.522 + pv = cv; 1.523 + X = x; 1.524 + Y = y; 1.525 + for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) { 1.526 + const int 1.527 + cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>(int)dx1)?(int)dx1:cx+1, 1.528 + cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>(int)dy1)?(int)dy1:cy+1; 1.529 + if (orientations_only) { 1.530 + const float 1.531 + curru = (float)W(cx,cy,0,0), 1.532 + currv = (float)W(cx,cy,0,1); 1.533 + pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py); 1.534 + pcimg_valign2d(px,cy); pcimg_valign2d(nx,cy); 1.535 + pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny); 1.536 + } 1.537 + const float 1.538 + u0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,0)), 1.539 + v0 = (float)(0.5f*dl*W._linear_atXY(X,Y,0,1)), 1.540 + u1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,0)), 1.541 + v1 = (float)(0.5f*dl*W._linear_atXY(X+u0,Y+v0,0,1)), 1.542 + u2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,0)), 1.543 + v2 = (float)(0.5f*dl*W._linear_atXY(X+u1,Y+v1,0,1)), 1.544 + u3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,0)), 1.545 + v3 = (float)(0.5f*dl*W._linear_atXY(X+u2,Y+v2,0,1)); 1.546 + float 1.547 + u = u0/6 + u1/3 + u2/3 + u3/6, 1.548 + v = v0/6 + v1/3 + v2/3 + v3/6; 1.549 + if (orientations_only && (pu*u+pv*v)<0) { u=-u; v=-v; } 1.550 + X-=(pu=u); Y-=(pv=v); 1.551 + tracking.insert(CImg<restype>::vector(X,Y),0); 1.552 + } 1.553 + } break; 1.554 + } 1.555 + 1.556 + } break; 1.557 + 1.558 + default: 1.559 + throw CImgInstanceException("CImg<%s>::get_integral_line() : Instance image must have dimv()=2 or 3 (current is %u).", 1.560 + pixel_type(),dim); 1.561 + break; 1.562 + } 1.563 + 1.564 + return tracking; 1.565 +} 1.566 + 1.567 +#endif