C++ CImgList::insert方法代码示例

void insert_fiber(const CImg<T>& fiber, const CImg<te>& eigen, const CImg<tc>& palette,
                  const int xm, const int ym, const int zm,
                  const float vx, const float vy, const float vz,
                  CImgList<tp>& points, CImgList<tf>& primitives, CImgList<tc>& colors) {
  const int N0 = points.size();
  float x0 = fiber(0,0), y0 = fiber(0,1), z0 = fiber(0,2), fa0 = eigen.linear_atXYZ(x0,y0,z0,12);
  points.insert(CImg<>::vector(vx*(x0  -xm),vy*(y0 - ym),vz*(z0 - zm)));
  for (int l = 1; l<fiber.width(); ++l) {
    float x1 = fiber(l,0), y1 = fiber(l,1), z1 = fiber(l,2), fa1 = eigen.linear_atXYZ(x1,y1,z1,12);
    points.insert(CImg<tp>::vector(vx*(x1 - xm),vy*(y1 - ym),vz*(z1 - zm)));
    primitives.insert(CImg<tf>::vector(N0 + l - 1,N0 + l));
    const unsigned char
      icol = (unsigned char)(fa0*255),
      r = palette(icol,0),
      g = palette(icol,1),
      b = palette(icol,2);
    colors.insert(CImg<unsigned char>::vector(r,g,b));
    x0 = x1; y0 = y1; z0 = z1; fa0 = fa1;
  }
}

void insert_ellipsoid(const CImg<t>& tensor,const float X,const float Y,const float Z,const float tfact,
                      const float vx, const float vy, const float vz,
                      CImgList<tp>& points, CImgList<tf>& faces, CImgList<tc>& colors,
                      const unsigned int res1 = 20, const unsigned int res2 = 20) {

  // Compute eigen elements
  const float l1 = tensor[0], l2 = tensor[1], l3 = tensor[2], fa = get_FA(l1,l2,l3);

  CImg<> vec = CImg<>::matrix(tensor[3],tensor[6],tensor[9],
                              tensor[4],tensor[7],tensor[10],
                              tensor[5],tensor[8],tensor[11]);
  const int
    r = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[3]),255.0f),
    g = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[4]),255.0f),
    b = (int)cimg::min(30+1.5f*cimg::abs(255*fa*tensor[5]),255.0f);

  // Define mesh points
  const unsigned int N0 = points.size;
  for (unsigned int v=1; v<res2; v++)
    for (unsigned int u=0; u<res1; u++) {
      const float
        alpha = (float)(u*2*cimg::valuePI/res1),
        beta = (float)(-cimg::valuePI/2 + v*cimg::valuePI/res2),
        x = (float)(tfact*l1*std::cos(beta)*std::cos(alpha)),
        y = (float)(tfact*l2*std::cos(beta)*std::sin(alpha)),
        z = (float)(tfact*l3*std::sin(beta));
      points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(x,y,z)).mul(CImg<tp>::vector(vx,vy,vz)));
    }
  const unsigned int N1 = points.size;
  points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(0,0,-l3*tfact)));
  points.insert((CImg<tp>::vector(X,Y,Z)+vec*CImg<tp>::vector(0,0,l3*tfact)));
  points[points.size-2](0)*=vx; points[points.size-2](1)*=vy;  points[points.size-2](2)*=vz;
  points[points.size-1](0)*=vx; points[points.size-1](1)*=vy;  points[points.size-1](2)*=vz;

  // Define mesh triangles
  for (unsigned int vv=0; vv<res2-2; vv++)
    for (unsigned int uu=0; uu<res1; uu++) {
      const int nv = (vv+1)%(res2-1), nu = (uu+1)%res1;
      faces.insert(CImg<tf>::vector(N0+res1*vv+nu,N0+res1*nv+uu,N0+res1*vv+uu));
      faces.insert(CImg<tf>::vector(N0+res1*vv+nu,N0+res1*nv+nu,N0+res1*nv+uu));
      colors.insert(CImg<tc>::vector(r,g,b));
      colors.insert(CImg<tc>::vector(r,g,b));
    }
  for (unsigned int uu=0; uu<res1; uu++) {
    const int nu = (uu+1)%res1;
    faces.insert(CImg<tf>::vector(N0+nu,N0+uu,N1));
    faces.insert(CImg<tf>::vector(N0+res1*(res2-2)+nu, N1+1,N0+res1*(res2-2)+uu));
    colors.insert(CImg<tc>::vector(r,g,b));
    colors.insert(CImg<tc>::vector(r,g,b));
  }
}

CImg<> get_fibertrack(CImg<T>& eigen,
                      const int X0, const int Y0, const int Z0, const float lmax=100,
                      const float dl=0.1f, const float FAmin=0.7f, const float cmin=0.5f) {
#define align_eigen(i,j,k) \
  { T &u = eigen(i,j,k,3), &v = eigen(i,j,k,4), &w = eigen(i,j,k,5); \
    if (u*cu + v*cv + w*cw<0) { u=-u; v=-v; w=-w; }}

  CImgList<> resf;

  // Forward tracking
  float normU = 0, normpU = 0, l = 0, X = (float)X0, Y = (float)Y0, Z = (float)Z0;
  T
    pu = eigen(X0,Y0,Z0,3),
    pv = eigen(X0,Y0,Z0,4),
    pw = eigen(X0,Y0,Z0,5);
  normpU = (float)std::sqrt(pu*pu + pv*pv + pw*pw);
  bool stopflag = false;

  while (!stopflag) {
    if (X<0 || X>eigen.width() - 1 || Y<0 || Y>eigen.height() - 1 || Z<0 || Z>eigen.depth() - 1 ||
        eigen((int)X,(int)Y,(int)Z,12)<FAmin || l>lmax) stopflag = true;
    else {
      resf.insert(CImg<>::vector(X,Y,Z));

      const int
        cx = (int)X, px = (cx - 1<0)?0:cx - 1, nx = (cx + 1>=eigen.width())?eigen.width() - 1:cx + 1,
        cy = (int)Y, py = (cy - 1<0)?0:cy - 1, ny = (cy + 1>=eigen.height())?eigen.height() - 1:cy + 1,
        cz = (int)Z, pz = (cz - 1<0)?0:cz - 1, nz = (cz + 1>=eigen.depth())?eigen.depth() - 1:cz + 1;
      const T cu = eigen(cx,cy,cz,3), cv = eigen(cx,cy,cz,4), cw = eigen(cx,cy,cz,5);

      align_eigen(px,py,pz); align_eigen(cx,py,pz); align_eigen(nx,py,pz);
      align_eigen(px,cy,pz); align_eigen(cx,cy,pz); align_eigen(nx,cy,pz);
      align_eigen(px,ny,pz); align_eigen(cx,ny,pz); align_eigen(nx,ny,pz);
      align_eigen(px,py,cz); align_eigen(cx,py,cz); align_eigen(nx,py,cz);
      align_eigen(px,cy,cz);                        align_eigen(nx,cy,cz);
      align_eigen(px,ny,cz); align_eigen(cx,ny,cz); align_eigen(nx,ny,cz);
      align_eigen(px,py,nz); align_eigen(cx,py,nz); align_eigen(nx,py,nz);
      align_eigen(px,cy,nz); align_eigen(cx,cy,nz); align_eigen(nx,cy,nz);
      align_eigen(px,ny,nz); align_eigen(cx,ny,nz); align_eigen(nx,ny,nz);

      const T
        u0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,3),
        v0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,4),
        w0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,5),
        u1 = 0.5f*dl*eigen.linear_atXYZ(X + u0,Y + v0,Z + w0,3),
        v1 = 0.5f*dl*eigen.linear_atXYZ(X + u0,Y + v0,Z + w0,4),
        w1 = 0.5f*dl*eigen.linear_atXYZ(X + u0,Y + v0,Z + w0,5),
        u2 = 0.5f*dl*eigen.linear_atXYZ(X + u1,Y + v1,Z + w1,3),
        v2 = 0.5f*dl*eigen.linear_atXYZ(X + u1,Y + v1,Z + w1,4),
        w2 = 0.5f*dl*eigen.linear_atXYZ(X + u1,Y + v1,Z + w1,5),
        u3 = 0.5f*dl*eigen.linear_atXYZ(X + u2,Y + v2,Z + w2,3),
        v3 = 0.5f*dl*eigen.linear_atXYZ(X + u2,Y + v2,Z + w2,4),
        w3 = 0.5f*dl*eigen.linear_atXYZ(X + u2,Y + v2,Z + w2,5);
      T
        u = u0/3 + 2*u1/3 + 2*u2/3 + u3/3,
        v = v0/3 + 2*v1/3 + 2*v2/3 + v3/3,
        w = w0/3 + 2*w1/3 + 2*w2/3 + w3/3;
      if (u*pu + v*pv + w*pw<0) { u = -u; v = -v; w = -w; }
      normU = (float)std::sqrt(u*u + v*v + w*w);
      const float scal = (u*pu + v*pv + w*pw)/(normU*normpU);
      if (scal<cmin) stopflag=true;

      X+=(pu=u); Y+=(pv=v); Z+=(pw=w);
      normpU = normU;
      l+=dl;
    }
  }

  // Backward tracking
  l = dl; X = (float)X0; Y = (float)Y0; Z = (float)Z0;
  pu = eigen(X0,Y0,Z0,3);
  pv = eigen(X0,Y0,Z0,4);
  pw = eigen(X0,Y0,Z0,5);
  normpU = (float)std::sqrt(pu*pu + pv*pv + pw*pw);
  stopflag = false;

  while (!stopflag) {
    if (X<0 || X>eigen.width() - 1 || Y<0 || Y>eigen.height() - 1 || Z<0 || Z>eigen.depth() - 1 ||
        eigen((int)X,(int)Y,(int)Z,12)<FAmin || l>lmax) stopflag = true;
    else {

      const int
        cx = (int)X, px = (cx - 1<0)?0:cx - 1, nx = (cx + 1>=eigen.width())?eigen.width() - 1:cx + 1,
        cy = (int)Y, py = (cy - 1<0)?0:cy - 1, ny = (cy + 1>=eigen.height())?eigen.height() - 1:cy + 1,
        cz = (int)Z, pz = (cz - 1<0)?0:cz - 1, nz = (cz + 1>=eigen.depth())?eigen.depth() - 1:cz + 1;
      const T cu = eigen(cx,cy,cz,3), cv = eigen(cx,cy,cz,4), cw = eigen(cx,cy,cz,5);

      align_eigen(px,py,pz); align_eigen(cx,py,pz); align_eigen(nx,py,pz);
      align_eigen(px,cy,pz); align_eigen(cx,cy,pz); align_eigen(nx,cy,pz);
      align_eigen(px,ny,pz); align_eigen(cx,ny,pz); align_eigen(nx,ny,pz);
      align_eigen(px,py,cz); align_eigen(cx,py,cz); align_eigen(nx,py,cz);
      align_eigen(px,cy,cz);                        align_eigen(nx,cy,cz);
      align_eigen(px,ny,cz); align_eigen(cx,ny,cz); align_eigen(nx,ny,cz);
      align_eigen(px,py,nz); align_eigen(cx,py,nz); align_eigen(nx,py,nz);
      align_eigen(px,cy,nz); align_eigen(cx,cy,nz); align_eigen(nx,cy,nz);
      align_eigen(px,ny,nz); align_eigen(cx,ny,nz); align_eigen(nx,ny,nz);

      const T
        u0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,3),
        v0 = 0.5f*dl*eigen.linear_atXYZ(X,Y,Z,4),
//.........这里部分代码省略.........

// Main procedure
//----------------
int main(int argc,char **argv) {

  // Read command line arguments.
  cimg_usage("Render an image as a surface");
  const char *file_i    = cimg_option("-i",cimg_imagepath "logo.bmp","Input image");
  const char *file_o    = cimg_option("-o",(char*)0,"Output 3D object");
  const float sigma     = cimg_option("-smooth",1.0f,"Amount of image smoothing");
  const float ratioz    = cimg_option("-z",0.25f,"Aspect ratio along z-axis");
  const unsigned int di = cimg_option("-di",10,"Step for isophote skipping");

  // Load 2D image file.
  std::fprintf(stderr,"\n- Load file '%s'",cimg::basename(file_i)); std::fflush(stderr);
  const CImg<unsigned char>
    img  = CImg<>(file_i).blur(sigma).resize(-100,-100,1,3),
    norm = img.get_norm().normalize(0,255);

  // Compute surface with triangles.
  std::fprintf(stderr,"\n- Create image surface"); std::fflush(stderr);
  CImgList<unsigned int> primitives;
  CImgList<unsigned char> colors;
  const CImg<> points = img.get_elevation3d(primitives,colors,norm*-ratioz);

  // Compute image isophotes.
  std::fprintf(stderr,"\n- Compute image isophotes"); std::fflush(stderr);
  CImgList<unsigned int> isoprimitives;
  CImgList<unsigned char> isocolors;
  CImg<> isopoints;
  for (unsigned int i = 0; i<255; i+=di) {
    CImgList<> prims;
    const CImg<> pts = norm.get_isoline3d(prims,(float)i);
    isopoints.append_object3d(isoprimitives,pts,prims);
  }
  cimglist_for(isoprimitives,l) {
    const unsigned int i0 = isoprimitives(l,0);
    const float x0 = isopoints(i0,0), y0 = isopoints(i0,1);
    const unsigned char
      r = (unsigned char)img.linear_atXY(x0,y0,0),
      g = (unsigned char)img.linear_atXY(x0,y0,1),
      b = (unsigned char)img.linear_atXY(x0,y0,2);
    isocolors.insert(CImg<unsigned char>::vector(r,g,b));
  }
  cimg_forX(isopoints,ll) isopoints(ll,2) = -ratioz*norm.linear_atXY(isopoints(ll,0),isopoints(ll,1));

  // Save object if necessary
  if (file_o) {
    std::fprintf(stderr,"\n- Save 3d object as '%s'",cimg::basename(file_o)); std::fflush(stderr);
    points.save_off(primitives,colors,file_o);
  }

  // Enter event loop
  std::fprintf(stderr,
               "\n- Enter interactive loop.\n\n"
               "Reminder : \n"
               " + Use mouse to rotate and zoom object\n"
               " + key 'F'          : Toggle fullscreen\n"
               " + key 'Q' or 'ESC' : Quit\n"
               " + Any other key    : Change rendering type\n\n"); std::fflush(stderr);
  const char *const title = "Image viewed as a surface";
  CImgDisplay disp(800,600,title,0);
  unsigned int rtype = 2;
  CImg<float> pose = CImg<float>::identity_matrix(4);

  while (!disp.is_closed()) {
    const unsigned char white[3]={ 255, 255, 255 };
    CImg<unsigned char> visu(disp.width(),disp.height(),1,3,0);
    visu.draw_text(10,10,"%s",white,0,1,24,
                rtype==0?"Points":(rtype==1?"Lines":(rtype==2?"Faces":(rtype==3?"Flat-shaded faces":
               (rtype==4?"Gouraud-shaded faces":(rtype==5?"Phong-shaded faces":"Isophotes"))))));
    static bool first_time = true;
    if (rtype==6) visu.display_object3d(disp,isopoints,isoprimitives,isocolors,first_time,1,-1,true,
                                        500.0f,0.0f,0.0f,-5000.0f,0.0f,0.0f,true,pose.data());
    else visu.display_object3d(disp,points,primitives,colors,first_time,rtype,-1,true,
                               500.0f,0.0f,0.0f,-5000.0f,0.0f,0.0f,true,pose.data());
    first_time = false;
    switch (disp.key()) {
    case 0: break;
    case cimg::keyBACKSPACE: rtype = (7 + rtype - 1)%7; break;
    case cimg::keyQ:
    case cimg::keyESC: disp.close(); break;
    case cimg::keyF:
      if (disp.is_fullscreen()) disp.resize(800,600); else disp.resize(disp.screen_width(),disp.screen_height());
      disp.toggle_fullscreen();
      break;
    default: rtype = (rtype + 1)%7; break;
    }
  }

  return 0;
}

// Main procedure
//----------------
int main (int argc, char **argv) {

  cimg_usage("Compute the skeleton of a shape, using Hamilton-Jacobi equations");

  // Read command line arguments
  cimg_help("Input/Output options\n"
            "--------------------");
  const char* file_i = cimg_option("-i",cimg_imagepath "milla.bmp","Input (black&white) image");
  const int median = cimg_option("-median",0,"Apply median filter");
  const bool invert = cimg_option("-inv",false,"Invert image values");
  const char* file_o = cimg_option("-o",(char*)0,"Output skeleton image");
  const bool display = cimg_option("-visu",true,"Display results");

  cimg_help("Skeleton computation parameters\n"
            "-------------------------------");
  const float thresh = cimg_option("-t",-0.3f,"Threshold");
  const bool curve = cimg_option("-curve",false,"Create medial curve");

  cimg_help("Torsello correction parameters\n"
            "------------------------------");
  const bool correction = cimg_option("-corr",false,"Torsello correction");
  const float dlt1 = 2;
  const float dlt2 = cimg_option("-dlt",1.0f,"Discrete step");

  // Load the image (forcing it to be scalar with 2 values { 0,1 }).
  CImg<unsigned int> image0(file_i), image = image0.get_norm().quantize(2).normalize(0.0f,1.0f);
  if (median) image.blur_median(median);
  if (invert) (image-=1)*=-1;
  if (display) (image0.get_normalize(0,255),image.get_normalize(0,255)).display("Input image - Binary image");

  // Compute distance map.
  CImgList<float> visu;
  CImg<float> distance = image.get_distance(0);
  if (display) visu.insert(distance);

  // Compute the gradient of the distance function, and the flux (divergence) of the gradient field.
  const CImgList<float> grad = distance.get_gradient("xyz");
  CImg<float> flux = image.get_flux(grad,1,1);
  if (display) visu.insert(flux);

  // Use the Torsello correction of the flux if necessary.
  if (correction) {
    CImg<float>
      logdensity = image.get_logdensity(distance,grad,flux,dlt1),
      nflux = image.get_corrected_flux(logdensity,grad,flux,dlt2);
    if (display) visu.insert(logdensity).insert(nflux);
    flux = nflux;
  }

  if (visu) {
    cimglist_apply(visu,normalize)(0,255);
    visu.display(visu.size()==2?"Distance function - Flux":"Distance function - Flux - Log-density - Corrected flux");
  }

  // Compute the skeleton
  const CImg<unsigned int> skel = image.get_skeleton(flux,distance,curve,thresh);
  if (display) {
    (image0.resize(-100,-100,1,3)*=0.7f).get_shared_channel(1)|=skel*255.0;
    image0.draw_image(0,0,0,0,image*255.0,0.5f).display("Image + Skeleton");
  }

  // Save output image if necessary.
  if (file_o) skel.save(file_o);

  return 0;
}

void display3D(CImg<T> image,
               const char *file_o,const float ratioz,const unsigned int di,
               const char *file_pose_i,const char *file_pose_o,
               unsigned int rtype,bool color_type)
{
    std::cout<<"display image as 3D surface"<<std::flush;
    const CImg<unsigned char> norm=image.get_norm().normalize(0,255);
    // Compute surface with triangles.
    std::fprintf(stderr,"\n- Create image surface");
    std::fflush(stderr);
    CImgList<unsigned int> primitives;
    ////image colors
    CImgList<unsigned char> colors;
    const CImg<> points = image.get_elevation3d(primitives,colors,norm*-ratioz);
    ////constant colors
    CImgList<unsigned char> colors2;
    colors2=colors;
    cimglist_for(colors2,l) colors2(l).fill(255);//white

    // Compute image isophotes.
    std::fprintf(stderr,"\n- Compute image isophotes");
    std::fflush(stderr);
    CImgList<unsigned int> isoprimitives;
    ////image colors
    CImgList<unsigned char> isocolors;
    CImg<> isopoints;
    for (unsigned int i = 0; i<255; i+=di) {
        CImgList<> prims;
        const CImg<> pts = norm.get_isoline3d(prims,(float)i);
        isopoints.append_object3d(isoprimitives,pts,prims);
    }
    cimglist_for(isoprimitives,l) {
        const unsigned int i0 = isoprimitives(l,0);
        const float x0 = isopoints(i0,0), y0 = isopoints(i0,1);
        const unsigned char
        r = (unsigned char)image.linear_atXY(x0,y0,0),
        g = (unsigned char)image.linear_atXY(x0,y0,1),
        b = (unsigned char)image.linear_atXY(x0,y0,2);
        isocolors.insert(CImg<unsigned char>::vector(r,g,b));
    }
    cimg_forX(isopoints,l) isopoints(l,2) = -ratioz*norm.linear_atXY(isopoints(l,0),isopoints(l,1));
    ////constant colors
    CImgList<unsigned char> isocolors2;
    isocolors2=isocolors;
    cimglist_for(isocolors2,l) isocolors2(l).fill(255);//white

    // Save object if necessary
    if (file_o)
    {
        std::fprintf(stderr,"\n- Save 3d object as '%s'",cimg::basename(file_o));
        std::fflush(stderr);
        points.save_off(primitives,colors,file_o);
    }

    //display GUI information
    std::fprintf(stderr,
                 "\n- Enter interactive loop.\n\n"
                 "GUI reminder: \n"
                 " + Use mouse to rotate and zoom object\n"
                 " + key 'F'          : toggle Fullscreen\n"
                 " + key 'Q' or 'ESC' : Quit (i.e. exit)\n"
                 " load or save file:\n"
                 " + key 'S'          : Save displayed image (i.e. 3D view)\n"
                 " + key 'O'          : save pOse (i.e. 3D view parameters)\n"
                 " + key 'R'          : Read pose (i.e. 3D view parameters)\n"
                 " render type:\n"
                 " + key 'C'          : color render (image or constant)\n"
                 " + key 'T'          : poinTs render\n"
                 " + key 'L'          : Lines render\n"
                 " + key 'A'          : fAces render\n"
                 " + key 'H'          : flat-sHaded faces render\n"
                 " + key 'G'          : Gouraud-shaded faces render\n"
                 " + key 'P'          : Phong-shaded faces render\n"
                 " + key 'I'          : Isophotes render\n"
                 " + key 'BackSpace'  : change rendering type (i.e. decrement type)\n"
                 " + Any other key    : change rendering type (i.e. increment type)\n\n"
                );
    std::fflush(stderr);
    const char *const title = "Image viewed as a surface";
    CImgDisplay disp(800,600,title,0);
    CImg<float> pose=CImg<float>::identity_matrix(4);
    //load pose if set
    if(file_pose_i) {
        std::cerr<<"- read pose from file \""<<file_pose_i<<"\"\n"<<std::flush;
        pose.load(file_pose_i);
    }
//pose.print("pose");std::cerr<<std::flush;

    //GUI loop
    while (!disp.is_closed())
    {
        const unsigned char text_color[3]= {123,234,234};
        CImg<unsigned char> visu(disp.width(),disp.height(),1,3,0);
        visu.draw_text(10,10,"%s",text_color,0,1,24,
                       rtype==0?"Points (0,T)":(rtype==1?"Lines (1,L)":(rtype==2?"Faces (2,A)":(rtype==3?"Flat-shaded faces (3,H)":
                                                (rtype==4?"Gouraud-shaded faces (4,G)":(rtype==5?"Phong-shaded faces (5,P)":"Isophotes (6,I)"))))));
        static bool first_time=(file_pose_i)?false:true;
        if (rtype==6) visu.display_object3d(disp,isopoints,isoprimitives,(color_type)?isocolors:isocolors2,first_time,1,-1,true,
                                                500.0f,0.0f,0.0f,-5000.0f,0.0f,0.0f,true,pose.data());
        else visu.display_object3d(disp,points,primitives,(color_type)?colors:colors2,first_time,rtype,-1,true,
//.........这里部分代码省略.........