C++ Cube::at方法代码示例

int main(int argc, char **argv) {

  if(argc!=5){
    printf("Usage: cubeNonMaxSup cube.nfo theta.nfo phi.nfo output\n");
    exit(0);
  }

  Cube< float, double>* orig  = new Cube<float, double>(argv[1]);
  Cube< float, double>* theta = new Cube<float, double>(argv[2]);
  Cube< float, double>* phi   = new Cube<float, double>(argv[3]);
  Cube< float, double>* output;
  // Cube< float, double>* outputOrient;

  string outputName       = argv[4];
  // string outputOrientName = argv[5];
  if(fileExists(outputName)){
    output = new Cube<float, double>(outputName);
  } else {
    output = orig->duplicate_clean(outputName);
  }
  // if(fileExists(outputOrientName)){
    // outputOrient = new Cube<float, double>(outputOrientName);
  // } else {
    // outputOrient = orig->duplicate_clean(outputOrientName);
  // }


  // Loop for all the voxels of the cube. For now we will ignore the borders

  printf("Performing non-maxima-supression[");
#pragma omp parallel for
  for(int z = 1; z < orig->cubeDepth -1; z++){
    int  currOrient = 0;
    bool isMax = false;
    for(int y = 1; y < orig->cubeHeight -1; y++){
      for(int x = 1; x < orig->cubeWidth -1; x++){
        isMax = true;
        currOrient = computeOrientation(theta->at(x,y,z), phi->at(x,y,z));
        // outputOrient->put(x,y,z,currOrient);
        for(int neigh = 0; neigh < 8; neigh++){
          if(orig->at(x+nbrToIdx[perpNeigh[currOrient][neigh]][0],
                      y+nbrToIdx[perpNeigh[currOrient][neigh]][1],
                      z+nbrToIdx[perpNeigh[currOrient][neigh]][2])
             > orig->at(x,y,z)){
            isMax = false;
            break;
          }
        } //Neighbors
        if(isMax){
          output->put(x,y,z, orig->at(x,y,z));
        }
      }
    }
    printf("#"); fflush(stdout);
  }
  printf("]\n");

}

int main(int argc, char **argv) {

  struct arguments arguments;
  /* Default values. */
  arguments.verbose = 0;
  arguments.outputFile = "out.cl";
  arguments.cloudName = "";
  arguments.volumeX = "";
  arguments.volumeY = "";
  arguments.volumeZ = "";

  argp_parse (&argp, argc, argv, 0, 0, &arguments);

  printf ("Cloud = %s\n OutputFile = %s\n",
          arguments.cloudName.c_str(), arguments.outputFile.c_str());

  //Code starts here
  Cloud<Point3D>* cl = new Cloud<Point3D>(arguments.cloudName);
  Cloud<Point3Do>* out = new Cloud<Point3Do>(arguments.outputFile);
  Cube< float, double>* vx = new Cube<float, double>(arguments.volumeX);
  Cube< float, double>* vy = new Cube<float, double>(arguments.volumeY);
  Cube< float, double>* vz = new Cube<float, double>(arguments.volumeZ);


  vector< int > idx(3);
  vector< float > micr(3);
  float theta, phi, r;

  for(int i = 0; i < cl->points.size(); i++){
    micr[0] = cl->points[i]->coords[0];
    micr[1] = cl->points[i]->coords[1];
    micr[2] = cl->points[i]->coords[2];
    vx->micrometersToIndexes(micr, idx);

    r = sqrt( vx->at(idx[0],idx[1],idx[2])*vx->at(idx[0],idx[1],idx[2]) +
              vy->at(idx[0],idx[1],idx[2])*vy->at(idx[0],idx[1],idx[2]) +
              vz->at(idx[0],idx[1],idx[2])*vz->at(idx[0],idx[1],idx[2]) );

    theta = atan2(vy->at(idx[0],idx[1],idx[2]), vx->at(idx[0],idx[1],idx[2]));
    phi   = acos(vz->at(idx[0],idx[1],idx[2])/r);

    Point3Do* pt = new Point3Do(micr[0], micr[1], micr[2],
                                theta, phi);

    out->points.push_back(pt);
  }
  out->saveToFile(arguments.outputFile);
}

int main(int argc, char **argv) {

  if(argc!=5)
    {
      printf("Usage: cubeAt volume_name x y z\n");
      exit(0);
    }

  string volume_str = argv[1];

  Cube<uchar,ulong>* cube_test = new Cube<uchar,ulong>(volume_str,false);
  if(cube_test->type == "uchar"){
    Cube<uchar,ulong>* cube = new Cube<uchar,ulong>(volume_str);
    printf("%s ->at(%i, %i,%i) = %u\n", argv[1], atoi(argv[2]), atoi(argv[3]), atoi(argv[4]),
           cube->at(atoi(argv[2]), atoi(argv[3]), atoi(argv[4])));
  }
  if(cube_test->type == "float"){
    Cube<float,double>* cube = new Cube<float,double>(volume_str);
    printf("%s ->at(%i, %i,%i) = %f\n", argv[1], atoi(argv[2]), atoi(argv[3]), atoi(argv[4]),
           cube->at(atoi(argv[2]), atoi(argv[3]), atoi(argv[4])));

  }
}

 arma_inline elem_type at         (const u32 row, const u32 col, const u32 slice) const { return Q.at(row, col, slice); }

inline
void
op_mean::apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type            eT;
  typedef typename get_pod_type<eT>::result  T;
  
  typedef typename ProxyCube<T1>::stored_type P_stored_type;
  
  const unwrap_cube<P_stored_type> U(P.Q);
  
  const Cube<eT>& X = U.M;
  
  const uword X_n_rows   = X.n_rows;
  const uword X_n_cols   = X.n_cols;
  const uword X_n_slices = X.n_slices;
  
  if(dim == 0)
    {
    out.set_size((X_n_rows > 0) ? 1 : 0, X_n_cols, X_n_slices);
    
    if(X_n_rows == 0)  { return; }
    
    for(uword slice=0; slice < X_n_slices; ++slice)
      {
      eT* out_mem = out.slice_memptr(slice);
      
      for(uword col=0; col < X_n_cols; ++col)
        {
        out_mem[col] = op_mean::direct_mean( X.slice_colptr(slice,col), X_n_rows );
        }
      }
    }
  else
  if(dim == 1)
    {
    out.zeros(X_n_rows, (X_n_cols > 0) ? 1 : 0, X_n_slices);
    
    if(X_n_cols == 0)  { return; }
    
    for(uword slice=0; slice < X_n_slices; ++slice)
      {
      eT* out_mem = out.slice_memptr(slice);
      
      for(uword col=0; col < X_n_cols; ++col)
        {
        const eT* col_mem = X.slice_colptr(slice,col);
        
        for(uword row=0; row < X_n_rows; ++row)
          {
          out_mem[row] += col_mem[row];
          }
        }
      
      const Mat<eT> tmp('j', X.slice_memptr(slice), X_n_rows, X_n_cols);
      
      for(uword row=0; row < X_n_rows; ++row)
        {
        out_mem[row] /= T(X_n_cols);
        
        if(arma_isfinite(out_mem[row]) == false)
          {
          out_mem[row] = op_mean::direct_mean_robust( tmp, row );
          }
        }
      }
    }
  else
  if(dim == 2)
    {
    out.zeros(X_n_rows, X_n_cols, (X_n_slices > 0) ? 1 : 0);
    
    if(X_n_slices == 0)  { return; }
    
    eT* out_mem = out.memptr();
    
    for(uword slice=0; slice < X_n_slices; ++slice)
      {
      arrayops::inplace_plus(out_mem, X.slice_memptr(slice), X.n_elem_slice );
      }
    
    out /= T(X_n_slices);
    
    podarray<eT> tmp(X_n_slices);
      
    for(uword col=0; col < X_n_cols; ++col)
    for(uword row=0; row < X_n_rows; ++row)
      {
      if(arma_isfinite(out.at(row,col,0)) == false)
        {
        for(uword slice=0; slice < X_n_slices; ++slice)
          {
          tmp[slice] = X.at(row,col,slice);
          }
        
        out.at(row,col,0) = op_mean::direct_mean_robust(tmp.memptr(), X_n_slices);
        }
      }
//.........这里部分代码省略.........

int main(int argc, char **argv) {

  if(argc!=7) {
    printf("Usage: kMSTanalyzeDir directory kInit kEnd kStep groundTruth.nfo pp.asc\n");
    exit(0);
  }

  //string directory = "/home/ggonzale/mount/bbpsg1scratch/n2/3d/tree/";
  string directory = argv[1];
  Cube<uchar, ulong>* mask = new Cube<uchar, ulong>(argv[5]);
  Cube<uchar, ulong>* rendered = new Cube<uchar, ulong>(argv[5], false);
  rendered->load_volume_data("", false);
  //Cube<uchar, ulong>* rendered = mask->duplicate_clean("rendered");
  int kInit = atoi(argv[2]);;
  int kEnd  = atoi(argv[3]);
  int kStep = atoi(argv[4]);
  Neuron* gtNeuron = new Neuron(argv[6]);
  char nameGraph[2048];

  vector< vector< double > > toSave;

  //Counts all the negative points in the ground truth
  int negativePoints = 0;
  for(int z = 0; z < mask->cubeDepth; z++)
    for(int y = 0; y < mask->cubeHeight; y++)
      for(int x = 0; x < mask->cubeWidth; x++)
        if(mask->at(x,y,z) == 255)
          negativePoints++;

  double tp = 0;
  double fp = 0;
  double fn = 0;

  if(0){
    Cube<uchar, ulong>* rendered = mask->duplicate_clean("rendered");
    sprintf(nameGraph, "%s/tree/kmsts/kmst%iw.gr", directory.c_str(), 725);
    printf("%s\n", nameGraph);
    Graph< Point3Dw, EdgeW<Point3Dw> >* gr =
      new Graph< Point3Dw, EdgeW<Point3Dw> >(nameGraph);
    evaluateGraph(directory, mask, rendered, gtNeuron, gr, tp, fp, fn);
  }

  if(0){
    for(int k = kInit; k <=kEnd; k+=kStep){
      tp = 0; fp = 0; fn = 0;
      sprintf(nameGraph, "%s/tree/kmsts/kmst%iw.gr", directory.c_str(), k);
      printf("%s\n", nameGraph);
      Graph< Point3Dw, EdgeW<Point3Dw> >* gr =
        new Graph< Point3Dw, EdgeW<Point3Dw> >(nameGraph);
      evaluateGraph(directory, mask, rendered, gtNeuron, gr, tp, fp, fn);
      vector< double > it(5);
      double tn = negativePoints - fp;
      it[0] = k; it[1] = tp; it[2] = fp; it[3] = tn; it[4] = fn;
      toSave.push_back(it);
    }
    saveMatrix(toSave, directory + "/tree/kmsts/evaluation.txt");
  }

  //Evaluatoin of the mst
  if(1){
    toSave.resize(0);
    sprintf(nameGraph, "%s/tree/mstFromCptGraphw.gr", directory.c_str());
    printf("Evaluating the graph: %s\n", nameGraph);
    tp = 0; fp = 0; fn = 0;
    Graph< Point3Dw, EdgeW<Point3Dw> >* gr =
      new Graph< Point3Dw, EdgeW<Point3Dw> >(nameGraph);
    printf("Loaded the graph: %s\n", nameGraph);
    evaluateGraph(directory, mask, rendered, gtNeuron, gr, tp, fp, fn);
    vector< double > it(5);
    double tn = negativePoints - fp;
    it[0] = 0; it[1] = tp; it[2] = fp; it[3] = tn; it[4] = fn;
    toSave.push_back(it);
    saveMatrix(toSave, directory + "/tree/mstFromCptGraphwE.txt");
  }

  if(0){
    toSave.resize(0);
    sprintf(nameGraph, "%s/tree/mst/mstFromCptGraphPrunedw.gr", directory.c_str());
    printf("Evaluating the graph: %s\n", nameGraph);
    tp = 0; fp = 0; fn = 0;
    Graph< Point3Dw, EdgeW<Point3Dw> >* gr =
      new Graph< Point3Dw, EdgeW<Point3Dw> >(nameGraph);
    printf("Loaded the graph: %s\n", nameGraph);
    evaluateGraph(directory, mask, rendered, gtNeuron, gr, tp, fp, fn);
    vector< double > it(5);
    double tn = negativePoints - fp;
    it[0] = 0; it[1] = tp; it[2] = fp; it[3] = tn; it[4] = fn;
    toSave.push_back(it);
    saveMatrix(toSave, directory + "/tree/mst/mstFromCptGraphwPrunedE.txt");
  }
}