C++ PCA::eigenvector方法代码示例

int main(int argc, char**argv) {
  GetPot cl(argc, argv);
  
  if(cl.search(2,"-h","--help")) {USAGE(); exit(0);}
  string pcafile=cl.follow("pca.pca",2,"-p","--pca");
  string imagefile=cl.follow("image",2,"-i","--img");

  VectorFeature vec;
  PCA pca; pca.load(pcafile);
  
  uint width, height, depth;
  if(cl.search(2,"-c","--color")) {
    depth=3;
  } else {
    depth=1;
  }
    
  if(!cl.search(2,"-w","--width") && !cl.search(2,"-h","--height")) {
    width=uint(sqrt(double(pca.dim())));
    height=width;
    if(int(height*width)!=pca.dim()) {
      ERR << "pca not for squared images, specify width or height" << endl
          << "height=" << height << "* width=" << width << "!= size=" << pca.dim() << endl;
      exit(20);
    } 
  } else {
    if(cl.search(2,"-w","--width") && !cl.search(2,"-h","--height")) {
      width=cl.follow(10,2,"-w","--width");
      height=pca.dim()/width;
      if(int(height*width)!=pca.dim()) {
        ERR << "pca images of this width, specify valid values" << endl
            << "height=" << height << "* width=" << width << "!= size=" << pca.dim() << endl;
        exit(20);
      }
    } else if(!cl.search(2,"-w","--width") && cl.search(2,"-h","--height")) {
      height=cl.follow(10,2,"-j","--height");
      width=pca.dim()/height;
      if(int(height*width)!=pca.dim()) {
        ERR << "pca images of this height, specify valid values" << endl
            << "height=" << height << "* width=" << width << "!= size=" << pca.dim() << endl;
        exit(20);
      }
    } else {
      height=cl.follow(10,2,"-j","--height");
      width=cl.follow(10,2,"-w","--width");
      if(int(height*width)!=pca.dim()) {
        ERR << "pca images of this height and width, specify valid values" << endl
            << "height=" << height << "* width=" << width << "!= size=" << pca.dim() << endl;
        exit(20);
      }
    }
  }



  vector<double> backtransformed;
  
  if(cl.search(2,"-v","--vec") && !cl.search(2,"-b","--base") && !cl.search(2,"-M","--mean")) {
    string vecfile=cl.follow("vec.vec",2,"-v","--vec");
    DBG(10) << "Loading Vectorfile " << vecfile << endl;
    vec.load(vecfile);    
    DBG(10) << "Vector to be backtransformed" ;
    for(uint i=0;i<vec.size();++i) {
      BLINK(10) << " "<< vec[i];
    }
    BLINK(10) << endl;

    if(cl.search("-n1st"))  vec[0]=0;

    if(cl.search("--nopca")) {
      backtransformed=vec.data();
    } else {
      backtransformed=pca.backTransform(vec.data());
    }
    
    DBG(10) << "Backtransformed Vector" ;
    for(uint i=0;i<backtransformed.size();++i) {
      BLINK(10) << " " <<backtransformed[i];
    }
    BLINK(10) << endl;

  } else if(cl.search(2,"-b","--base") && !cl.search(2,"-v","--vec") && !cl.search(2,"-M","--mean")) {
    uint base=cl.follow(0,2,"-b","--base");
    backtransformed=pca.eigenvector(base);
  } else if(!cl.search(2,"-b","--base") && !cl.search(2,"-v","--vec") && cl.search(2,"-M","--mean")) {
    backtransformed=pca.mean();
  } else {
    USAGE();
    exit(20);
  }
  
  ImageFeature image(width,height,depth);
  
  if(cl.search(2,"-m","--minMean")) {
    DBG(10) << "Subtracting mean" << endl;
    for(uint i=0;i<backtransformed.size();++i) {
      backtransformed[i]-pca.mean()[i];
    }
  }
  
//.........这里部分代码省略.........

ImageFeature detect(const ImageFeature& img, const PCA& pca,uint dim, uint w, uint h) {
  uint imgdim=img.xsize()*img.ysize();
  uint padx=img.xsize()+w; uint pady=img.ysize()+h;
  uint paddim=padx*pady;
  
#ifdef HAVE_FFT_LIBRARY

  // get memory for calculations
  fftw_complex *FIMG=NULL, **FFILTER=NULL, **FMULT=NULL;
  FIMG=new fftw_complex[paddim]; 
  for(uint i=0;i<paddim;++i) {FIMG[i].re=0.0;FIMG[i].im=0.0;}
  FFILTER=new fftw_complex*[dim];
  FMULT=new fftw_complex*[dim];
  for(uint i=0;i<dim;++i) {
    FFILTER[i]=new fftw_complex[paddim];
    for(uint j=0;j<imgdim;++j) {FFILTER[i][j].re=0.0;FFILTER[i][j].im=0.0;}
    
    FMULT[i]=new fftw_complex[paddim];
    for(uint j=0;j<imgdim;++j) {FMULT[i][j].re=0.0;FMULT[i][j].im=0.0;}
  }
  
  DBG(10) << "Memory allocated" << endl;
  
  vector<double> meanTrans=pca.transform(pca.mean(),dim);
  
  // create strategy for fft
  fftwnd_plan plan = fftw2d_create_plan(padx,pady, FFTW_FORWARD, FFTW_ESTIMATE  | FFTW_IN_PLACE); 
  fftwnd_plan planb = fftw2d_create_plan(padx,pady, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE); 
  
  DBG(10) << "Strategies for FFT created" << endl;
  
  //copy image into fourier transform data structure
  for(uint x=0;x<img.xsize();++x) { for(uint y=0;y<img.ysize();++y) { FIMG[y*padx+x].re=img(x,y,0); } }
  
  // fourier transform the image
  fftwnd_one(plan,FIMG,NULL);
  
  DBG(10) << "Image Transformed" << endl;
  
  
  // fourier transform the filters
  for(uint d=0;d<dim;++d) {
    for(uint x=0;x<w;++x) {
      for(uint y=0;y<h;++y) {
        uint i=y*padx+x;
        FFILTER[d][i].re=pca.eigenvector(d)[y*w+x];
      }
    }
    fftwnd_one(plan,FFILTER[d],NULL);
    DBG(10) << "Filter " << d << " transformed." << endl;
  }
  
  // multiplication in fourier domain
  for(uint d=0;d<dim;++d) {
    for(uint i=0;i<paddim;++i) {
      FMULT[d][i].re=FIMG[i].re*FFILTER[d][i].re-FIMG[i].im*FFILTER[d][i].im;
      FMULT[d][i].im=FIMG[i].re*FFILTER[d][i].im+FIMG[i].im*FFILTER[d][i].re;
    }
    DBG(10) << "Filter " << d << " applied." ;
    //fourier back transform
    fftwnd_one(planb,FMULT[d],NULL);
    BLINK(10) << "... backtransformed.." ;
    
    
    // subtract transformed mean
    for(uint i=0;i<paddim;++i) {
      FMULT[d][i].re=FMULT[d][i].re/paddim-meanTrans[d];
    }
    BLINK(10) << ". Mean subtracted." << endl;
  }
  
  
  double energyTrans, energyImg;
  
  ImageFeature fpm(img.xsize(),img.ysize(),1);
  for(uint x=0;x<img.xsize();++x) {
    for(uint y=0;y<img.ysize();++y) {
      uint i=y*padx+x;
      energyTrans=energTrans(FMULT,i,dim);
      energyImg=energImg(img,x,y,w,h,pca.mean());
      DBG(25) << VAR(x) << " " << VAR(y) << " " << VAR(energyTrans) << " " << VAR(energyImg) << endl;
      fpm(x,y,0)=energyImg-energyTrans;
    }
  }

  DBG(10) << "fpm generation." << endl;
  return fpm;
  #endif

}