void RKCell::preUpdate(double epsilon[], BaseCell* neighbors[], Grid *grid, Vector3i position) {
    int model = Shared::instance()->getGridConfig()->getModel();
    MyVector3D colorGradient = MyVector3D();
    for (int a = 0; a < model; a++) {
        Vector3i neighbor = position + LBUtil::C[model][a];
        RKCell *cell = dynamic_cast<RKCell*>(grid->getGrid(neighbor.getY(), neighbor.getX(), neighbor.getZ()));
        if (cell != 0) {
            colorGradient = colorGradient + (LBUtil::C[model][a] ^ cell->getColor());
        }
    }
    for (int i = 0; i < model; i++) {
        double fi = red[i] + blue[i];
        double p = redP + blueP;
        double nextF = fi + 1.0 / epsilon[0] * (LBUtil::f_eq(u, p, model, i) - fi);
        double cos = (LBUtil::C[model][i] * colorGradient);
        if (colorGradient.norm() > 1e-10 && LBUtil::C[model][i].norm() > 1e-10) {
            cos = cos / (LBUtil::C[model][i].norm() * colorGradient.norm());
        } else {
            //cos = std::sqrt(2.0) / 2;
            //cos = 1;
            cos = 0;
        }
        GridConfig *config = Shared::instance()->getGridConfig();
        nextF += config->getRkA() * colorGradient.norm() * (2 * cos * cos - 1);
        double diff = config->getRkBeta() * redP * blueP / (p * p) * LBUtil::W[model][i] * p * cos;
        //double diff = 0;
        double newNextRed = redP / p * nextF + diff;
        double newNextBlue = blueP / p * nextF - diff;
        if (neighbors[i] != 0) {
            neighbors[i]->setNextF(i, newNextRed, 0);
            neighbors[i]->setNextF(i, newNextBlue, 1);
        }
    }
}

//----------------------------------------------------------------------------
void Texture::set(void* iPtr, const Vector3i& iS,
    GLenum iInternalFormat,
    GLenum iFormat,
    GLenum iDataType)
{
    vector<int> s(3, 0); s[0] = iS.x(); s[1] = iS.y(); s[2] = iS.z();
  set( iPtr, s, iInternalFormat, iFormat, iDataType );
}

coord_t Box::volume() const {
    if (!valid()) {
        return 0;
    }

    Vector3i diagonal = this->diagonal();

    return diagonal.x() * diagonal.y() * diagonal.z();
}

void DepositionCell::postUpdate(Grid *grid, Vector3i position) {
    for (int i = 0; i < 9; i++) {
        Vector3i neighbor = position + LBUtil::C[9][i];
        nextG[i] = g[i] + deposited * (grid->getGrid(neighbor.getY(), neighbor.getX(), neighbor.getZ())->getF(LBUtil::OPPOSITE[9][i], 1) - g[i]);
    }
    for (int i = 0; i < 9; i++) {
        g[i] = nextG[i];
    }
}

double Cube::value(const Vector3i &pos) const
{
  unsigned int index = pos.x() * m_points.y() * m_points.z()
      + pos.y() * m_points.z() + pos.z();
  if (index < m_data.size())
    return m_data[index];
  else
    return 6969.0;
}

int main(int, char**)
{
  cout.precision(3);
  Vector3i v = Vector3i::Random();
cout << "Here is the vector v:" << endl << v << endl;
cout << "v.rowwise().replicate(5) = ..." << endl;
cout << v.rowwise().replicate(5) << endl;

  return 0;
}

//------------------------------------------------------------------------------
IntAABB Octree::calculateTotalKeyBounds()
{
    IntAABB bounds;
    for (auto it = m_roots.begin(); it != m_roots.end(); ++it)
    {
        OctreeNode & node = *(it->second);
        Vector3i pos = it->first;
        bounds.addPoint(pos.x(), pos.y());
    }
    return bounds;
}

uint getOctreeDepthForBounding(const Vector3i& maxXYZ) {
    coord_t max = ::std::max(maxXYZ.x(), ::std::max(maxXYZ.y(), maxXYZ.z()));

    if (max < 0) {
        throw ::std::runtime_error("Invalid bounding (all components must not be negative)");
    }

    // how many bits are required to store numbers from 0 to max
    int requiredBits = getMostSignigicantSetBitPos(max) + 1;

    return static_cast<uint>(requiredBits);
}

void BotBotCollisionGrid::Setup(Vector3i Dimensions)
{
    assert(Dimensions.x() > 0 && "Attempting to create a world grid with zero x dimensions");
    assert(Dimensions.y() > 0 && "Attempting to create a world grid with zero y dimensions");
    
    SetupLevel(Grid1, Dimensions, 64);
    SetupLevel(Grid2, Dimensions, 128);
    SetupLevel(Grid3, Dimensions, 256);
    SetupLevel(Grid4, Dimensions, 512);

    MasterGrid.clear();    
}

void BotBotCollisionGrid::SetupLevel(vector< vector< list< Robot* > > > &Grid, Vector3i Dimensions, int GridSize)
{
    Grid.resize( (Dimensions.x() + GridSize - 1) / GridSize);//round up the number of dimensions needed
    for(unsigned int x = 0; x < Grid.size(); x++)
    {
        Grid[x].resize((Dimensions.y() + GridSize - 1) / GridSize);
        for(unsigned int y = 0; y < Grid[x].size(); y++)
        {
            Grid[x][y].clear();
        }
    }
}

void AbstractTexture::storageImplementationDefault(GLenum target, GLsizei levels, AbstractTexture::InternalFormat internalFormat, const Vector3i& size) {
    bindInternal();
    /** @todo Re-enable when extension wrangler is available for ES2 */
    #ifndef MAGNUM_TARGET_GLES2
    glTexStorage3D(target, levels, GLenum(internalFormat), size.x(), size.y(), size.z());
    #else
    //glTexStorage3DEXT(target, levels, GLenum(internalFormat), size.x(), size.y(), size.z());
    static_cast<void>(target);
    static_cast<void>(levels);
    static_cast<void>(internalFormat);
    static_cast<void>(size);
    #endif
}

bool Cube::setLimits(const Vector3 &min_, const Vector3i &dim,
                     const Vector3 &spacing_)
{
  Vector3 max_ = Vector3(min_.x() + (dim.x()-1) * spacing_[0],
                         min_.y() + (dim.y()-1) * spacing_[1],
                         min_.z() + (dim.z()-1) * spacing_[2]);
  m_min = min_;
  m_max = max_;
  m_points = dim;
  m_spacing = spacing_;
  m_data.resize(m_points.x() * m_points.y() * m_points.z());
  return true;
}

bool Cube::setLimits(const Vector3d &min_, const Vector3i &dim,
                     double spacing_)
{
  Vector3d max_ = Vector3d(min_.x() + (dim.x()-1) * spacing_,
                          min_.y() + (dim.y()-1) * spacing_,
                          min_.z() + (dim.z()-1) * spacing_);
  m_min = min_;
  m_max = max_;
  m_points = dim;
  m_spacing = Vector3d(spacing_, spacing_, spacing_);
  m_data.resize(m_points.x() * m_points.y() * m_points.z());
  return true;
}

void Map::Draw(Gdiplus::Graphics* g)
{
	Vector3i northWestTile = mMapViewport->GetNorthWestTileCoordinate();
	Vector3i southEastTile = mMapViewport->GetSouthEastTileCoordinate();
	Vector2i origin = mMapViewport->GetTileOrigin(northWestTile);
	int xTileCount = southEastTile.GetX() - northWestTile.GetX() + 1;
	int yTileCount = southEastTile.GetY() - northWestTile.GetY() + 1;
	int tileCount = xTileCount*yTileCount;

	for(int i = 0; i<xTileCount; i++)
	{
		for(int j = 0; j<yTileCount; j++)
		{
			Vector3i coord(northWestTile.GetX() + i, northWestTile.GetY() + j, mMapViewport->GetZoom());
			Tile* tile = GetTile(coord);
			if(!tile->IsLoaded())
			{
				tile->SignalReady += [this](Tile* tile) {
					std::lock_guard<std::mutex> lock(signal_mutex);
					SignalNewTile.emit();
				};
				continue;
			}
			Gdiplus::Image* im = tile->GetImage();
			if(im)
				g->DrawImage(im, origin.GetX() + i*mMapSource->GetTileSize(), origin.GetY() + j*mMapSource->GetTileSize());
		}
	}
}

void AbstractTexture::subImageImplementationDefault(GLenum target, GLint level, const Vector3i& offset, const Vector3i& size, AbstractImage::Format format, AbstractImage::Type type, const GLvoid* data) {
    bindInternal();
    /** @todo Get some extension wrangler instead to avoid linker errors to glTexSubImage3D() on ES2 */
    #ifndef MAGNUM_TARGET_GLES2
    glTexSubImage3D(target, level, offset.x(), offset.y(), offset.z(), size.x(), size.y(), size.z(), static_cast<GLenum>(format), static_cast<GLenum>(type), data);
    #else
    static_cast<void>(target);
    static_cast<void>(level);
    static_cast<void>(offset);
    static_cast<void>(size);
    static_cast<void>(format);
    static_cast<void>(type);
    static_cast<void>(data);
    #endif
}

MatrixX3f Surface::compute_normals(const MatrixX3f& rr, const MatrixX3i& tris)
{
    printf("\tcomputing normals\n");
    // first, compute triangle normals
    MatrixX3f r1(tris.rows(),3); MatrixX3f r2(tris.rows(),3); MatrixX3f r3(tris.rows(),3);

    for(qint32 i = 0; i < tris.rows(); ++i)
    {
        r1.row(i) = rr.row(tris(i, 0));
        r2.row(i) = rr.row(tris(i, 1));
        r3.row(i) = rr.row(tris(i, 2));
    }

    MatrixX3f x = r2 - r1;
    MatrixX3f y = r3 - r1;
    MatrixX3f tri_nn(x.rows(),y.cols());
    tri_nn.col(0) = x.col(1).cwiseProduct(y.col(2)) - x.col(2).cwiseProduct(y.col(1));
    tri_nn.col(1) = x.col(2).cwiseProduct(y.col(0)) - x.col(0).cwiseProduct(y.col(2));
    tri_nn.col(2) = x.col(0).cwiseProduct(y.col(1)) - x.col(1).cwiseProduct(y.col(0));

    //   Triangle normals and areas
    MatrixX3f tmp = tri_nn.cwiseProduct(tri_nn);
    VectorXf normSize = tmp.rowwise().sum();
    normSize = normSize.cwiseSqrt();

    for(qint32 i = 0; i < normSize.size(); ++i)
        if(normSize(i) != 0)
            tri_nn.row(i) /= normSize(i);

    MatrixX3f nn = MatrixX3f::Zero(rr.rows(), 3);

    for(qint32 p = 0; p < tris.rows(); ++p)
    {
        Vector3i verts = tris.row(p);
        for(qint32 j = 0; j < verts.size(); ++j)
            nn.row(verts(j)) = tri_nn.row(p);
    }

    tmp = nn.cwiseProduct(nn);
    normSize = tmp.rowwise().sum();
    normSize = normSize.cwiseSqrt();

    for(qint32 i = 0; i < normSize.size(); ++i)
        if(normSize(i) != 0)
            nn.row(i) /= normSize(i);

    return nn;
}

void setScalarAsFloat(T* data, Vector3i position, Vector3i size, float value, uchar channel, uchar nrOfChannels) {

    if(position.x() < 0 || position.y() < 0 || position.z() < 0 ||
            position.x() > size.x()-1 || position.y() > size.y()-1 || position.z() > size.z()-1 || channel >= nrOfChannels)
        throw OutOfBoundsException();

    data[(position.x() + position.y()*size.x() + position.z()*size.x()*size.y())*nrOfChannels + channel] = value;

}

bool Box::contains(const Vector3i& point) const {
    return point.x() >= m_llf.x() && point.y() >= m_llf.y() && point.z() >= m_llf.z() && m_urb.x() >= point.x() && m_urb.y() >= point.y() &&
           m_urb.z() >= point.z();
}

InitialParticle::InitialParticle(Vector3i position, Vector3i color, int id) : position(position), color(color), id(id) {
    streakLines = new std::list<Particle*>();
    streakLines->push_back(new Particle(position.toMyVector3D(), color, id, position.toMyVector3D()));
}

Array VoxelMesher::build(const VoxelBuffer &buffer, unsigned int channel, Vector3i min, Vector3i max) {
	uint64_t time_before = OS::get_singleton()->get_ticks_usec();

	ERR_FAIL_COND_V(_library.is_null(), Array());
	ERR_FAIL_COND_V(channel >= VoxelBuffer::MAX_CHANNELS, Array());

	const VoxelLibrary &library = **_library;

	for (unsigned int i = 0; i < MAX_MATERIALS; ++i) {
		Arrays &a = _arrays[i];
		a.positions.clear();
		a.normals.clear();
		a.uvs.clear();
		a.colors.clear();
		a.indices.clear();
	}

	float baked_occlusion_darkness;
	if (_bake_occlusion)
		baked_occlusion_darkness = _baked_occlusion_darkness / 3.0;

	// The technique is Culled faces.
	// Could be improved with greedy meshing: https://0fps.net/2012/06/30/meshing-in-a-minecraft-game/
	// However I don't feel it's worth it yet:
	// - Not so much gain for organic worlds with lots of texture variations
	// - Works well with cubes but not with any shape
	// - Slower
	// => Could be implemented in a separate class?

	// Data must be padded, hence the off-by-one
	Vector3i::sort_min_max(min, max);
	const Vector3i pad(1, 1, 1);
	min.clamp_to(pad, max);
	max.clamp_to(min, buffer.get_size() - pad);

	int index_offset = 0;

	// Iterate 3D padded data to extract voxel faces.
	// This is the most intensive job in this class, so all required data should be as fit as possible.

	// The buffer we receive MUST be dense (i.e not compressed, and channels allocated).
	// That means we can use raw pointers to voxel data inside instead of using the higher-level getters,
	// and then save a lot of time.

	uint8_t *type_buffer = buffer.get_channel_raw(Voxel::CHANNEL_TYPE);
	//       _
	//      | \
	//     /\ \\
	//    / /|\\\
	//    | |\ \\\
	//    | \_\ \\|
	//    |    |  )
	//     \   |  |
	//      \    /
	CRASH_COND(type_buffer == NULL);


	//CRASH_COND(memarr_len(type_buffer) != buffer.get_volume() * sizeof(uint8_t));

	// Build lookup tables so to speed up voxel access.
	// These are values to add to an address in order to get given neighbor.

	int row_size = buffer.get_size().y;
	int deck_size = buffer.get_size().x * row_size;

	int side_neighbor_lut[Cube::SIDE_COUNT];
	side_neighbor_lut[Cube::SIDE_LEFT] = row_size;
	side_neighbor_lut[Cube::SIDE_RIGHT] = -row_size;
	side_neighbor_lut[Cube::SIDE_BACK] = -deck_size;
	side_neighbor_lut[Cube::SIDE_FRONT] = deck_size;
	side_neighbor_lut[Cube::SIDE_BOTTOM] = -1;
	side_neighbor_lut[Cube::SIDE_TOP] = 1;

	int edge_neighbor_lut[Cube::EDGE_COUNT];
	edge_neighbor_lut[Cube::EDGE_BOTTOM_BACK] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_BACK];
	edge_neighbor_lut[Cube::EDGE_BOTTOM_FRONT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_FRONT];
	edge_neighbor_lut[Cube::EDGE_BOTTOM_LEFT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_LEFT];
	edge_neighbor_lut[Cube::EDGE_BOTTOM_RIGHT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_RIGHT];
	edge_neighbor_lut[Cube::EDGE_BACK_LEFT] = side_neighbor_lut[Cube::SIDE_BACK] + side_neighbor_lut[Cube::SIDE_LEFT];
	edge_neighbor_lut[Cube::EDGE_BACK_RIGHT] = side_neighbor_lut[Cube::SIDE_BACK] + side_neighbor_lut[Cube::SIDE_RIGHT];
	edge_neighbor_lut[Cube::EDGE_FRONT_LEFT] = side_neighbor_lut[Cube::SIDE_FRONT] + side_neighbor_lut[Cube::SIDE_LEFT];
	edge_neighbor_lut[Cube::EDGE_FRONT_RIGHT] = side_neighbor_lut[Cube::SIDE_FRONT] + side_neighbor_lut[Cube::SIDE_RIGHT];
	edge_neighbor_lut[Cube::EDGE_TOP_BACK] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_BACK];
	edge_neighbor_lut[Cube::EDGE_TOP_FRONT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_FRONT];
	edge_neighbor_lut[Cube::EDGE_TOP_LEFT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_LEFT];
	edge_neighbor_lut[Cube::EDGE_TOP_RIGHT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_RIGHT];

	int corner_neighbor_lut[Cube::CORNER_COUNT];
	corner_neighbor_lut[Cube::CORNER_BOTTOM_BACK_LEFT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_BACK] + side_neighbor_lut[Cube::SIDE_LEFT];
	corner_neighbor_lut[Cube::CORNER_BOTTOM_BACK_RIGHT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_BACK] + side_neighbor_lut[Cube::SIDE_RIGHT];
	corner_neighbor_lut[Cube::CORNER_BOTTOM_FRONT_RIGHT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_FRONT] + side_neighbor_lut[Cube::SIDE_RIGHT];
	corner_neighbor_lut[Cube::CORNER_BOTTOM_FRONT_LEFT] = side_neighbor_lut[Cube::SIDE_BOTTOM] + side_neighbor_lut[Cube::SIDE_FRONT] + side_neighbor_lut[Cube::SIDE_LEFT];
	corner_neighbor_lut[Cube::CORNER_TOP_BACK_LEFT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_BACK] + side_neighbor_lut[Cube::SIDE_LEFT];
	corner_neighbor_lut[Cube::CORNER_TOP_BACK_RIGHT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_BACK] + side_neighbor_lut[Cube::SIDE_RIGHT];
	corner_neighbor_lut[Cube::CORNER_TOP_FRONT_RIGHT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_FRONT] + side_neighbor_lut[Cube::SIDE_RIGHT];
	corner_neighbor_lut[Cube::CORNER_TOP_FRONT_LEFT] = side_neighbor_lut[Cube::SIDE_TOP] + side_neighbor_lut[Cube::SIDE_FRONT] + side_neighbor_lut[Cube::SIDE_LEFT];

	uint64_t time_prep = OS::get_singleton()->get_ticks_usec() - time_before;
	time_before = OS::get_singleton()->get_ticks_usec();

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