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

void wotreplay::read_png(std::istream &is, boost::multi_array<uint8_t, 3> &image) {
    png_structp png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr);
    png_infop info_ptr = png_create_info_struct(png_ptr);
    png_infop end_info = png_create_info_struct(png_ptr);
    if (setjmp(png_jmpbuf(png_ptr))) {
        png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
        return;
    }
    
    png_set_read_fn(png_ptr, &is, &user_read_data);

    png_read_info(png_ptr, info_ptr);
    int width = png_get_image_width(png_ptr, info_ptr);
    int height = png_get_image_height(png_ptr, info_ptr);
    png_byte color_type = png_get_color_type(png_ptr, info_ptr);
    if (color_type == PNG_COLOR_TYPE_RGB) {
        // add alpha channel
        png_set_add_alpha(png_ptr, 0xff, PNG_FILLER_AFTER);
    } else if (color_type == PNG_COLOR_TYPE_RGB_ALPHA) {
        // do nothing extra
    } else {
        throw std::runtime_error("Unsupported image type");
    }

    image.resize(boost::extents[height][width][4]);
    std::vector<png_bytep> row_pointers;
    get_row_pointers(image, row_pointers);
    png_read_image(png_ptr, &row_pointers[0]);
    png_read_end(png_ptr, info_ptr);
}

//-----------------------------------------------------------------------------
void DofMap::tabulate_coordinates(boost::multi_array<double, 2>& coordinates,
                                  const ufc::cell& ufc_cell) const
{
  // FIXME: This is a hack because UFC wants a double pointer for coordinates
  dolfin_assert(_ufc_dofmap);

  // Check dimensions
  if (coordinates.shape()[0] != cell_dimension(ufc_cell.index) ||
      coordinates.shape()[1] != _ufc_dofmap->geometric_dimension())
  {
    boost::multi_array<double, 2>::extent_gen extents;
    const std::size_t cell_dim = cell_dimension(ufc_cell.index);
    coordinates.resize(extents[cell_dim][_ufc_dofmap->geometric_dimension()]);
  }

  // Set vertex coordinates
  const std::size_t num_points = coordinates.size();
  std::vector<double*> coords(num_points);
  for (std::size_t i = 0; i < num_points; ++i)
    coords[i] = &(coordinates[i][0]);

  // Tabulate coordinates
  _ufc_dofmap->tabulate_coordinates(coords.data(),
                                    &ufc_cell.vertex_coordinates[0]);
}

    void load(input_archive& ar, boost::multi_array<T, N, Allocator>& marray, unsigned)
    {
        boost::array<std::size_t, N> shape;
        ar & shape;

        marray.resize(shape);
        ar & make_array(marray.data(), marray.num_elements());
    }

//-----------------------------------------------------------------------------
void HexahedronCell::create_entities(boost::multi_array<unsigned int, 2>& e,
                                        std::size_t dim, const unsigned int* v) const
{
  // We need to know how to create edges and faces
  switch (dim)
  {
  case 1:
    // Resize data structure
    e.resize(boost::extents[12][2]);

    // Create the 12 edges
    e[0][0]  = v[0]; e[0][1]  = v[1];
    e[1][0]  = v[2]; e[1][1]  = v[3];
    e[2][0]  = v[4]; e[2][1]  = v[5];
    e[3][0]  = v[6]; e[3][1]  = v[7];
    e[4][0]  = v[0]; e[4][1]  = v[2];
    e[5][0]  = v[1]; e[5][1]  = v[3];
    e[6][0]  = v[4]; e[6][1]  = v[6];
    e[7][0]  = v[5]; e[7][1]  = v[7];
    e[8][0]  = v[0]; e[8][1]  = v[4];
    e[9][0]  = v[1]; e[9][1]  = v[5];
    e[10][0] = v[2]; e[10][1] = v[6];
    e[11][0] = v[3]; e[11][1] = v[7];
    break;
  case 2:
    // Resize data structure
    e.resize(boost::extents[6][4]);

    // Create the 6 faces
    e[0][0] = v[0]; e[0][1] = v[2]; e[0][2] = v[4]; e[0][3] = v[6];
    e[1][0] = v[1]; e[1][1] = v[3]; e[1][2] = v[5]; e[1][3] = v[7];
    e[2][0] = v[0]; e[2][1] = v[1]; e[2][2] = v[4]; e[2][3] = v[5];
    e[3][0] = v[2]; e[3][1] = v[3]; e[3][2] = v[6]; e[3][3] = v[7];
    e[4][0] = v[0]; e[4][1] = v[1]; e[4][2] = v[2]; e[4][3] = v[3];
    e[5][0] = v[4]; e[5][1] = v[5]; e[5][2] = v[6]; e[5][3] = v[7];
    break;
  default:
    dolfin_error("HexahedronCell.cpp",
                 "create entities of tetrahedron cell",
                 "Don't know how to create entities of topological dimension %d", dim);
  }
}

//-----------------------------------------------------------------------------
void IntervalCell::create_entities(boost::multi_array<unsigned int, 2>& e,
                                   std::size_t dim, const unsigned int* v) const
{
  // For completeness, IntervalCell has two 'edges'
  dolfin_assert(dim == 0);

  // Resize data structure
  e.resize(boost::extents[2][1]);
  // Create the three edges
  e[0][0] = v[0]; e[1][0] = v[1];
}

void compute_Tnl_legendre(int n_matsubara, int n_legendre, boost::multi_array<std::complex<double>,2> &Tnl) {
  double sign_tmp = 1.0;
  Tnl.resize(boost::extents[n_matsubara][n_legendre]);
  for (int im = 0; im < n_matsubara; ++im) {
    std::complex<double> ztmp(0.0, 1.0);
    for (int il = 0; il < n_legendre; ++il) {
      Tnl[im][il] = sign_tmp * ztmp * std::sqrt(2 * il + 1.0) * boost::math::sph_bessel(il, 0.5 * (2 * im + 1) * M_PI);
      ztmp *= std::complex<double>(0.0, 1.0);
    }
    sign_tmp *= -1;
  }
}

    void precompute(double distance, double scalex, double scaley)
    {
        double const radius2 = distance * distance;

        width  = 2* ceil(distance / scalex) + 1;
        height = 2* ceil(distance / scaley) + 1;
        in_distance.resize(boost::extents[height][width]);
        std::fill(in_distance.data(), in_distance.data() + in_distance.num_elements(), false);
        parents.resize(boost::extents[height][width]);
        std::fill(parents.data(), parents.data() + parents.num_elements(), std::make_pair(-1, -1));

        centerx = width  / 2;
        centery = height / 2;
        parents[centery][centerx] = std::make_pair(-1, -1);

        for (unsigned int y = 0; y < height; ++y)
        {
            for (unsigned int x = 0; x < width; ++x)
            {
                int dx = (centerx - x);
                int dy = (centery - y);
                if (dx == 0 && dy == 0) continue;

                double d2 = dx * dx * scalex * scalex + dy * dy * scaley * scaley;
                in_distance[y][x] = (d2 < radius2);
                if (abs(dx) > abs(dy))
                {
                    int parentx = x + dx / abs(dx);
                    int parenty = y + rint(static_cast<double>(dy) / abs(dx));
                    parents[y][x] = std::make_pair(parentx, parenty);
                }
                else
                {
                    int parentx = x + rint(static_cast<double>(dx) / abs(dy));
                    int parenty = y + dy / abs(dy);
                    parents[y][x] = std::make_pair(parentx, parenty);
                }
            }
        }
    }

void load( Archive & ar,
           boost::multi_array<double,2> & t,
           const unsigned int file_version )
{
    typedef boost::multi_array<double,2> multi_array_;
    typedef typename multi_array_::size_type size_;
    size_ n0;
    ar >> BOOST_SERIALIZATION_NVP( n0 );
    size_ n1;
    ar >> BOOST_SERIALIZATION_NVP( n1 );
    t.resize( boost::extents[n0][n1] );
    ar >> make_array( t.data(), t.num_elements() );
}

void
resize_from_MultiArray(boost::multi_array<T,N>& marray, U& other) {

  // U must be a model of MultiArray
  boost::function_requires<
  boost::multi_array_concepts::ConstMultiArrayConcept<U,U::dimensionality> >();
  // U better have U::dimensionality == N
  BOOST_STATIC_ASSERT(U::dimensionality == N);

  boost::array<typename boost::multi_array<T,N>::size_type, N> shape;

  std::copy(other.shape(), other.shape()+N, shape.begin());

  marray.resize(shape);
  
}

//-----------------------------------------------------------------------------
void TriangleCell::create_entities(boost::multi_array<unsigned int, 2>&  e,
                                   std::size_t dim, const unsigned int* v) const
{
  // We only need to know how to create edges
  if (dim != 1)
  {
    dolfin_error("TriangleCell.cpp",
                 "create entities of triangle cell",
                 "Don't know how to create entities of topological dimension %d", dim);
  }

  // Resize data structure
  e.resize(boost::extents[3][2]);

  // Create the three edges
  e[0][0] = v[1]; e[0][1] = v[2];
  e[1][0] = v[0]; e[1][1] = v[2];
  e[2][0] = v[0]; e[2][1] = v[1];
}

void wotreplay::resize(boost::multi_array<uint8_t, 3> &original, int width, int height, boost::multi_array<uint8_t, 3> &result) {
    const size_t *shape = original.shape();
    result.resize(boost::extents[height][width][shape[2]]);

    float ty = ((float) shape[0] / height);
    float tx = ((float) shape[1] / width);

    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            int xx = (int) (tx * x);
            int yy = (int) (ty * y);
            float dx = (tx * x) - xx;
            float dy = (ty * y) - yy;
            
            for (int c = 0; c < shape[2]; ++c) {
                result[y][x][c] = original[yy][xx][c]*(1-dy)*(1-dx) +
                                    original[std::min(yy + 1, (int) shape[0] - 1)][xx][c]*dy*(1-dx) +
                                    original[std::min(yy + 1, (int) shape[0] - 1)][std::min(xx + 1, (int) shape[1] - 1)][c]*dy*dx +
                                    original[yy][std::min(xx + 1, (int) shape[1]  - 1)][c]*(1-dy)*dx;
            }
        }
    }
}

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

	ParseCommandLine(argc, argv);
	ParseConfigFile(config_filename);

	volume.resize(boost::extents[D][H][W]);

	FILE* in = fopen(input_filename.c_str(), "r+b");
	short *buff = new short[W];
	for(int k=0; k<D; k++){
		for(int j=0; j<H; j++){
			fread(buff, W, sizeof(short), in);
			for(int i=0; i<W; i++){
				volume[k][j][i] = (float)buff[i];// -(float)log((double)buff[i]);
			}
		}
	}
	fclose(in);


	string basename = split(input_filename, ".").at(0);

	std::ofstream ofs;
	ofs.open((basename + ".inr").c_str(), std::ios::binary);

	stringstream header;
	std::string headerEND = "##}\n";
	//Header
	header << "#INRIMAGE-4#{" << endl;
	header << "XDIM="<< W << endl;
	header << "YDIM="<< H << endl;
	header << "ZDIM="<< D << endl;
	header << "VDIM=" << 1 << endl;
	header << "VX=" << VX << endl;
	header << "VY=" << VY << endl;
	header << "VZ=" << VZ << endl;
	header << "TYPE=float" <<endl;
	header << "PIXSIZE=" << pixel_size << " bits" << endl;
	header << "CPU=decm" << endl;

	int hlen = 256 - header.str().length() - headerEND.length();
	for(int i=0; i<hlen; i++){
		header << "\n";
	}
	header << headerEND;

	ofs.write(header.str().c_str(), header.str().size());

	float *buf = new float[W];

	for(int i=0; i<D; i++){
		for(int j=0; j<H; j++){
			for(int k=0; k<W; k++){
				//ofs << volume[i][j][k];
				//buf[k] = volume[i][j][k];
				ofs.write((char*)&volume[i][j][k], sizeof(float));
			}

		}
	}

	ofs.close();
}

//-----------------------------------------------------------------------------
void MeshPartitioning::distribute_vertices(const LocalMeshData& mesh_data,
                    const boost::multi_array<std::size_t, 2>& cell_vertices,
                    std::vector<std::size_t>& vertex_indices,
                    std::map<std::size_t, std::size_t>& vertex_global_to_local,
                    boost::multi_array<double, 2>& vertex_coordinates)
{
  // This function distributes all vertices (coordinates and
  // local-to-global mapping) according to the cells that are stored on
  // each process. This happens in several stages: First each process
  // figures out which vertices it needs (by looking at its cells)
  // and where those vertices are located. That information is then
  // distributed so that each process learns where it needs to send
  // its vertices.

  // Get number of processes
  const std::size_t num_processes = MPI::num_processes();

  // Get geometric dimension
  const std::size_t gdim = mesh_data.gdim;

  // Compute which vertices we need
  std::set<std::size_t> needed_vertex_indices;
  boost::multi_array<std::size_t, 2>::const_iterator vertices;
  for (vertices = cell_vertices.begin(); vertices != cell_vertices.end();
       ++vertices)
  {
    needed_vertex_indices.insert(vertices->begin(), vertices->end());
  }

  // Compute where (process number) the vertices we need are located
  std::vector<std::vector<std::size_t> > send_vertex_indices(num_processes);
  std::vector<std::vector<std::size_t> > vertex_location(num_processes);
  std::set<std::size_t>::const_iterator required_vertex;
  for (required_vertex = needed_vertex_indices.begin();
       required_vertex != needed_vertex_indices.end(); ++required_vertex)
  {
    // Get process that has required vertex
    const std::size_t location
      = MPI::index_owner(*required_vertex, mesh_data.num_global_vertices);
    send_vertex_indices[location].push_back(*required_vertex);
    vertex_location[location].push_back(*required_vertex);
  }

  // Send required vertices to other processes, and receive back vertices
  // required by other processes.
  std::vector<std::vector<std::size_t> > received_vertex_indices;
  MPI::all_to_all(send_vertex_indices, received_vertex_indices);

  // Distribute vertex coordinates
  std::vector<std::vector<double> > send_vertex_coordinates(num_processes);
  const std::pair<std::size_t, std::size_t> local_vertex_range
    = MPI::local_range(mesh_data.num_global_vertices);
  for (std::size_t p = 0; p < num_processes; ++p)
  {
    send_vertex_coordinates[p].reserve(received_vertex_indices[p].size()*gdim);
    for (std::size_t i = 0; i < received_vertex_indices[p].size(); ++i)
    {
      dolfin_assert(received_vertex_indices[p][i] >= local_vertex_range.first
                 && received_vertex_indices[p][i] < local_vertex_range.second);
      const std::size_t location
        = received_vertex_indices[p][i] - local_vertex_range.first;
      for (std::size_t j = 0; j < gdim; ++j)
        send_vertex_coordinates[p].push_back(mesh_data.vertex_coordinates[location][j]);
    }
  }
  std::vector<std::vector<double> > received_vertex_coordinates;
  MPI::all_to_all(send_vertex_coordinates, received_vertex_coordinates);

  // Set index counters to first position in receive buffers
  std::vector<std::size_t> index_counters(num_processes, 0);

  // Clear data
  vertex_indices.clear();
  vertex_global_to_local.clear();

  // Count number of local vertices
  std::size_t num_local_vertices = 0;
  for (std::size_t p = 0; p < num_processes; ++p)
    num_local_vertices += received_vertex_coordinates[p].size()/gdim;

  // Store coordinates and construct global to local mapping
  vertex_coordinates.resize(boost::extents[num_local_vertices][gdim]);
  vertex_indices.resize(num_local_vertices);
  std::size_t v = 0;
  for (std::size_t p = 0; p < num_processes; ++p)
  {
    for (std::size_t i = 0; i < received_vertex_coordinates[p].size();
         i += gdim)
    {
      for (std::size_t j = 0; j < gdim; ++j)
        vertex_coordinates[v][j] = received_vertex_coordinates[p][i + j];

      const std::size_t global_vertex_index
        = vertex_location[p][index_counters[p]++];
      vertex_global_to_local[global_vertex_index] = v;
      vertex_indices[v] = global_vertex_index;

      ++v;
    }
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
void  MeshPartitioning::distribute_cells(const LocalMeshData& mesh_data,
                            const std::vector<std::size_t>& cell_partition,
                            std::vector<std::size_t>& global_cell_indices,
                            boost::multi_array<std::size_t, 2>& cell_vertices)
{
  // This function takes the partition computed by the partitioner
  // (which tells us to which process each of the local cells stored in
  // LocalMeshData on this process belongs. We use MPI::all_to_all to
  // redistribute all cells (the global vertex indices of all cells).

  // Number of MPI processes
  const std::size_t num_processes = MPI::num_processes();

  // Get dimensions of local mesh_data
  const std::size_t num_local_cells = mesh_data.cell_vertices.size();
  dolfin_assert(mesh_data.global_cell_indices.size() == num_local_cells);
  const std::size_t num_cell_vertices = mesh_data.num_vertices_per_cell;
  if (!mesh_data.cell_vertices.empty())
  {
    if (mesh_data.cell_vertices[0].size() != num_cell_vertices)
    {
      dolfin_error("MeshPartitioning.cpp",
                   "distribute cells",
                   "Mismatch in number of cell vertices (%d != %d) on process %d",
                   mesh_data.cell_vertices[0].size(), num_cell_vertices,
                   MPI::process_number());
    }
  }

  // Build array of cell-vertex connectivity and partition vector
  // Distribute the global cell number as well
  std::vector<std::vector<std::size_t> > send_cell_vertices(num_processes);
  for (std::size_t i = 0; i < num_local_cells; i++)
  {
    const std::size_t dest = cell_partition[i];
    send_cell_vertices[dest].push_back(mesh_data.global_cell_indices[i]);
    for (std::size_t j = 0; j < num_cell_vertices; j++)
      send_cell_vertices[dest].push_back(mesh_data.cell_vertices[i][j]);
  }

  // Distribute cell-vertex connectivity
  std::vector<std::vector<std::size_t> > received_cell_vertices(num_processes);
  MPI::all_to_all(send_cell_vertices, received_cell_vertices);

  // Count number of received cells
  std::size_t num_new_local_cells = 0;
  for (std::size_t p = 0; p < received_cell_vertices.size(); ++p)
  {
    num_new_local_cells
      += received_cell_vertices[p].size()/(num_cell_vertices + 1);
  }

  // Put mesh_data back into mesh_data.cell_vertices
  cell_vertices.resize(boost::extents[num_new_local_cells][num_cell_vertices]);
  global_cell_indices.resize(num_new_local_cells);

  // Loop over new cells
  std::size_t c = 0;
  for (std::size_t p = 0; p < num_processes; ++p)
  {
    for (std::size_t i = 0; i < received_cell_vertices[p].size();
         i += (num_cell_vertices + 1))
    {
      global_cell_indices[c] = received_cell_vertices[p][i];
      for (std::size_t j = 0; j < num_cell_vertices; ++j)
        cell_vertices[c][j] = received_cell_vertices[p][i + 1 + j];

      ++c;
    }
  }
}

void Octtree::find_cell_ranks( const boost::multi_array<Real,2>& coordinates, std::vector<Uint>& ranks )
{
  ranks.resize(coordinates.size());

  Handle< Elements > element_component;
  Uint element_idx;
  std::deque<Uint> missing_cells;

  RealVector dummy(m_dim);

  for(Uint i=0; i<coordinates.size(); ++i)
  {
    for (Uint d=0; d<m_dim; ++d)
      dummy[d] = coordinates[i][d];
    if( find_element(dummy,element_component,element_idx) )
    {
      ranks[i] = Comm::instance().rank();
    }
    else
    {
      ranks[i] = math::Consts::uint_max();
      missing_cells.push_back(i);
    }
  }

  std::vector<Real> send_coords(m_dim*missing_cells.size());
  std::vector<Real> recv_coords;

  Uint c(0);
  boost_foreach(const Uint i, missing_cells)
  {
    for(Uint d=0; d<m_dim; ++d)
      send_coords[c++]=coordinates[i][d];
  }

  for (Uint root=0; root<PE::Comm::instance().size(); ++root)
  {

    recv_coords.resize(0);
    PE::Comm::instance().broadcast(send_coords,recv_coords,root,m_dim);

    // size is only because it doesn't get resized for this rank
    std::vector<Uint> send_found(missing_cells.size(),math::Consts::uint_max());

    if (root!=Comm::instance().rank())
    {
      std::vector<RealVector> recv_coordinates(recv_coords.size()/m_dim) ;
      boost_foreach(RealVector& realvec, recv_coordinates)
          realvec.resize(m_dim);

      c=0;
      for (Uint i=0; i<recv_coordinates.size(); ++i)
      {
        for(Uint d=0; d<m_dim; ++d)
          recv_coordinates[i][d]=recv_coords[c++];
      }

      send_found.resize(recv_coordinates.size());
      for (Uint i=0; i<recv_coordinates.size(); ++i)
      {
        if( find_element(recv_coordinates[i],element_component,element_idx) )
        {
          send_found[i] = Comm::instance().rank();
        }
        else
          send_found[i] = math::Consts::uint_max();
      }
    }

    std::vector<Uint> recv_found(missing_cells.size()*Comm::instance().size());
    PE::Comm::instance().gather(send_found,recv_found,root);

    if( root==Comm::instance().rank())
    {
      const Uint stride = missing_cells.size();
      for (Uint i=0; i<missing_cells.size(); ++i)
      {
        for(Uint p=0; p<Comm::instance().size(); ++p)
        {
          ranks[missing_cells[i]] = std::min(recv_found[i+p*stride] , ranks[missing_cells[i]]);
        }
      }
    }
  }
}