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

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

}

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];
    }
  }
  
//.........这里部分代码省略.........

int main(int argc , char **argv) {
  
  GetPot cl(argc,argv);
  
  if(cl.search("-h")) {USAGE(); exit(0);}

  string  line;
  if(cl.search("-E")) {
    cout << "Estimating PCA" << endl;
    
    ifstream filelist(cl.follow("filelist","-E"));
    ImageFeature img;    
    getline(filelist,line);

    img.load(line,true);
    cout << line << endl;
    
    PCA pca(img.size());
    cout << img.size() << endl;
    
    pca.putData(img.layerVector(0));
    
    while(getline(filelist,line)) {
      cout << line << endl;
      img.load(line,true);
      pca.putData(img.layerVector(0));
    }
    pca.dataEnd();
    pca.save(cl.follow("covariance.pca","-c"));
    pca.calcPCA();
    pca.save(cl.follow("transformation.pca","-t"));
    filelist.close();
    
  } else if(cl.search("-T")) {
    PCA pca;
    pca.load(cl.follow("transformation.pca","-t"));
    int dim=cl.follow(20,"-d");
    
    ifstream filelist(cl.follow("filelist","-T"));
    vector<double> tmp;
    ImageFeature img;
    while(getline(filelist,line)) {
      cout << line << endl;
      img.load(line,true);
      tmp=pca.transform(img.layerVector(0),img.size());
      tmp.resize(dim);
      VectorFeature tmpvec(tmp);
      tmpvec.save(line+".pca.vec.gz");
    }
    filelist.close();
  } else if(cl.search("-B")) {
    PCA pca;
    pca.load(cl.follow("transformation.pca","-t"));
    int x=cl.follow(16,"-x");
    int y=cl.follow(16,"-y");
    ifstream filelist(cl.follow("filelist","-B"));
    vector<double> tmp;
    VectorFeature tmpvec;
    while(getline(filelist,line)) {
      cout << line << endl;
      tmpvec.load(line);
      vector<double> tmp;
      tmp=pca.backTransform(tmpvec.data());
      ImageFeature img(tmp,x,y);
      cutoff(img);
      img.save(line+".backpca.png");
    }
    filelist.close();
  } else if(cl.search("-M")) {
    ImageFeature img; img.load(cl.follow("image.png","-M"),true);
    PCA pca;
    vector<double> backproj,vec;
    pca.load(cl.follow("transformation.pca","-t"));
    int w=cl.follow(16,"-x");
    int h=cl.follow(16,"-y");
    double scalefac=cl.follow(0.8333,"-s");
    int dim=cl.follow(20,"-d");
    
    ImageFeature scimg(img);
    ImageFeature patch;
    uint minX=100000, minY=100000; 
    uint maxX=100000, maxY=100000; 
   
    while(int(scimg.xsize())>=w and int(scimg.ysize())>=h) {
      DBG(10) << VAR(scimg.xsize()) << " x " << VAR(scimg.ysize()) << endl;
      ImageFeature faceprobmap(scimg.xsize(),scimg.ysize(),1);
      double maxDist=0.0;
      double minDist=numeric_limits<double>::max();
      
      vector<double> tmpvec;

      for(uint x=0;x<scimg.xsize();++x) {
        DBG(10) << VAR(x) << endl;
        for(uint y=0;y<scimg.ysize();++y) {
          patch=getPatch(scimg,x,y,x+w,y+h);
          vec=patch.layerVector(0);
          
          vector<double> imgMinMean=vec;
          for(uint i=0;i<imgMinMean.size();++i) {
            imgMinMean[i]-=pca.mean()[i];
//.........这里部分代码省略.........