pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);        const float          u0 = (float)(0.5f*dl*W.linear_pix3d(X,Y,Z,0)),          v0 = (float)(0.5f*dl*W.linear_pix3d(X,Y,Z,1)),          w0 = (float)(0.5f*dl*W.linear_pix3d(X,Y,Z,2)),          u1 = (float)(0.5f*dl*W.linear_pix3d(X+u0,Y+v0,Z+w0,0)),          v1 = (float)(0.5f*dl*W.linear_pix3d(X+u0,Y+v0,Z+w0,1)),          w1 = (float)(0.5f*dl*W.linear_pix3d(X+u0,Y+v0,Z+w0,2)),          u2 = (float)(0.5f*dl*W.linear_pix3d(X+u1,Y+v1,Z+w1,0)),          v2 = (float)(0.5f*dl*W.linear_pix3d(X+u1,Y+v1,Z+w1,1)),          w2 = (float)(0.5f*dl*W.linear_pix3d(X+u1,Y+v1,Z+w1,2)),          u3 = (float)(0.5f*dl*W.linear_pix3d(X+u2,Y+v2,Z+w2,0)),          v3 = (float)(0.5f*dl*W.linear_pix3d(X+u2,Y+v2,Z+w2,1)),          w3 = (float)(0.5f*dl*W.linear_pix3d(X+u2,Y+v2,Z+w2,2));        float          u = u0/6 + u1/3 + u2/3 + u3/6,          v = v0/6 + v1/3 + v2/3 + v3/6,          w = w0/6 + w1/3 + w2/3 + w3/6;        if ((pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }        X+=(pu=u); Y+=(pv=v); Z+=(pw=w);      }      pu = cu;      pv = cv;      pw = cw;      X  = x;      Y  = y;      Z  = z;      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1 && Z>=0 && Z<=dz1; l+=dl) {        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1,          cz = (int)Z, pz = (cz-1<0)?0:cz-1, nz = (cz+1>dz1)?dz1:cz+1;        const float          curru = (float)W(cx,cy,cz,0),          currv = (float)W(cx,cy,cz,1),          currw = (float)W(cx,cy,cz,2);        pcimg_valign3d(px,py,pz); pcimg_valign3d(cx,py,pz); pcimg_valign3d(nx,py,pz);        pcimg_valign3d(px,cy,pz); pcimg_valign3d(cx,cy,pz); pcimg_valign3d(nx,cy,pz);        pcimg_valign3d(px,ny,pz); pcimg_valign3d(cx,ny,pz); pcimg_valign3d(nx,ny,pz);        pcimg_valign3d(px,py,cz); pcimg_valign3d(cx,py,cz); pcimg_valign3d(nx,py,cz);        pcimg_valign3d(px,cy,cz);                           pcimg_valign3d(nx,cy,cz);        pcimg_valign3d(px,ny,cz); pcimg_valign3d(cx,ny,cz); pcimg_valign3d(nx,ny,cz);        pcimg_valign3d(px,py,nz); pcimg_valign3d(cx,py,nz); pcimg_valign3d(nx,py,nz);        pcimg_valign3d(px,cy,nz); pcimg_valign3d(cx,cy,nz); pcimg_valign3d(nx,cy,nz);        pcimg_valign3d(px,ny,nz); pcimg_valign3d(cx,ny,nz); pcimg_valign3d(nx,ny,nz);        const float          u0 = (float)(0.5f*dl*W.linear_pix3d(X,Y,Z,0)),          v0 = (float)(0.5f*dl*W.linear_pix3d(X,Y,Z,1)),          w0 = (float)(0.5f*dl*W.linear_pix3d(X,Y,Z,2)),          u1 = (float)(0.5f*dl*W.linear_pix3d(X+u0,Y+v0,Z+w0,0)),          v1 = (float)(0.5f*dl*W.linear_pix3d(X+u0,Y+v0,Z+w0,1)),          w1 = (float)(0.5f*dl*W.linear_pix3d(X+u0,Y+v0,Z+w0,2)),          u2 = (float)(0.5f*dl*W.linear_pix3d(X+u1,Y+v1,Z+w1,0)),          v2 = (float)(0.5f*dl*W.linear_pix3d(X+u1,Y+v1,Z+w1,1)),          w2 = (float)(0.5f*dl*W.linear_pix3d(X+u1,Y+v1,Z+w1,2)),          u3 = (float)(0.5f*dl*W.linear_pix3d(X+u2,Y+v2,Z+w2,0)),          v3 = (float)(0.5f*dl*W.linear_pix3d(X+u2,Y+v2,Z+w2,1)),          w3 = (float)(0.5f*dl*W.linear_pix3d(X+u2,Y+v2,Z+w2,2));        float          u = u0/6 + u1/3 + u2/3 + u3/6,          v = v0/6 + v1/3 + v2/3 + v3/6,          w = w0/6 + w1/3 + w2/3 + w3/6;        if ((pu*u+pv*v+pw*w)<0) { u=-u; v=-v; w=-w; }        X-=(pu=u); Y-=(pv=v); Z-=(pw=w);        tracking.insert(CImg<restype>::vector(X,Y,Z),0);      }    } break;    }  } break;  // 2D integral lines  //-------------------  case 2: {    switch (interpolation) {    case 0: { // Nearest neighbor      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        tracking.insert(CImg<restype>::vector(X,Y));        const int          cx = (int)(X+0.5f),          cy = (int)(Y+0.5f);        float          u = (float)(dl*W(cx,cy,0,0)),          v = (float)(dl*W(cx,cy,0,1));        if ((pu*u + pv*v)<0) { u=-u; v=-v; }        X+=(pu=u); Y+=(pv=v);      }      pu = cu;      pv = cv;      X  = x;      Y  = y;      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        const int          cx = (int)(X+0.5f),          cy = (int)(Y+0.5f);        float          u = (float)(dl*W(cx,cy,0,0)),          v = (float)(dl*W(cx,cy,0,1));        if ((pu*u + pv*v)<0) { u=-u; v=-v; }        X-=(pu=u); Y-=(pv=v);        tracking.insert(CImg<restype>::vector(X,Y),0);      }    } break;    case 1: { // Linear      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        tracking.insert(CImg<restype>::vector(X,Y));        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1;        const float          curru = (float)W(cx,cy,0,0),          currv = (float)W(cx,cy,0,1);        pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);        pcimg_valign2d(px,cy);                        pcimg_valign2d(nx,cy);        pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);        float          u = (float)(dl*W.linear_pix2d(X,Y,0,0)),          v = (float)(dl*W.linear_pix2d(X,Y,0,1));        if ((pu*u + pv*v)<0) { u=-u; v=-v; }        X+=(pu=u); Y+=(pv=v);      }      pu = cu;      pv = cv;      X  = x;      Y  = y;      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1;        const float          curru = (float)W(cx,cy,0,0),          currv = (float)W(cx,cy,0,1);        pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);        pcimg_valign2d(px,cy);                        pcimg_valign2d(nx,cy);        pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);        float          u = (float)(dl*W.linear_pix2d(X,Y,0,0)),          v = (float)(dl*W.linear_pix2d(X,Y,0,1));        if ((pu*u+pv*v)<0) { u=-u; v=-v; }        X-=(pu=u); Y-=(pv=v);        tracking.insert(CImg<restype>::vector(X,Y),0);      }    } break;    case 2: {  // 2nd order Runge Kutta      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        tracking.insert(CImg<restype>::vector(X,Y));        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1;        const float          curru = (float)W(cx,cy,0,0),          currv = (float)W(cx,cy,0,1);        pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);        pcimg_valign2d(px,cy);                        pcimg_valign2d(nx,cy);        pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);        const float          u0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,0)),          v0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,1));        float          u = (float)(dl*W.linear_pix2d(X+u0,Y+v0,0,0)),          v = (float)(dl*W.linear_pix2d(X+u0,Y+v0,0,1));        if ((pu*u+pv*v)<0) { u=-u; v=-v; }        X+=(pu=u); Y+=(pv=v);      }      pu = cu;      pv = cv;      X  = x;      Y  = y;      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1;        const float          curru = (float)W(cx,cy,0,0),          currv = (float)W(cx,cy,0,1);        pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);        pcimg_valign2d(px,cy);                        pcimg_valign2d(nx,cy);        pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);        const float          u0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,0)),          v0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,1));        float          u = (float)(dl*W.linear_pix2d(X+u0,Y+v0,0,0)),          v = (float)(dl*W.linear_pix2d(X+u0,Y+v0,0,1));        if ((pu*u+pv*v)<0) { u=-u; v=-v; }        X-=(pu=u); Y-=(pv=v);        tracking.insert(CImg<restype>::vector(X,Y),0);      }    } break;    case 3: {  // 4nd order Runge Kutta      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        tracking.insert(CImg<restype>::vector(X,Y));        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1;        const float          curru = (float)W(cx,cy,0,0),          currv = (float)W(cx,cy,0,1);        pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);        pcimg_valign2d(px,cy);                        pcimg_valign2d(nx,cy);        pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);        const float          u0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,0)),          v0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,1)),          u1 = (float)(0.5f*dl*W.linear_pix2d(X+u0,Y+v0,0,0)),          v1 = (float)(0.5f*dl*W.linear_pix2d(X+u0,Y+v0,0,1)),          u2 = (float)(0.5f*dl*W.linear_pix2d(X+u1,Y+v1,0,0)),          v2 = (float)(0.5f*dl*W.linear_pix2d(X+u1,Y+v1,0,1)),          u3 = (float)(0.5f*dl*W.linear_pix2d(X+u2,Y+v2,0,0)),          v3 = (float)(0.5f*dl*W.linear_pix2d(X+u2,Y+v2,0,1));        float          u = u0/6 + u1/3 + u2/3 + u3/6,          v = v0/6 + v1/3 + v2/3 + v3/6;        if ((pu*u+pv*v)<0) { u=-u; v=-v; }        X+=(pu=u); Y+=(pv=v);      }      pu = cu;      pv = cv;      X  = x;      Y  = y;      for (float l=0; l<L2 && X>=0 && X<=dx1 && Y>=0 && Y<=dy1; l+=dl) {        const int          cx = (int)X, px = (cx-1<0)?0:cx-1, nx = (cx+1>dx1)?dx1:cx+1,          cy = (int)Y, py = (cy-1<0)?0:cy-1, ny = (cy+1>dy1)?dy1:cy+1;        const float          curru = (float)W(cx,cy,0,0),          currv = (float)W(cx,cy,0,1);        pcimg_valign2d(px,py); pcimg_valign2d(cx,py); pcimg_valign2d(nx,py);        pcimg_valign2d(px,cy);                        pcimg_valign2d(nx,cy);        pcimg_valign2d(px,ny); pcimg_valign2d(cx,ny); pcimg_valign2d(nx,ny);        const float          u0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,0)),          v0 = (float)(0.5f*dl*W.linear_pix2d(X,Y,0,1)),          u1 = (float)(0.5f*dl*W.linear_pix2d(X+u0,Y+v0,0,0)),          v1 = (float)(0.5f*dl*W.linear_pix2d(X+u0,Y+v0,0,1)),          u2 = (float)(0.5f*dl*W.linear_pix2d(X+u1,Y+v1,0,0)),          v2 = (float)(0.5f*dl*W.linear_pix2d(X+u1,Y+v1,0,1)),          u3 = (float)(0.5f*dl*W.linear_pix2d(X+u2,Y+v2,0,0)),          v3 = (float)(0.5f*dl*W.linear_pix2d(X+u2,Y+v2,0,1));        float          u = u0/6 + u1/3 + u2/3 + u3/6,          v = v0/6 + v1/3 + v2/3 + v3/6;        if ((pu*u+pv*v)<0) { u=-u; v=-v; }        X-=(pu=u); Y-=(pv=v);        tracking.insert(CImg<restype>::vector(X,Y),0);      }    } break;    }  } break;  default:    throw CImgInstanceException("CImg<%s>::get_integral_line() : Instance image must have dimv()=2 or 3 (current is %u).",                                pixel_type(),dim);    break;  }  return tracking;}#endif