integral_line.h

this a image processing program

        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