int HMatrix::ParaDividing( 
		const array2d &vec2,
		array2d &product
	) const
{
	if(nColumns != (int)vec2[0].size())
	{
//		cout << "MatrixParaProduct has different columns for vec1 and vec2!" << endl;
		return -1;
	}

	product.resize( vec2.size(), array1d(nRows) );

	for (int i=0; i<(int)vec2.size(); i++)
	{
		for (int j=0; j<nRows; j++)
		{
			product[i][j] = 0.0;
			for (int k=0; k<nColumns; k++)
			{
				product[i][j] += data[j][k] / vec2[i][k];
			}
			if(product[i][j] == 0.0) product[i][j] = EPS;
		}
	}

	return 1;
}

    reference operator=(const array2d<float> arr){
      if(arr.width()!=width() || arr.height()!=height())
        throw "Error! Grids were not of the same size!";

      #pragma omp parallel for collapse(2)
      for(int x=0;x<width();x++)
      for(int y=0;y<height();y++)
        operator()(x,y)=arr(x,y);
    }

inline double
total(const array2d& u)
{
  double s = 0;
  const int nx = u.size_x();
  const int ny = u.size_y();
  for (int y = 1; y < ny - 1; y++)
    for (int x = 1; x < nx - 1; x++)
      s += u(x, y);
  return s;
}

inline void
time_step_indexed(array2d& u, const Constants& c)
{
  // compute du/dt
  array2d du(c.nx, c.ny, u.rate(), 0, u.cache_size());
  for (int y = 1; y < c.ny - 1; y++) {
    for (int x = 1; x < c.nx - 1; x++) {
      double uxx = (u(x - 1, y) - 2 * u(x, y) + u(x + 1, y)) / (c.dx * c.dx);
      double uyy = (u(x, y - 1) - 2 * u(x, y) + u(x, y + 1)) / (c.dy * c.dy);
      du(x, y) = c.dt * c.k * (uxx + uyy);
    }
  }
  // take forward Euler step
  for (uint i = 0; i < u.size(); i++)
    u[i] += du[i];
}

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

  diagnostic("Setting up the flow directions matrix...");
  flowdirs.copyprops(elevations);
  flowdirs.init(NO_FLOW);
  flowdirs.no_data=d8_NO_DATA;
  diagnostic("succeeded.\n");

  diagnostic("%%Calculating D8 flow directions...\n");
  progress.start( elevations.width()*elevations.height() );
  #pragma omp parallel for
  for(int x=0;x<elevations.width();x++){
    progress.update( x*elevations.height() );
    for(int y=0;y<elevations.height();y++)
      if(elevations(x,y)==elevations.no_data)
        flowdirs(x,y)=flowdirs.no_data;
      else
        flowdirs(x,y)=d8_FlowDir(elevations,x,y);
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());
}

inline bool SaveGrayscaleToImageFile( const array2d<U8>& texel, const std::string& filepath )
{
	int x,y;
	int width  = texel.size_x();
	int height = texel.size_y();
	const int depth = 24;

	BitmapImage img( width, height, depth );

	for( y=0; y<height ; y++ )
	{
		for( x=0; x<width; x++ )
		{
			img.SetGrayscalePixel( x, y, texel(x,y) );
		}
	}

	return img.SaveToFile( filepath );
}

inline void
time_step_iterated(array2d& u, const Constants& c)
{
  // compute du/dt
  array2d du(c.nx, c.ny, u.rate(), 0, u.cache_size());
  for (typename array2d::iterator p = du.begin(); p != du.end(); p++) {
    int x = p.i();
    int y = p.j();
    if (1 <= x && x <= c.nx - 2 &&
        1 <= y && y <= c.ny - 2) {
      double uxx = (u(x - 1, y) - 2 * u(x, y) + u(x + 1, y)) / (c.dx * c.dx);
      double uyy = (u(x, y - 1) - 2 * u(x, y) + u(x, y + 1)) / (c.dy * c.dy);
      *p = c.dt * c.k * (uxx + uyy);
    }
  }
  // take forward Euler step
  for (typename array2d::iterator p = u.begin(), q = du.begin(); p != u.end(); p++, q++)
    *p += *q;
}

inline void GetPerlinTexture( const PerlinNoiseParams& params, array2d<float>& dest )
{
	Perlin pn( params.octaves, params.freq, params.amp, params.seed );

	float (*pPerlinFunc) (Perlin&,float,float,float,float);
	pPerlinFunc = params.tilable ? TilablePerlin : StdPerlin;

	const int w = dest.size_x();
	const int h = dest.size_y();
	float min_val = FLT_MAX, max_val = -FLT_MAX;
	for( int y=0; y<h; y++ )
	{
		for( int x=0; x<w; x++ )
		{
			float fx = (float)x / (float)w;
			float fy = (float)y / (float)h;
//			float f = pn.Get( fx, fy );
			float f = pPerlinFunc( pn, fx, fy, 1.0f, 1.0f );
			min_val = take_min( min_val, f );
			max_val = take_max( max_val, f );

			dest(x,y) = f;
		}
	}

	float val_range = max_val - min_val;

//	printf( "(min,max) = (%f,%f)\n", min_val, max_val );

	float dest_min   = params.min_value;
	float dest_range = params.max_value - params.min_value;
	for( int y=0; y<h; y++ )
	{
		for( int x=0; x<w; x++ )
		{
			float normalized_val = ( dest(x,y) - min_val ) / val_range; // normalized_val is [0,1]
			dest(x,y) = dest_min + normalized_val * dest_range;
		}
	}
}

int VMatrix::Dividing(
		const array2d &scaler,
		array2d &result
	) const
{
	result.resize( nRows, array1d(nColumns) );

	for (int i=0; i<nRows; i++)
	{
		for (int j=0; j<nColumns; j++)
		{
			result[i][j] = data[i][j] / scaler[i][j];
		}
	}

	return 1;
}

bool toDLib(const ofPixels& inPix, array2d<rgb_pixel>& outPix){
    
    int width = inPix.getWidth();
    int height = inPix.getHeight();
    outPix.set_size( height, width );
    int chans = inPix.getNumChannels();
    const unsigned char* data = inPix.getData();
    for ( unsigned n = 0; n < height;n++ )
    {
        const unsigned char* v =  &data[n * width *  chans];
        for ( unsigned m = 0; m < width;m++ )
        {
            if ( chans==1 )
            {
                unsigned char p = v[m];
                assign_pixel( outPix[n][m], p );
            }
            else{
                rgb_pixel p;
                p.red = v[m*3];
                p.green = v[m*3+1];
                p.blue = v[m*3+2];
                assign_pixel( outPix[n][m], p );
            }
        }
    }
//    if(inPix.getNumChannels() == 3){
//        int h = inPix.getHeight();
//        int w = inPix.getWidth();
//        outPix.clear();
//        outPix.set_size(h,w);
//        for (int i = 0; i < h; i++) {
//            for (int j = 0; j < w; j++) {
//                
//                outPix[i][j].red = inPix.getColor(j, i).r; //inPix[i*w + j];
//                outPix[i][j].green = inPix.getColor(j, i).g; //inPix[i*w + j + 1];
//                outPix[i][j].blue = inPix.getColor(j, i).b; //inPix[i*w + j + 2];
//            }
//        }
//        return true;
//    }else{
//        return  false;
//    }
    return true;
}

static int d8_FlowDir(const array2d<T> &elevations, const int x, const int y){
  T minimum_elevation=elevations(x,y);
  int flowdir=NO_FLOW;

  if (elevations.edge_grid(x,y)){
    if(x==0 && y==0)
      return 2;
    else if(x==0 && y==elevations.height()-1)
      return 8;
    else if(x==elevations.width()-1 && y==0)
      return 4;
    else if(x==elevations.width()-1 && y==elevations.height()-1)
      return 6;
    else if(x==0)
      return 1;
    else if(x==elevations.width()-1)
      return 5;
    else if(y==0)
      return 3;
    else if(y==elevations.height()-1)
      return 7;
  }

  /*NOTE: Since the very edges of the DEM are defined to always flow outwards,
  if they have defined elevations, it is not necessary to check if a neighbour
  is IN_GRID in the following
  NOTE: It is assumed that the no_data datum is an extremely negative
  number, such that all water which makes it to the edge of the DEM's region
  of defined elevations is sucked directly off the grid, rather than piling up
  on the edges.*/
  for(int n=1;n<=8;n++)
    if(
      elevations(x+dx[n],y+dy[n])<minimum_elevation
      || (elevations(x+dx[n],y+dy[n])==minimum_elevation
            && flowdir>0 && flowdir%2==0 && n%2==1)
    ){
      minimum_elevation=elevations(x+dx[n],y+dy[n]);
      flowdir=n;
    }

  return flowdir;
}

int output_ascii_data(
  const std::string filename,
  const array2d<T> &output_grid,
  int precision=8
){
  std::ofstream fout;
  std::string outputsep=" ";
  int output_type=OUTPUT_DEM;
  Timer write_time;
  ProgressBar progress;

  write_time.start();

  diagnostic_arg("Opening ASCII output file \"%s\"...",filename.c_str());
  fout.open(filename.c_str());
  if(!fout.is_open()){
    diagnostic("failed!\n");
    exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
  }
  diagnostic("succeeded.\n");

  //OmniGlyph output
  if(filename.substr(filename.length()-4)==".omg"){
    outputsep="|";
    output_type=OUTPUT_OMG;
    diagnostic("Writing OmniGlyph file header...");
    fout<<"Contents: Pixel array"<<std::endl;
    fout<<std::endl;
    fout<<"Width:    "<<output_grid.width()<<std::endl;
    fout<<"Height:   "<<output_grid.height()<<std::endl;
    fout<<std::endl;
    fout<<"Spectral bands:   1"<<std::endl;
    fout<<"Bits per band:   32"<<std::endl;
    fout<<"Range of values:   "<<output_grid.min()<<","<<output_grid.max()<<std::endl;
    fout<<"Actual range:   "<<output_grid.no_data<<","<<output_grid.max()<<std::endl;  //TODO: Assumes no_data is a small negative value
    fout<<"Gamma exponent:   0."<<std::endl;
    fout<<"Resolution:   100 pixels per inch"<<std::endl;
    fout<<std::endl;
    fout<<"|"<<std::endl;
  } else {
    diagnostic("Writing ArcGrid ASCII file header...");
    fout<<"ncols\t\t"<<output_grid.width()<<std::endl;
    fout<<"nrows\t\t"<<output_grid.height()<<std::endl;
    fout<<"xllcorner\t"<<std::fixed<<std::setprecision(precision)<<output_grid.xllcorner<<std::endl;
    fout<<"yllcorner\t"<<std::fixed<<std::setprecision(precision)<<output_grid.yllcorner<<std::endl;
    fout<<"cellsize\t"<<std::fixed<<std::setprecision(precision)<<output_grid.cellsize<<std::endl;
    fout<<"NODATA_value\t"<<std::fixed<<std::setprecision(precision);
    if(sizeof(T)==1)  //TODO: Crude way of detecting chars and bools
      fout<<(int)output_grid.no_data<<std::endl;
    else
      fout<<output_grid.no_data<<std::endl;
  }
  diagnostic("succeeded.\n");

  diagnostic("%%Writing ArcGrid ASCII file data...\n");
  progress.start( output_grid.width()*output_grid.height() );
  fout.precision(precision);
  fout.setf(std::ios::fixed);
  for(int y=0;y<output_grid.height();y++){
    progress.update( y*output_grid.width() );
    if(output_type==OUTPUT_OMG)
      fout<<"|";
    for(int x=0;x<output_grid.width();x++)
      if(sizeof(T)==1)  //TODO: Crude way of detecting chars and bools
        fout<<(int)output_grid(x,y)<<outputsep;
      else
        fout<<output_grid(x,y)<<outputsep;
    fout<<std::endl;
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  fout.close();

  write_time.stop();
  diagnostic_arg("Write time was: %lf\n", write_time.accumulated());

  return 0;
}

namespace amorphous
{

using boost::shared_ptr;


class CStageLightAttributeHolder
{
	enum DefaultBaseEntityHandleIndex
	{
		DIRECTIONAL_LIGHT,
		POINT_LIGHT,
		SPOTLIGHT,
		HS_DIRECTIONAL_LIGHT,
		HS_POINT_LIGHT,
		NUM_DEFAULT_BASE_ENTITY_HANDLES
	};

public:

	boost::weak_ptr<CStage> m_pStage;

	BaseEntityHandle m_aBaseEntityHandle[NUM_DEFAULT_BASE_ENTITY_HANDLES];

	/// target light entity
	EntityHandle<LightEntity> m_TargetLightEntity;
};


CStageLightAttributeHolder gs_StageLightAttribute;


CCopyEntity *CreateEntityFromDesc( CCopyEntityDesc& desc )
{
	CStage *pStage = GetStageForScriptCallback();

	if( pStage )
		return pStage->CreateEntity( desc );
	else
		return NULL;
};



/// row:    holds light groups
/// column: holds object groups (lighting groups)
array2d<char> m_LightToObject;

int init_light_groups()
{
	const int num_light_groups = 16;
	const int num_ligting_groups = 16;
	m_LightToObject.resize( num_light_groups, num_ligting_groups );

	return 0;
}


inline static shared_ptr<LightEntity> GetTargetLightEntity()
{
	return gs_StageLightAttribute.m_TargetLightEntity.Get();
}


inline static void GetInvalidDesc( LightEntityDesc& desc )
{
	// set all colors to invalid values
	desc.aColor[0] = CBE_Light::ms_InvalidColor;
	desc.aColor[1] = CBE_Light::ms_InvalidColor;
	desc.aColor[2] = CBE_Light::ms_InvalidColor;
	desc.LightGroup = CBE_Light::ms_InvalidLightGroup;

	// attenuation (for point lights)
	desc.afAttenuation[0] = desc.afAttenuation[0];
	desc.afAttenuation[1] = desc.afAttenuation[1];
	desc.afAttenuation[2] = desc.afAttenuation[2];
}



using namespace py::light;


static CCopyEntity *GetEntityByName( const char* entity_name )
{
	if( GetStageForScriptCallback() )
        return GetStageForScriptCallback()->GetEntitySet()->GetEntityByName(entity_name);
	else
        return NULL;
}


PyObject* py::light::CreateDirectionalLight( PyObject* self, PyObject* args, PyObject *keywords )
{
	LightEntityDesc desc( Light::DIRECTIONAL );
	char *base_name = "";
	char *light_name = "";
	int shadow_for_light = 1; // true(1) by default
	Vector3 dir = Vector3(0,-1,0); // default direction = vertically down
	SFloatRGBAColor color = SFloatRGBAColor(1,1,1,1);
//.........这里部分代码省略.........

void d8_upslope_cells(
  int x0, int y0, int x1, int y1,
  const array2d<T> &flowdirs,array2d<U> &upslope_cells
){
  diagnostic("Setting up the upslope_cells matrix...");
  upslope_cells.copyprops(flowdirs);
  upslope_cells.init(d8_NO_DATA);
  upslope_cells.no_data=d8_NO_DATA;
  diagnostic("succeeded.\n");
  ProgressBar progress;

  std::queue<grid_cell> expansion;

  if(x0>x1){
    std::swap(x0,x1);
    std::swap(y0,y1);
  }

  //Modified Bresenham Line-Drawing Algorithm
  int deltax=x1-x0;
  int deltay=y1-y0;
  float error=0;
  float deltaerr=(float)deltay/(float)deltax;
  if (deltaerr<0)
    deltaerr=-deltaerr;
  diagnostic_arg("Line slope is %f\n",deltaerr);
  int y=y0;
  for(int x=x0;x<=x1;x++){
    expansion.push(grid_cell(x,y));
    upslope_cells(x,y)=2;
    error+=deltaerr;
    if (error>=0.5) {
      expansion.push(grid_cell(x+1,y));
      upslope_cells(x+1,y)=2;
      y+=sgn(deltay);
      error-=1;
    }
  }

  progress.start(flowdirs.data_cells);
  long int ccount=0;
  while(expansion.size()>0){
    grid_cell c=expansion.front();
    expansion.pop();

    progress.update(ccount++);

    for(int n=1;n<=8;n++)
      if(!flowdirs.in_grid(c.x+dx[n],c.y+dy[n]))
        continue;
      else if(flowdirs(c.x+dx[n],c.y+dy[n])==NO_FLOW)
        continue;
      else if(flowdirs(c.x+dx[n],c.y+dy[n])==flowdirs.no_data)
        continue;
      else if(upslope_cells(c.x+dx[n],c.y+dy[n])==upslope_cells.no_data && n==inverse_flow[flowdirs(c.x+dx[n],c.y+dy[n])]){
        expansion.push(grid_cell(c.x+dx[n],c.y+dy[n]));
        upslope_cells(c.x+dx[n],c.y+dy[n])=1;
      }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());
  diagnostic_arg("Found %ld up-slope cells.\n",ccount);
}

int write_floating_data(
  const std::string basename,
  const array2d<T> &output_grid
){
  Timer write_time;
  ProgressBar progress;
  std::string fn_header(basename), fn_data(basename);

  //TODO: The section below should work, but is something of an abomination
  if(typeid(T)==typeid(float)){
    fn_header+=".hdr";
    fn_data+=".flt";
  } else if (typeid(T)==typeid(double)){
    fn_header+=".hdr";
    fn_data+=".dflt";
  } else {
    std::cerr<<"Cannot write floating type data into this format!"<<std::endl;
    exit(-1);
  }

  write_time.start();


  {
    diagnostic_arg("Opening floating-point header file \"%s\" for writing...",fn_header.c_str());
    std::ofstream fout;
    fout.open(fn_header.c_str());
    if(!fout.is_open()){
      diagnostic("failed!\n");
      exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
    }
    diagnostic("succeeded.\n");

    diagnostic("Writing floating-point header file...");
    fout<<"ncols\t\t"<<output_grid.width()<<std::endl;
    fout<<"nrows\t\t"<<output_grid.height()<<std::endl;
    fout<<"xllcorner\t"<<std::fixed<<std::setprecision(10)<<output_grid.xllcorner<<std::endl;
    fout<<"yllcorner\t"<<std::fixed<<std::setprecision(10)<<output_grid.yllcorner<<std::endl;
    fout<<"cellsize\t"<<std::fixed<<std::setprecision(10)<<output_grid.cellsize<<std::endl;
    fout<<"NODATA_value\t"<<std::fixed<<std::setprecision(10)<<output_grid.no_data<<std::endl;
    fout<<"BYTEORDER\tLSBFIRST"<<std::endl; //TODO
    fout.close();
    diagnostic("succeeded.\n");
  }


  diagnostic_arg("Opening floating-point data file \"%s\" for writing...",fn_data.c_str());

  {
    std::ofstream fout(fn_data.c_str(), std::ios::binary | std::ios::out);
    if(!fout.is_open()){
      diagnostic("failed!\n");
      exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
    }
    diagnostic("succeeded.\n");

    diagnostic("%%Writing floating-point data file...\n");
    progress.start( output_grid.width()*output_grid.height() );
    for(int y=0;y<output_grid.height();++y){
      progress.update( y*output_grid.width() );
      for(int x=0;x<output_grid.width();++x)
        fout.write(reinterpret_cast<const char*>(&output_grid(x,y)), std::streamsize(sizeof(T)));
    }
    fout.close();
    write_time.stop();
    diagnostic_arg(SUCCEEDED_IN,progress.stop());
  }

  diagnostic_arg("Write time was: %lf\n", write_time.accumulated());

  return 0;
}

int read_floating_data(
  const std::string basename,
  array2d<T> &grid
){
  Timer io_time;
  ProgressBar progress;
  std::string fn_header(basename), fn_data(basename);

  //TODO: The section below should work, but is something of an abomination
  if(typeid(T)==typeid(float)){
    fn_header+=".hdr";
    fn_data+=".flt";
  } else if (typeid(T)==typeid(double)){
    fn_header+=".hdr";
    fn_data+=".dflt";
  } else {
    std::cerr<<"Cannot read floating type data into this format!"<<std::endl;
    exit(-1);
  }

  int columns, rows;
  std::string byteorder;

  io_time.start();


  {
    std::ifstream fin;
    diagnostic_arg("Opening floating-point header file \"%s\" for reading...",fn_header.c_str());
    fin.open(fn_header.c_str());
    if(fin==NULL){
      diagnostic("failed!\n");
      exit(-1);
    }
    diagnostic("succeeded.\n");


    diagnostic("Reading DEM header...");
    fin>>must_be("ncols")         >>columns;
    fin>>must_be("nrows")         >>rows;
    fin>>must_be("xllcorner")     >>grid.xllcorner;
    fin>>must_be("yllcorner")     >>grid.yllcorner;
    fin>>must_be("cellsize")      >>grid.cellsize;
    fin>>must_be("NODATA_value")  >>grid.no_data;
    fin>>must_be("BYTEORDER")     >>byteorder;
    diagnostic("succeeded.\n");
    fin.close();
  }

  diagnostic_arg("The loaded DEM will require approximately %ldMB of RAM.\n",columns*rows*((long)sizeof(float))/1024/1024);

  diagnostic("Resizing grid...");  //TODO: Consider abstracting this block
  grid.resize(columns,rows);
  diagnostic("succeeded.\n");



  diagnostic_arg("Opening floating-point data file \"%s\" for reading...",fn_data.c_str());

  {
    std::ifstream fin(fn_data.c_str(), std::ios::binary | std::ios::in);
    if(!fin.is_open()){
      diagnostic("failed!\n");
      exit(-1);  //TODO: Need to make this safer! Don't just close after all that work!
    }
    diagnostic("succeeded.\n");


    diagnostic("%%Reading data...\n");
    progress.start(columns*rows);
    grid.data_cells=0;
    for(int y=0;y<rows;++y){
      progress.update(y*columns); //Todo: Check to see if ftell fails here?
      for(int x=0;x<columns;++x){
        fin.read(reinterpret_cast<char*>(&grid(x,y)), std::streamsize(sizeof(T)));
        if(grid(x,y)!=grid.no_data)
          grid.data_cells++;
      }
    }
    io_time.stop();
    diagnostic_arg(SUCCEEDED_IN,progress.stop());

  }

  diagnostic_arg(
    "Read %ld cells, of which %ld contained data (%ld%%).\n",
    grid.width()*grid.height(), grid.data_cells,
    grid.data_cells*100/grid.width()/grid.height()
  );

  diagnostic_arg("Read time was: %lf\n", io_time.accumulated());

  return 0;
}

int load_ascii_data(std::string filename, array2d<T> &elevations){
  std::ifstream fin;
  size_t file_size;
  int rows,columns;
  Timer load_time;
  ProgressBar progress;

  load_time.start();

  diagnostic_arg("Opening input ASCII-DEM file \"%s\"...",filename.c_str());
  fin.open(filename.c_str());
  if(!fin.good()){
    diagnostic("failed!\n");
    exit(-1);
  }
  diagnostic("succeeded.\n");

  diagnostic("Calculating file size...");
  fin.seekg(0, fin.end);
  file_size=fin.tellg();
  fin.seekg(0, fin.beg);
  diagnostic("succeeded.\n");

//  posix_fadvise(fileno(fin),0,0,POSIX_FADV_SEQUENTIAL);

  diagnostic("Reading DEM header...");
  fin>>must_be("ncols")         >>columns;
  fin>>must_be("nrows")         >>rows;
  fin>>must_be("xllcorner")     >>elevations.xllcorner;
  fin>>must_be("yllcorner")     >>elevations.yllcorner;
  fin>>must_be("cellsize")      >>elevations.cellsize;
  try {
    fin>>must_be("NODATA_value")  >>elevations.no_data;
  } catch (std::string e) {
    std::cerr<<e<<std::endl;
    std::cerr<<"Continuing without a NoData value!"<<std::endl;
  }
  diagnostic("succeeded.\n");

  diagnostic_arg("The loaded DEM will require approximately %ldMB of RAM.\n",columns*rows*((long)sizeof(float))/1024/1024);

  diagnostic("Resizing elevation matrix...");  //TODO: Consider abstracting this block
  elevations.resize(columns,rows);
  diagnostic("succeeded.\n");

  diagnostic("%%Reading elevation matrix...\n");
  progress.start(file_size);

  elevations.data_cells=0;
  for(int y=0;y<rows;y++){
    progress.update(fin.tellg()); //Todo: Check to see if ftell fails here?
    for(int x=0;x<columns;x++){
      fin>>elevations(x,y);
      if(elevations(x,y)!=elevations.no_data)
        elevations.data_cells++;
    }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  fin.close();

  diagnostic_arg(
    "Read %ld cells, of which %ld contained data (%ld%%).\n",
    elevations.width()*elevations.height(), elevations.data_cells,
    elevations.data_cells*100/elevations.width()/elevations.height()
  );

  load_time.stop();
  diagnostic_arg("Read time was: %lfs\n", load_time.accumulated());

  return 0;
}

/// Sets up char rects as well
void RenderUTF8TextToBuffer( const FT_Face& face,
						const std::string &text,
						int char_height,
						array2d<U8>& dest_bitmap_buffer,                 ///< [out] buffer to render the text to
						vector<TrueTypeTextureFont::CharRect>& char_rect
						)
{
	LOG_FUNCTION_SCOPE();

	float top, left, bottom, right;
	int _top, _left, _bottom, _right;
	FT_GlyphSlot slot = face->glyph; // a small shortcut
//	FT_UInt glyph_index;
	int pen_x, pen_y, n;

	const int num_chars = (int)text.size();

	int error = 0;

	vector<U32> utf8_code_points;
	string::const_iterator itr = text.begin();
	while( itr != text.end() )
	{
		U32 cp = utf8::next( itr, text.end() );
		utf8_code_points.push_back( cp );
	}
	const int num_utf8_chars = (int)utf8_code_points.size();

	char_rect.resize( num_utf8_chars );

	/// Shouldn't we change the character codes?
//	for( size_t i=0; i<text.size(); i++ )
//	{}

	int img_width = 1600;
	int img_height= 256;
	dest_bitmap_buffer.resize( img_width, img_height, 0 );

	int sx = 0;
	int sy = char_height;
	pen_x = sx;
	pen_y = sy;
	int margin = 4;
	for ( n = 0; n < num_utf8_chars; n++ )
	{
		/* load glyph image into the slot (erase previous one) */
		FT_ULong char_code = utf8_code_points[n];
//		FT_ULong char_code = text[n];
		error = FT_Load_Char( face, char_code, FT_LOAD_RENDER );

		if ( error )
			continue;

		if( slot->bitmap_left < 0 )
		{
			// some chars have negative left offset
			// - add the offset to make sure that x coords of their bounding box 
			//   do not overlap with those of other chars
			pen_x += (int)slot->bitmap_left * (-1);
		}

		// ignore errors

		_left = pen_x + slot->bitmap_left;
		_top  = pen_y - slot->bitmap_top;
		clamp( _top, 0, 4096 );

		left = (float)pen_x / (float)img_width;
		top  = (float)(pen_y - slot->bitmap_top) / (float)img_width;


		// now, draw to our target surface
		DrawBitmap( &slot->bitmap, pen_x + slot->bitmap_left, pen_y - slot->bitmap_top, dest_bitmap_buffer );

		// increment pen position
		pen_x += slot->bitmap_left + slot->bitmap.width;//slot->advance.x >> 6;

		right  = (float)pen_x / (float)img_width;
		bottom = (float)pen_y / (float)img_width;

		_right  = _left + slot->bitmap.width;
		_bottom = _top  + slot->bitmap.rows;

		pen_x += margin;

		if( img_width - char_height * 2 < pen_x )
		{
			// line feed
			pen_x = sx;
			pen_y += (int)(char_height * 1.5f); // Not sure if x1.5 is always sufficient
		}

//		DrawRect( dest_bitmap_buffer, RectLTRB( _left, _top, _right, _bottom ), 0x90 );

		char_rect[n].tex_min = TEXCOORD2( (float)_left,  (float)_top )    / (float)img_width;
		char_rect[n].tex_max = TEXCOORD2( (float)_right, (float)_bottom ) / (float)img_width;

		char_rect[n].rect.vMin = Vector2( (float)slot->bitmap_left,                      (float)char_height - slot->bitmap_top );
		char_rect[n].rect.vMax = Vector2( (float)slot->bitmap_left + slot->bitmap.width, (float)char_height - slot->bitmap_top + slot->bitmap.rows );
//.........这里部分代码省略.........

void d8_upslope_area(const array2d<T> &flowdirs, array2d<U> &area){
  char_2d dependency;
  std::queue<grid_cell> sources;
  ProgressBar progress;

  diagnostic("\n###D8 Upslope Area\n");

  diagnostic_arg(
    "The sources queue will require at most approximately %ldMB of RAM.\n",
    flowdirs.width()*flowdirs.height()*((long)sizeof(grid_cell))/1024/1024
  );

  diagnostic("Resizing dependency matrix...");
  dependency.copyprops(flowdirs);
  diagnostic("succeeded.\n");

  diagnostic("Setting up the area matrix...");
  area.copyprops(flowdirs);
  area.init(0);
  area.no_data=d8_NO_DATA;
  diagnostic("succeeded.\n");

  diagnostic("%%Calculating dependency matrix & setting no_data cells...\n");
  progress.start( flowdirs.width()*flowdirs.height() );
  #pragma omp parallel for
  for(int x=0;x<flowdirs.width();x++){
    progress.update( x*flowdirs.height() );
    for(int y=0;y<flowdirs.height();y++){
      dependency(x,y)=0;
      if(flowdirs(x,y)==flowdirs.no_data){
        area(x,y)=area.no_data;
        continue;
      }
      for(int n=1;n<=8;n++)
        if(!flowdirs.in_grid(x+dx[n],y+dy[n]))
          continue;
        else if(flowdirs(x+dx[n],y+dy[n])==NO_FLOW)
          continue;
        else if(flowdirs(x+dx[n],y+dy[n])==flowdirs.no_data)
          continue;
        else if(n==inverse_flow[(int)flowdirs(x+dx[n],y+dy[n])])
          ++dependency(x,y);
    }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  diagnostic("%%Locating source cells...\n");
  progress.start( flowdirs.width()*flowdirs.height() );
  for(int x=0;x<flowdirs.width();x++){
    progress.update( x*flowdirs.height() );
    for(int y=0;y<flowdirs.height();y++)
      if(flowdirs(x,y)==flowdirs.no_data)
        continue;
      else if(dependency(x,y)==0)
        sources.push(grid_cell(x,y));
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());

  diagnostic("%%Calculating up-slope areas...\n");
  progress.start(flowdirs.data_cells);
  long int ccount=0;
  while(sources.size()>0){
    grid_cell c=sources.front();
    sources.pop();

    ccount++;
    progress.update(ccount);

    area(c.x,c.y)+=1;

    if(flowdirs(c.x,c.y)==NO_FLOW)
      continue;

    int nx=c.x+dx[(int)flowdirs(c.x,c.y)];
    int ny=c.y+dy[(int)flowdirs(c.x,c.y)];
    if(flowdirs.in_grid(nx,ny) && area(nx,ny)!=area.no_data){
      area(nx,ny)+=area(c.x,c.y);
      if((--dependency(nx,ny))==0)
        sources.push(grid_cell(nx,ny));
    }
  }
  diagnostic_arg(SUCCEEDED_IN,progress.stop());
}

int init_light_groups()
{
	const int num_light_groups = 16;
	const int num_ligting_groups = 16;
	m_LightToObject.resize( num_light_groups, num_ligting_groups );

	return 0;
}