C++ Array2D::resize方法代码示例

void Garbrecht_GradientTowardsLower(const Array2D<T> &elevations, const Array2D<d8_flowdir_t> &flowdirs, garbrecht_flat_type &flats, Array2D<int32_t> &inc1){
  int loops     = 0;
  int number_incremented = -1;
  std::cerr<<"Setting up the inc1 matrix..."<<std::endl;
  inc1.resize(elevations);
  inc1.setAll(0);

  while(number_incremented!=0){
    number_incremented=0;
    for(int i=0;i<(int)flats.size();i++){
      bool increment_elevation=true;
      int x=flats[i].x;
      int y=flats[i].y;
      for(int n=1;n<=8;n++){
        if( elevations(x+dx[n],y+dy[n])<elevations(x,y) &&
          flowdirs(x+dx[n],y+dy[n])!=NO_FLOW && !flowdirs.isNoData(x+dx[n],y+dy[n])
        ){
          increment_elevation=false;
          break;
        }
        else if(inc1(x+dx[n],y+dy[n])<loops &&
            elevations(x+dx[n],y+dy[n])==elevations(x,y)
        ){
          increment_elevation=false;
          break;
        }
      }
      if(increment_elevation){
        inc1(x,y)++;
        number_incremented++;
      }
    }
    loops++;
  }
}

void TA_CTI(
  const Array2D<T> &flow_accumulation,
  const Array2D<U> &riserun_slope,
        Array2D<V> &result
){
  Timer timer;

  RDLOG_ALG_NAME<<"d8_CTI";

  if(flow_accumulation.width()!=riserun_slope.width() || flow_accumulation.height()!=riserun_slope.height())
    throw std::runtime_error("Couldn't calculate CTI! The input matricies were of unequal dimensions!");

  RDLOG_PROGRESS<<"Setting up the CTI matrix..."<<std::flush;
  result.resize(flow_accumulation);
  result.setNoData(-1);  //Log(x) can't take this value of real inputs, so we're good
  RDLOG_PROGRESS<<"succeeded.";

  RDLOG_PROGRESS<<"Calculating CTI..."<<std::flush;
  timer.start();
  #pragma omp parallel for collapse(2)
  for(int x=0;x<flow_accumulation.width();x++)
    for(int y=0;y<flow_accumulation.height();y++)
      if(flow_accumulation(x,y)==flow_accumulation.noData() || riserun_slope(x,y)==riserun_slope.noData())
        result(x,y)=result.noData();
      else
        result(x,y)=log( (flow_accumulation(x,y)/flow_accumulation.getCellArea()) / (riserun_slope(x,y)+0.001) );
  RDLOG_TIME_USE<<"succeeded in "<<timer.stop()<<"s.";
}

void d8_flow_directions(
  const Array2D<T> &elevations,
        Array2D<U> &flowdirs
){
  ProgressBar progress;

  std::cerr<<"A D8 Flow Directions"<<std::endl;
  std::cerr<<"C TODO"<<std::endl;

  std::cerr<<"p Setting up the flow directions matrix..."<<std::endl;
  flowdirs.resize(elevations);
  flowdirs.setAll(NO_FLOW);
  flowdirs.setNoData(FLOWDIR_NO_DATA);

  std::cerr<<"p Calculating D8 flow directions..."<<std::endl;
  progress.start( elevations.width()*elevations.height() );
  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++){
    progress.update( y*elevations.width() );
    for(int x=0;x<elevations.width();x++)
      if(elevations(x,y)==elevations.noData())
        flowdirs(x,y) = flowdirs.noData();
      else
        flowdirs(x,y) = d8_FlowDir(elevations,x,y);
  }
  std::cerr<<"t Succeeded in = "<<progress.stop()<<" s"<<std::endl;
}

void RobotWithGeometry::GetSelfCollisionPairs(Array2D<bool>& collision) const
{
  collision.resize(geometry.size(),geometry.size(),false);
  for(int i=0;i<collision.m;i++)
    for(int j=0;j<collision.n;j++)
      if(selfCollisions(i,j) != NULL)
	collision(i,j) = true;

}

int Application::main(int argc,char *argv[])
  {
    // Process command line
    CommandLine cmd(argc,argv,"");
    if(cmd.numArgs()!=2)
      throw String("train-3state <alignment-file> <outfile>");
    String infile=cmd.arg(0);
    String outfile=cmd.arg(1);

    // Initialization
    transCounts.resize(4,4);
    emitCounts.resize(4,nAlpha,nAlpha);
    transCounts.setAllTo(0);
    emitCounts.setAllTo(0);

    // Process all alignments (pairs of sequences in fasta file)
    FastaReader reader(infile,alphabet);
    String def1, seq1, def2, seq2;
    while(reader.nextSequence(def1,seq1) &&
	  reader.nextSequence(def2,seq2)) {
      Sequence S1(seq1,alphabet), S2(seq2,alphabet);
      updateCounts(S1,S2);
    }

    // Compute HMM parameters
    PairHMM hmm(alphabet,gapSymbol,4);
    hmm.setStateType(INSERTION_STATE,PHMM_INSERT);
    hmm.setStateType(DELETION_STATE,PHMM_DELETE);
    hmm.setStateType(MATCH_STATE,PHMM_MATCH);
    hmm.setStateType(0,PHMM_OTHER);
    for(STATE i=0 ; i<4 ; ++i) {
      double sum=0;
      for(STATE j=0 ; j<4 ; ++j)
	sum+=transCounts[i][j];
      for(STATE j=0 ; j<4 ; ++j)
	hmm.setTransP(i,j,transCounts[i][j]/sum);
    }
    for(STATE k=0 ; k<4 ; ++k) {
      double sum=0;
      for(Symbol i=0 ; i<nAlpha ; ++i)
	for(Symbol j=0 ; j<nAlpha ; ++j)
	  sum+=emitCounts(k,i,j);
      if(sum==0) sum=1;
      for(Symbol i=0 ; i<nAlpha ; ++i)
	for(Symbol j=0 ; j<nAlpha ; ++j)
	  hmm.setEmitP(k,i,j,emitCounts(k,i,j)/sum);
    }      

    // Save HMM
    hmm.save(outfile);

    return 0;
  }

  ACO_Likelihooder(int numNodes,int numAlpha,Alphabet &alpha,
		   MultSeqAlignment &A,int firstCol,int lastCol,
		   RootNode *root,AlphabetMap &alphabetMap,
		   MolecularSequenceType seqType)
    : A(A), numNodes(numNodes), numAlpha(numAlpha),alpha(alpha), 
      firstCol(firstCol), lastCol(lastCol), alphabetMap(alphabetMap),
      seqType(seqType), gapSymbols(A.getGapSymbols())
  {
    numCols=lastCol-firstCol+1;
    numNmers=pow((float)alphabetMap.getRangeSize(),(float)numCols);
    L.resize(numNodes,numNmers);
  }

DLLEXPORT void rasterize_bezigon(double * alpha_arr,
                                 const double * data, const int n, const int w, const int h)
{
    //-- filter
    typedef FilterBox<1> FilterType;
    FilterType filter;

    //-- curves
    Array<Curve<2,2> > lines;
    Array<Curve<3,2> > quads;
    Array<Curve<4,2> > cubics;
    Array<Curve<3,3> > rquads;

    for (int i = 0; i < n; i += 6) {
        Curve<4,2> cubic;
        cubic.p[0].set(data[i],			data[i+1]);
        cubic.p[1].set(data[i+2],		data[i+3]);
        cubic.p[2].set(data[i+4],		data[i+5]);
        cubic.p[3].set(data[(i+6)%n],	data[(i+7)%n]);
        cubics.push_back(cubic);
    }

    //-- image
    Array2D<vect1d> img;
    int extra = (filter.W >> 1) << 1;
    img.resize(w + extra, h + extra);
    vect1d zero;
    zero = 0;
    img = zero;

    //-- color
    ColorFunction<1,1> color;
    color.comp[0].c[0] = 1;

    //-- rasterize
    RasterizerInstance<FilterType::W, 1, 1> inst(
        img, lines.s, quads.s, cubics.s, rquads.s);
    rasterize(inst, color, &lines, &filter.filter22[0],
              &quads, &filter.filter32[0], &cubics, &filter.filter42[0],
              &rquads, &filter.filter33[0]);

    //-- output
    for (int i = 0; i < h; i++) {
        for (int j = 0; j < w; j++) {
            int idx = i*w + j;
            alpha_arr[idx] = img.data[idx].getitem(0);
        }
    }
}

static inline void TerrainProcessor(F func, const Array2D<T> &elevations, const float zscale, Array2D<float> &output){
  if(elevations.getCellLengthX()!=elevations.getCellLengthY())
    RDLOG_WARN<<"Cell X and Y dimensions are not equal!";

  output.resize(elevations);
  ProgressBar progress;

  progress.start(elevations.size());
  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++){
    progress.update(y*elevations.width());
    for(int x=0;x<elevations.width();x++)
      if(elevations.isNoData(x,y))
        output(x,y) = output.noData();
      else
        output(x,y) = func(elevations,x,y,zscale);
  }
  RDLOG_TIME_USE<<"Wall-time = "<<progress.stop();
}

void FindFlats(
  const Array2D<T>   &elevations,
  Array2D<int8_t>    &flats
){
  flats.resize(elevations);
  flats.setNoData(FLAT_NO_DATA);

  ProgressBar progress;

  progress.start( elevations.size() );

  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++)
  for(int x=0;x<elevations.width();x++){
    if(elevations.isNoData(x,y)){
      flats(x,y) = FLAT_NO_DATA;
      continue;
    }

    if(elevations.isEdgeCell(x,y)){
      flats(x,y) = NOT_A_FLAT;
      continue;
    }

    //We'll now assume that the cell is a flat unless proven otherwise
    flats(x,y) = IS_A_FLAT;

    for(int n=1;n<=8;n++){
      const int nx = x+dx[n];
      const int ny = y+dy[n];
      if(elevations(nx,ny)<elevations(x,y) || elevations.isNoData(nx,ny)){
        flats(x,y) = NOT_A_FLAT;
        break;
      }
    }

    //We handled the base case just above the for loop
  }

  RDLOG_TIME_USE<<"Succeeded in = "<<progress.stop()<<" s";
}

void Garbrecht_GradientAwayFromHigher(const Array2D<T> &elevations, const Array2D<d8_flowdir_t> &flowdirs, garbrecht_flat_type &flats, Array2D<int32_t> &inc2){
  int loops              = 0;
  int number_incremented = 0;
  std::cerr<<"Setting up the inc2 matrix..."<<std::endl;
  inc2.resize(elevations);
  inc2.setAll(0);

  while(number_incremented<(int)flats.size()){
    for(int i=0;i<(int)flats.size();i++){
      int x=flats[i].x;
      int y=flats[i].y;
      if(inc2(x,y)>0){
        inc2(x,y)++;
        continue;
      }
    }
    for(int i=0;i<(int)flats.size();i++){
      bool has_higher=false,has_lower=false;
      int x=flats[i].x;
      int y=flats[i].y;
      if(inc2(x,y)>0) continue;
      for(int n=1;n<=8;n++){
        if( !has_higher &&
          (elevations(x+dx[n],y+dy[n])>elevations(x,y) ||
          inc2(x+dx[n],y+dy[n])==2)
        )
          has_higher=true;
        else if( !has_lower &&
            elevations(x+dx[n],y+dy[n])<elevations(x,y)
        )
          has_lower=true;
      }
      if(has_higher && !has_lower){
        inc2(x,y)++;
        number_incremented++;
      }
    }
    loops++;
  }
}

// static
bool PNGIO::write( const std::string& filename,
    Array2DReadView< uint8x4 > image )
{
    Array2D< uint8x4 > tmpImage;
    const uint8_t* srcPointer;

    if( !image.packed() )
    {
        tmpImage.resize( image.size() );
        copy< uint8x4 >( image, tmpImage );
        srcPointer = reinterpret_cast< const uint8_t* >( tmpImage.pointer() );
    }
    else
    {
        srcPointer = reinterpret_cast< const uint8_t* >( image.pointer() );
    }

    unsigned int errVal = lodepng::encode(
        filename, srcPointer, image.width(), image.height(), LCT_RGBA );
    bool succeeded = ( errVal == 0 );
    return succeeded;
}

void dinf_flow_directions(const Array2D<T> &elevations, Array2D<float> &flowdirs){
  ProgressBar progress;

  std::cerr<<"\nA Dinf Flow Directions"<<std::endl;
  std::cerr<<"C Tarboton, D.G. 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research. Vol. 33. pp 309-319."<<std::endl;

  std::cerr<<"p Setting up the Dinf flow directions matrix..."<<std::endl;
  flowdirs.resize(elevations);
  flowdirs.setNoData(dinf_NO_DATA);
  flowdirs.setAll(NO_FLOW);

  std::cerr<<"p Calculating Dinf flow directions..."<<std::endl;
  progress.start( elevations.size() );
  #pragma omp parallel for
  for(int y=0;y<elevations.height();y++){
    progress.update( y*elevations.width() );
    for(int x=0;x<elevations.width();x++)
      if(elevations(x,y)==elevations.noData())
        flowdirs(x,y) = flowdirs.noData();
      else
        flowdirs(x,y) = dinf_FlowDir(elevations,x,y);
  }
  std::cerr<<"t Succeeded in = "<<progress.stop()<<" s"<<std::endl;
}

//******************** M A I N ************************
void get_font(string FontName_base, char letter, int ptsize, vect2d *times)
{
    _mkdir("font_rasters");
    std::string FontName = "fonts/" + FontName_base + ".ttf";
    char buf[256];
    string letter_name;
    letter_name += letter;
    if (letter >= 'A' && letter <= 'Z')
        letter_name += letter;

    //-----------------------------------------
    // Load contour of glyph
    //-----------------------------------------

    lines.clear();
    bez2s.clear();

    FT_Face face;
    FT_Library    library;
    FT_Error error;

    error = FT_Init_FreeType( &library );
    Check(error, "", "error initializing FT lib");

    error = FT_New_Face( library, FontName.c_str(), 0, &face );
    Check(error, "",  "error loading font");

    FT_F26Dot6 font_size = ptsize*64;
    error = FT_Set_Char_Size( face, font_size, font_size, 72, 72 );
    Check(error, "", "error setting char size");

    FT_UInt   glyph_index = FT_Get_Char_Index( face, letter );
    FT_Int32  load_flags = FT_LOAD_NO_HINTING | FT_LOAD_NO_BITMAP;
    error = FT_Load_Glyph( face,  glyph_index, load_flags );
    Check(error, "", "error loading glyph");

    FT_Glyph glyph;
    error = FT_Get_Glyph( face->glyph, &glyph );
    Check(error, "", "error getting glyph");

    FT_OutlineGlyph Outg;
    error = GetOutLine(glyph, &Outg);
    Check(error,"", "error getting outline");

    // use my own callcacks to walk over the contour
    FT_Outline_Funcs func_interface;
    func_interface.shift = 0;
    func_interface.delta = 0;
    func_interface.move_to = move_to;
    func_interface.line_to = line_to;
    func_interface.conic_to = conic_to;
    func_interface.cubic_to = cubic_to;
    FT_Outline_Decompose(&Outg->outline, &func_interface, 0);

    //-----------------------------------------
    // Rasterize using FreeType
    //-----------------------------------------

    // rasterize using FreeType so that a comparison to my code can be made
    double t_ft_start = get_time();
    FT_Render_Glyph(face->glyph, FT_RENDER_MODE_NORMAL);
    double t_ft_stop = get_time();

    // save timing
    FILE *fa = fopen("timings.txt", "a");
    fprintf(fa, "time FreeType = %f\n", t_ft_stop - t_ft_start);
    fclose(fa);

    // create grid / calculate resolution
    g.max_depth = 0;
    g.grid_res = 2;

    int maxsize = max(face->glyph->bitmap.width, face->glyph->bitmap.rows);
    while (g.grid_res < maxsize)
    {
        g.grid_res *= 2;
        g.max_depth++;
    }

    vect2i offset;
    offset.set((g.grid_res - face->glyph->bitmap.width) / 2, (g.grid_res - face->glyph->bitmap.rows) / 2);

    // copy bitmap into a buffer
    Array2D<vect3ub> ftgrid;
    ftgrid.resize(g.grid_res, g.grid_res);

    for (int k = 0; k < g.grid_res*g.grid_res; k++)
        ftgrid.data[k].set(0);

    for (int j = 0; j < face->glyph->bitmap.rows; j++)
    {
        unsigned char *row = face->glyph->bitmap.buffer + (face->glyph->bitmap.rows-j-1)*face->glyph->bitmap.pitch;

        for (int i = 0; i < face->glyph->bitmap.width; i++)
        {
            ftgrid(i + offset[0], j + offset[1]) = row[i];
        }
    }

//.........这里部分代码省略.........

void pit_mask(const Array2D<elev_t> &elevations, Array2D<int32_t> &pit_mask)
{
    grid_cellz_pq<elev_t> open;
    std::queue<grid_cellz<elev_t> > pit;
    unsigned long processed_cells = 0;
    unsigned long pitc = 0;
    ProgressBar progress;

    std::cerr << "\n###Pit Mask" << std::endl;
    std::cerr << "Setting up boolean flood array matrix..." << std::flush;
    Array2D<int8_t> closed(elevations.viewWidth(), elevations.viewHeight(), false);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "Setting up the pit mask matrix..." << std::endl;
    pit_mask.resize(elevations.viewWidth(), elevations.viewHeight());
    pit_mask.setNoData(3);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "The priority queue will require approximately "
              << (elevations.viewWidth() * 2 + elevations.viewHeight() * 2)*((long)sizeof(grid_cellz<elev_t>)) / 1024 / 1024
              << "MB of RAM."
              << std::endl;
    std::cerr << "Adding cells to the priority queue..." << std::flush;
    for (int x = 0; x < elevations.viewWidth(); x++)
    {
        open.push_cell(x, 0, elevations(x, 0));
        open.push_cell(x, elevations.viewHeight() - 1, elevations(x, elevations.viewHeight() - 1));
        closed(x, 0) = true;
        closed(x, elevations.viewHeight() - 1) = true;
    }
    for (int y = 1; y < elevations.viewHeight() - 1; y++)
    {
        open.push_cell(0, y, elevations(0, y));
        open.push_cell(elevations.viewWidth() - 1, y, elevations(elevations.viewWidth() - 1, y));
        closed(0, y) = true;
        closed(elevations.viewWidth() - 1, y) = true;
    }
    std::cerr << "succeeded." << std::endl;

    std::cerr << "%%Performing the pit mask..." << std::endl;
    progress.start(elevations.viewWidth()*elevations.viewHeight());
    while (open.size() > 0 || pit.size()>0)
    {
        grid_cellz<elev_t> c;
        if (pit.size() > 0)
        {
            c = pit.front();
            pit.pop();
        }
        else
        {
            c = open.top();
            open.pop();
        }
        processed_cells++;

        for (int n = 1; n <= 8; n++)
        {
            int nx = c.x + dx[n];
            int ny = c.y + dy[n];
            if (!elevations.in_grid(nx, ny)) continue;
            if (closed(nx, ny))
                continue;

            closed(nx, ny) = true;
            if (elevations(nx, ny) <= c.z)
            {
                if (elevations(nx, ny) < c.z)
                    pit_mask(nx, ny) = 1;
                pit.push(grid_cellz<elev_t>(nx, ny, c.z));
            }
            else
            {
                pit_mask(nx, ny) = 0;
                open.push_cell(nx, ny, elevations(nx, ny));
            }
        }

        if (elevations(c.x, c.y) == elevations.noData())
            pit_mask(c.x, c.y) = pit_mask.noData();

        progress.update(processed_cells);
    }
    std::cerr << "\t\033[96msucceeded in " << progress.stop() << "s.\033[39m" << std::endl;
    std::cerr << processed_cells << " cells processed. " << pitc << " in pits." << std::endl;
}

void priority_flood_watersheds(
    Array2D<elev_t> &elevations, Array2D<int32_t> &labels, bool alter_elevations
)
{
    grid_cellz_pq<elev_t> open;
    std::queue<grid_cellz<elev_t> > pit;
    unsigned long processed_cells = 0;
    unsigned long pitc = 0, openc = 0;
    int clabel = 1;  //TODO: Thought this was more clear than zero in the results.
    ProgressBar progress;

    std::cerr << "\n###Priority-Flood+Watershed Labels" << std::endl;
    std::cerr << "Setting up boolean flood array matrix..." << std::flush;
    Array2D<int8_t> closed(elevations.viewWidth(), elevations.viewHeight(), false);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "Setting up watershed label matrix..." << std::flush;
    labels.resize(elevations.viewWidth(), elevations.viewHeight(), -1);
    labels.setNoData(-1);
    std::cerr << "succeeded." << std::endl;

    std::cerr << "The priority queue will require approximately "
              << (elevations.viewWidth() * 2 + elevations.viewHeight() * 2)*((long)sizeof(grid_cellz<elev_t>)) / 1024 / 1024
              << "MB of RAM."
              << std::endl;
    std::cerr << "Adding cells to the priority queue..." << std::endl;
    for (int x = 0; x < elevations.viewWidth(); x++)
    {
        open.push_cell(x, 0, elevations(x, 0));
        open.push_cell(x, elevations.viewHeight() - 1, elevations(x, elevations.viewHeight() - 1));
        closed(x, 0) = true;
        closed(x, elevations.viewHeight() - 1) = true;
    }
    for (int y = 1; y < elevations.viewHeight() - 1; y++)
    {
        open.push_cell(0, y, elevations(0, y));
        open.push_cell(elevations.viewWidth() - 1, y, elevations(elevations.viewWidth() - 1, y));
        closed(0, y) = true;
        closed(elevations.viewWidth() - 1, y) = true;
    }
    std::cerr << "succeeded." << std::endl;

    std::cerr << "%%Performing Priority-Flood+Watershed Labels..." << std::endl;
    progress.start(elevations.viewWidth()*elevations.viewHeight());
    while (open.size() > 0 || pit.size()>0)
    {
        grid_cellz<elev_t> c;
        if (pit.size() > 0)
        {
            c = pit.front();
            pit.pop();
            pitc++;
        }
        else
        {
            c = open.top();
            open.pop();
            openc++;
        }
        processed_cells++;

        //Since all interior cells will be flowing into a cell which has already
        //been processed, the following line identifies only the edge cells of the
        //DEM. Each edge cell seeds its own watershed/basin. The result of this will
        //be many small watersheds/basins around the edge of the DEM.
        if (labels(c.x, c.y) == labels.noData() && elevations(c.x, c.y) != elevations.noData())  //Implies a cell without a label which borders the edge of the DEM or a region of no_data
            labels(c.x, c.y) = clabel++;

        for (int n = 1; n <= 8; n++)
        {
            int nx = c.x + dx[n];
            int ny = c.y + dy[n];
            if (!elevations.in_grid(nx, ny)) continue;
            if (closed(nx, ny))
                continue;

            //Since the neighbouring cell is not closed, its flow is directed to this
            //cell. Therefore, it is part of the same watershed/basin as this cell.
            labels(nx, ny) = labels(c.x, c.y);

            closed(nx, ny) = true;
            if (elevations(nx, ny) <= c.z)
            {
                if (alter_elevations)
                    elevations(nx, ny) = c.z;
                pit.push(grid_cellz<elev_t>(nx, ny, c.z));
            }
            else
                open.push(grid_cellz<elev_t>(nx, ny, elevations(nx, ny)));
        }
        progress.update(processed_cells);
    }

    std::cerr << "\t\033[96msucceeded in " << progress.stop() << "s.\033[39m" << std::endl;

    std::cerr << processed_cells << " cells processed. "
              << pitc << " in pits, "
              << openc << " not in pits."
              << std::endl;
}