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