diff --git a/pybindings/src/bounds.cpp b/pybindings/src/bounds.cpp deleted file mode 100644 index ef97472..0000000 --- a/pybindings/src/bounds.cpp +++ /dev/null @@ -1,144 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "bounds.hpp" - -Bounds::Bounds(const WorldDimension& worldDimension, const FloatPoint& position, - const Dimension& dimension) - : _worldDimension(worldDimension), - _position(position), - _dimension(dimension) -{ - if (_dimension.getWidth() > _worldDimension.getWidth()) { - throw runtime_error("(Bounds::Bounds) Plate is larger than the world containing it"); - } - if (_dimension.getHeight() > _worldDimension.getHeight()) { - throw runtime_error("(Bounds::Bounds) Plate is taller than the world containing it"); - } -} - -uint32_t Bounds::index(uint32_t x, uint32_t y) const -{ - if (x >= _dimension.getWidth()) { - throw runtime_error("Bounds::Index: unvalid x coordinate"); - } - if (y >= _dimension.getHeight()) { - throw runtime_error("Bounds::Index: unvalid y coordinate"); - } - return y * _dimension.getWidth() + x; -} - -uint32_t Bounds::area() const -{ - return _dimension.getArea(); -} - -uint32_t Bounds::width() const -{ - return _dimension.getWidth(); -} - -uint32_t Bounds::height() const -{ - return _dimension.getHeight(); -} - -uint32_t Bounds::leftAsUint() const -{ - p_assert(_position.getX() >= 0, "_position.getX() should be not negative"); - return _position.getX(); -} - -uint32_t Bounds::topAsUint() const -{ - p_assert(_position.getY() >= 0, "_position.getY() should be not negative"); - return _position.getY(); -} - -uint32_t Bounds::rightAsUintNonInclusive() const -{ - return leftAsUint() + width() - 1; -} - -uint32_t Bounds::bottomAsUintNonInclusive() const -{ - return topAsUint() + height() - 1; -} - -bool Bounds::containsWorldPoint(uint32_t x, uint32_t y) const -{ - return asRect().contains(x, y); -} - -bool Bounds::isInLimits(float x, float y) const -{ - if (x<0) return false; - if (y<0) return false; - uint32_t ux = x; - uint32_t uy = y; - return ux < _dimension.getWidth() && uy < _dimension.getHeight(); -} - -void Bounds::shift(float dx, float dy) { - _position.shift(dx, dy, _worldDimension); - p_assert(_worldDimension.contains(_position), ""); -} - -void Bounds::grow(int dx, int dy) -{ - if (dx<0) throw runtime_error("negative value"); - if (dy<0) throw runtime_error("negative value"); - _dimension.grow(dx, dy); - - if (_dimension.getWidth() > _worldDimension.getWidth()) { - throw runtime_error("(Bounds::grow) Plate is larger than the world containing it"); - } - if (_dimension.getHeight() > _worldDimension.getHeight()) { - string s("(Bounds::grow) Plate is taller than the world containing it:"); - s += " delta=" + Platec::to_string(dy); - s += " resulting plate height=" + Platec::to_string(_dimension.getHeight()); - s += " world height=" + Platec::to_string(_worldDimension.getHeight()); - throw runtime_error(s); - } -} - -Platec::Rectangle Bounds::asRect() const -{ - p_assert(_position.getX() > 0.0f && _position.getY() >= 0.0f, "(Bounds::asRect) Left and top must be positive"); - const uint32_t ilft = leftAsUint(); - const uint32_t itop = topAsUint(); - const uint32_t irgt = ilft + _dimension.getWidth(); - const uint32_t ibtm = itop + _dimension.getHeight(); - - return Platec::Rectangle(_worldDimension, ilft, irgt, itop, ibtm); -} - -uint32_t Bounds::getMapIndex(uint32_t* px, uint32_t* py) const -{ - return asRect().getMapIndex(px, py); -} - -uint32_t Bounds::getValidMapIndex(uint32_t* px, uint32_t* py) const -{ - uint32_t res = asRect().getMapIndex(px, py); - if (res == BAD_INDEX) { - throw runtime_error("BAD INDEX found"); - } - return res; -} diff --git a/pybindings/src/geometry.cpp b/pybindings/src/geometry.cpp deleted file mode 100644 index a567dfb..0000000 --- a/pybindings/src/geometry.cpp +++ /dev/null @@ -1,272 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include -#include "geometry.hpp" - -// -// IntPoint -// - -IntPoint::IntPoint(int x, int y) - : _x(x), _y(y) -{ } - -int IntPoint::getX() const -{ - return _x; -} - -int IntPoint::getY() const -{ - return _y; -} - -IntPoint operator-(const IntPoint& a, const IntPoint& b) { - return IntPoint(a.getX() - b.getX(), a.getY() - b.getY()); -} - -// -// FloatPoint -// - -FloatPoint::FloatPoint(float x, float y) - : _x(x), _y(y) -{ } - -float FloatPoint::getX() const -{ - return _x; -} - -float FloatPoint::getY() const -{ - return _y; -} - -void FloatPoint::shift(float dx, float dy, const WorldDimension& _worldDimension) -{ - _x += dx; - _x += _x > 0 ? 0 : _worldDimension.getWidth(); - _x -= _x < _worldDimension.getWidth() ? 0 : _worldDimension.getWidth(); - - _y += dy; - _y += _y > 0 ? 0 : _worldDimension.getHeight(); - _y -= _y < _worldDimension.getHeight() ? 0 : _worldDimension.getHeight(); - - p_assert(_worldDimension.contains(*this), "(FloatPoint::shift)"); -} - -IntPoint FloatPoint::toInt() const { - return IntPoint((int)_x, (int)_y); -} - -// -// Dimension -// - -Dimension::Dimension(uint32_t width, uint32_t height) : - _width(width), _height(height) -{ -} - -Dimension::Dimension(const Dimension& original) : - _width(original.getWidth()), _height(original.getHeight()) -{ -} - -uint32_t Dimension::getWidth() const -{ - return _width; -} - -uint32_t Dimension::getHeight() const -{ - return _height; -} - -uint32_t Dimension::getArea() const -{ - return _width * _height; -} - -bool Dimension::contains(const uint32_t x, const uint32_t y) const -{ - return (x >= 0 && x < _width && y >= 0 && y < _height); -} - -bool Dimension::contains(const float x, const float y) const -{ - return (x >= 0 && x < _width && y >= 0 && y < _height); -} - -bool Dimension::contains(const FloatPoint& p) const -{ - return (p.getX() >= 0 && p.getX() < _width && p.getY() >= 0 && p.getY() < _height); -} - -void Dimension::grow(uint32_t amountX, uint32_t amountY) -{ - _width += amountX; - _height += amountY; -} - -// -// WorldDimension -// - -WorldDimension::WorldDimension(uint32_t width, uint32_t height) : Dimension(width, height) -{ -}; - -WorldDimension::WorldDimension(const WorldDimension& original) : Dimension(original) -{ -}; - -uint32_t WorldDimension::getMax() const -{ - return _width > _height ? _width : _height; -} - -uint32_t WorldDimension::xMod(int x) const -{ - return (x + _width) % getWidth(); -} - -uint32_t WorldDimension::yMod(int y) const -{ - return (y + _height) % getHeight(); -} - -uint32_t WorldDimension::xMod(uint32_t x) const -{ - return (x + _width) % getWidth(); -} - -uint32_t WorldDimension::yMod(uint32_t y) const -{ - return (y + _height) % getHeight(); -} - -void WorldDimension::normalize(uint32_t& x, uint32_t& y) const -{ - x %= _width; - y %= _height; -} - -uint32_t WorldDimension::indexOf(const uint32_t x, const uint32_t y) const -{ - return y * getWidth() + x; -} - -uint32_t WorldDimension::lineIndex(const uint32_t y) const -{ - if (y<0 || y>=_height){ - throw invalid_argument("WorldDimension::line_index: y is not valid"); - } - return indexOf(0, y); -} - -uint32_t WorldDimension::yFromIndex(const uint32_t index) const -{ - return index / _width; -} - -uint32_t WorldDimension::xFromIndex(const uint32_t index) const -{ - const uint32_t y = yFromIndex(index); - return index - y * _width; -} - -uint32_t WorldDimension::normalizedIndexOf(const uint32_t x, const uint32_t y) const -{ - return indexOf(xMod(x), yMod(y)); -} - -uint32_t WorldDimension::xCap(const uint32_t x) const -{ - return x < _width ? x : (_width-1); -} - -uint32_t WorldDimension::yCap(const uint32_t y) const -{ - return y < _height ? y : (_height-1); -} - -uint32_t WorldDimension::largerSize() const -{ - return _width > _height ? _width : _height; -} - - -namespace Platec { - -IntVector::IntVector(int x, int y) - : _x(x), _y(y) -{ } - -float IntVector::length() const { - float nx = _x; - float ny = _y; - return sqrt(nx * nx + ny * ny); -} - -IntVector IntVector::fromDistance(const IntPoint& a, const IntPoint& b){ - return IntVector(a.getX() - b.getX(), a.getY() - b.getY()); -} - -IntVector operator-(const IntVector& a, const IntVector& b) { - return IntVector(a.x() - b.x(), a.y() - b.y()); -} - -FloatVector IntVector::toUnitVector() const { - return FloatVector(_x/length(), _y/length()); -} - -bool operator==(const FloatVector& a, const FloatVector& b) { - return a.x() == b.x() && a.y() == b.y(); -} - -FloatVector::FloatVector(float x, float y) - : _x(x), _y(y) -{ } - -IntVector FloatVector::toIntVector() const { - return IntVector((int)_x, (int)_y); -} - -FloatVector operator-(const FloatVector& a, const FloatVector& b) { - return FloatVector(a.x() - b.x(), a.y() - b.y()); -} - -FloatVector operator*(const FloatVector& v, float f) { - return FloatVector(v.x() * f, v.y() * f); -} - -float FloatVector::dotProduct(const FloatVector& other) const { - return x() * other.x() + y() * other.y(); -} - -float FloatVector::length() const { - float nx = _x; - float ny = _y; - return sqrt(nx * nx + ny * ny); -} - -} diff --git a/pybindings/src/heightmap.cpp b/pybindings/src/heightmap.cpp deleted file mode 100644 index a7a0c52..0000000 --- a/pybindings/src/heightmap.cpp +++ /dev/null @@ -1,24 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "heightmap.hpp" - -#include // std::invalid_argument - -using namespace std; diff --git a/pybindings/src/lithosphere.cpp b/pybindings/src/lithosphere.cpp deleted file mode 100644 index 86e822f..0000000 --- a/pybindings/src/lithosphere.cpp +++ /dev/null @@ -1,879 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "lithosphere.hpp" -#include "plate.hpp" -#include "sqrdmd.hpp" -#include "simplexnoise.hpp" -#include "noise.hpp" - -#include -#include -#include -#include -#include -#include -#include - -#define BOOL_REGENERATE_CRUST 1 - -using namespace std; - -/** - * Wrapper for growing plate from a seed. Contains plate's dimensions. - * - * Used exclusively in plate creation phase. - */ -class plateArea -{ - public: - vector border; ///< Plate's unprocessed border pixels. - uint32_t btm; ///< Most bottom pixel of plate. - uint32_t lft; ///< Most left pixel of plate. - uint32_t rgt; ///< Most right pixel of plate. - uint32_t top; ///< Most top pixel of plate. - uint32_t wdt; ///< Width of area in pixels. - uint32_t hgt; ///< Height of area in pixels. -}; - -static const float SUBDUCT_RATIO = 0.5f; - -static const float BUOYANCY_BONUS_X = 3; -static const uint32_t MAX_BUOYANCY_AGE = 20; -static const float MULINV_MAX_BUOYANCY_AGE = 1.0f / (float)MAX_BUOYANCY_AGE; - -static const float RESTART_ENERGY_RATIO = 0.15; -static const float RESTART_SPEED_LIMIT = 2.0; -static const uint32_t NO_COLLISION_TIME_LIMIT = 10; - -uint32_t findBound(const uint32_t* map, uint32_t length, uint32_t x0, uint32_t y0, - int dx, int dy); -uint32_t findPlate(plate** plates, float x, float y, uint32_t num_plates); - -WorldPoint lithosphere::randomPosition() -{ - return WorldPoint( - _randsource.next() % _worldDimension.getWidth(), - _randsource.next() % _worldDimension.getHeight(), - _worldDimension); -} - -void lithosphere::createNoise(float* tmp, const WorldDimension& tmpDim, bool useSimplex) -{ - ::createNoise(tmp, tmpDim, _randsource, useSimplex); -} - -void lithosphere::createSlowNoise(float* tmp, const WorldDimension& tmpDim) -{ - ::createSlowNoise(tmp, tmpDim, _randsource); -} - -lithosphere::lithosphere(long seed, uint32_t width, uint32_t height, float sea_level, - uint32_t _erosion_period, float _folding_ratio, uint32_t aggr_ratio_abs, - float aggr_ratio_rel, uint32_t num_cycles) throw(invalid_argument) : - hmap(width, height), - amap(width, height), - imap(width, height), - plates(0), - aggr_overlap_abs(aggr_ratio_abs), - aggr_overlap_rel(aggr_ratio_rel), - cycle_count(0), - erosion_period(_erosion_period), - folding_ratio(_folding_ratio), - iter_count(0), - max_cycles(num_cycles), - max_plates(0), - num_plates(0), - _worldDimension(width, height), - _randsource(seed), - _steps(0) -{ - if (width < 5 || height < 5){ - throw runtime_error("Width and height should be >=5"); - } - - WorldDimension tmpDim = WorldDimension(width+1, height+1); - const uint32_t A = tmpDim.getArea(); - float* tmp = new float[A]; - memset(tmp, 0, A * sizeof(float)); - - createSlowNoise(tmp, tmpDim); - - float lowest = tmp[0], highest = tmp[0]; - for (uint32_t i = 1; i < A; ++i) - { - lowest = lowest < tmp[i] ? lowest : tmp[i]; - highest = highest > tmp[i] ? highest : tmp[i]; - } - - for (uint32_t i = 0; i < A; ++i) // Scale to [0 ... 1] - tmp[i] = (tmp[i] - lowest) / (highest - lowest); - - float sea_threshold = 0.5; - float th_step = 0.5; - - // Find the actual value in height map that produces the continent-sea - // ratio defined be "sea_level". - while (th_step > 0.01) - { - uint32_t count = 0; - for (uint32_t i = 0; i < A; ++i) - count += (tmp[i] < sea_threshold); - - th_step *= 0.5; - if (count / (float)A < sea_level) - sea_threshold += th_step; - else - sea_threshold -= th_step; - } - - sea_level = sea_threshold; - for (uint32_t i = 0; i < A; ++i) // Genesis 1:9-10. - { - tmp[i] = (tmp[i] > sea_level) * - (tmp[i] + CONTINENTAL_BASE) + - (tmp[i] <= sea_level) * OCEANIC_BASE; - } - - // Scalp the +1 away from map side to get a power of two side length! - // Practically only the redundant map edges become removed. - for (uint32_t y = 0; y < _worldDimension.getHeight(); ++y) { - memcpy(&hmap[_worldDimension.lineIndex(y)], - &tmp[ tmpDim.lineIndex(y)], - _worldDimension.getWidth()*sizeof(float)); - } - - delete[] tmp; -} - -lithosphere::~lithosphere() throw() -{ - delete[] plates; plates = 0; -} - -void lithosphere::growPlates(plateArea*& area, IndexMap& owner) -{ - // "Grow" plates from their origins until surface is fully populated. - uint32_t max_border = 1; - uint32_t i; - while (max_border) { - for (max_border = i = 0; i < num_plates; ++i) { - const uint32_t N = area[i].border.size(); - max_border = max_border > N ? max_border : N; - - if (N == 0) - continue; - - const uint32_t j = _randsource.next() % N; - const uint32_t p = area[i].border[j]; - const uint32_t cy = _worldDimension.yFromIndex(p); - const uint32_t cx = _worldDimension.xFromIndex(p); - - const uint32_t lft = cx > 0 ? cx - 1 : _worldDimension.getWidth() - 1; - const uint32_t rgt = cx < _worldDimension.getWidth() - 1 ? cx + 1 : 0; - const uint32_t top = cy > 0 ? cy - 1 : _worldDimension.getHeight() - 1; - const uint32_t btm = cy < _worldDimension.getHeight() - 1 ? cy + 1 : 0; - - const uint32_t n = top * _worldDimension.getWidth() + cx; // North. - const uint32_t s = btm * _worldDimension.getWidth() + cx; // South. - const uint32_t w = cy * _worldDimension.getWidth() + lft; // West. - const uint32_t e = cy * _worldDimension.getWidth() + rgt; // East. - - if (owner[n] >= num_plates) - { - owner[n] = i; - area[i].border.push_back(n); - - if (area[i].top == _worldDimension.yMod(top + 1)) - { - area[i].top = top; - area[i].hgt++; - } - } - - if (owner[s] >= num_plates) - { - owner[s] = i; - area[i].border.push_back(s); - - if (btm == _worldDimension.yMod(area[i].btm + 1)) - { - area[i].btm = btm; - area[i].hgt++; - } - } - - if (owner[w] >= num_plates) - { - owner[w] = i; - area[i].border.push_back(w); - - if (area[i].lft == _worldDimension.xMod(lft + 1)) - { - area[i].lft = lft; - area[i].wdt++; - } - } - - if (owner[e] >= num_plates) - { - owner[e] = i; - area[i].border.push_back(e); - - if (rgt == _worldDimension.xMod(area[i].rgt + 1)) - { - area[i].rgt = rgt; - area[i].wdt++; - } - } - - // Overwrite processed point with unprocessed one. - area[i].border[j] = area[i].border.back(); - area[i].border.pop_back(); - } - } -} - -void lithosphere::createPlates(uint32_t num_plates) -{ -try { - const uint32_t map_area = _worldDimension.getArea(); - this->max_plates = this->num_plates = num_plates; - - std::vector vec; - vec.reserve(_worldDimension.largerSize()*4); // == map's circumference. - - collisions.reserve(num_plates); - subductions.reserve(num_plates); - - for (uint32_t i = 0; i < num_plates; ++i) - { - collisions.push_back(vec); - subductions.push_back(vec); - } - - // Initialize "Free plate center position" lookup table. - // This way two plate centers will never be identical. - for (uint32_t i = 0; i < map_area; ++i) - imap[i] = i; - - // Select N plate centers from the global map. - plateArea* area = new plateArea[num_plates]; - for (uint32_t i = 0; i < num_plates; ++i) - { - // Randomly select an unused plate origin. - const uint32_t p = imap[(uint32_t)_randsource.next() % (map_area - i)]; - const uint32_t y = _worldDimension.yFromIndex(p); - const uint32_t x = _worldDimension.xFromIndex(p); - - area[i].lft = area[i].rgt = x; // Save origin... - area[i].top = area[i].btm = y; - area[i].wdt = area[i].hgt = 1; - - area[i].border.reserve(8); - area[i].border.push_back(p); // ...and mark it as border. - - // Overwrite used entry with last unused entry in array. - imap[p] = imap[map_area - i - 1]; - } - - imap.set_all(0xFFFFFFFF); - - growPlates(area, imap); - - // check all the points of the map are owned - for (int i=0; i < imap.area(); i++){ - p_assert(imap[i]addCollision(x_mod, y_mod); - uint32_t prev_area = plates[imap[k]]->addCollision(x_mod, y_mod); - - if (this_area < prev_area) - { - plateCollision coll(imap[k], x_mod, y_mod, - this_map[j] * folding_ratio); - - // Give some... - hmap[k] += coll.crust; - plates[imap[k]]->setCrust(x_mod, y_mod, hmap[k], - this_age[j]); - - // And take some. - plates[i]->setCrust(x_mod, y_mod, this_map[j] * - (1.0 - folding_ratio), this_age[j]); - - // Add collision to the earlier plate's list. - collisions[i].push_back(coll); - ++continental_collisions; - } - else - { - plateCollision coll(i, x_mod, y_mod, - hmap[k] * folding_ratio); - - plates[i]->setCrust(x_mod, y_mod, - this_map[j]+coll.crust, amap[k]); - - plates[imap[k]]->setCrust(x_mod, y_mod, hmap[k] - * (1.0 - folding_ratio), amap[k]); - - collisions[imap[k]].push_back(coll); - ++continental_collisions; - - // Give the location to the larger plate. - hmap[k] = this_map[j]; - imap[k] = i; - amap[k] = this_age[j]; - } -} - -// Update height and plate index maps. -// Doing it plate by plate is much faster than doing it index wise: -// Each plate's map's memory area is accessed sequentially and only -// once as opposed to calculating "num_plates" indices within plate -// maps in order to find out which plate(s) own current location. -void lithosphere::updateHeightAndPlateIndexMaps(const uint32_t& map_area, - uint32_t& oceanic_collisions, - uint32_t& continental_collisions) -{ - hmap.set_all(0); - imap.set_all(0xFFFFFFFF); - for (uint32_t i = 0; i < num_plates; ++i) - { - const uint32_t x0 = plates[i]->getLeftAsUint(); - const uint32_t y0 = plates[i]->getTopAsUint(); - const uint32_t x1 = x0 + plates[i]->getWidth(); - const uint32_t y1 = y0 + plates[i]->getHeight(); - - const float* this_map; - const uint32_t* this_age; - plates[i]->getMap(&this_map, &this_age); - - // Copy first part of plate onto world map. - for (uint32_t y = y0, j = 0; y < y1; ++y) - { - for (uint32_t x = x0; x < x1; ++x, ++j) - { - const uint32_t x_mod = _worldDimension.xMod(x); - const uint32_t y_mod = _worldDimension.yMod(y); - - const uint32_t k = _worldDimension.indexOf(x_mod, y_mod); - - if (this_map[j] < 2 * FLT_EPSILON) // No crust here... - continue; - - if (imap[k] >= num_plates) // No one here yet? - { - // This plate becomes the "owner" of current location - // if it is the first plate to have crust on it. - hmap[k] = this_map[j]; - imap[k] = i; - amap[k] = this_age[j]; - - continue; - } - - // DO NOT ACCEPT HEIGHT EQUALITY! Equality leads to subduction - // of shore that 's barely above sea level. It's a lot less - // serious problem to treat very shallow waters as continent... - const bool prev_is_oceanic = hmap[k] < CONTINENTAL_BASE; - const bool this_is_oceanic = this_map[j] < CONTINENTAL_BASE; - - const uint32_t prev_timestamp = plates[imap[k]]-> - getCrustTimestamp(x_mod, y_mod); - const uint32_t this_timestamp = this_age[j]; - const uint32_t prev_is_bouyant = (hmap[k] > this_map[j]) | - ((hmap[k] + 2 * FLT_EPSILON > this_map[j]) & - (hmap[k] < 2 * FLT_EPSILON + this_map[j]) & - (prev_timestamp >= this_timestamp)); - - // Handle subduction of oceanic crust as special case. - if (this_is_oceanic & prev_is_bouyant) - { - // This plate will be the subducting one. - // The level of effect that subduction has - // is directly related to the amount of water - // on top of the subducting plate. - const float sediment = SUBDUCT_RATIO * OCEANIC_BASE * - (CONTINENTAL_BASE - this_map[j]) / - CONTINENTAL_BASE; - - // Save collision to the receiving plate's list. - plateCollision coll(i, x_mod, y_mod, sediment); - subductions[imap[k]].push_back(coll); - ++oceanic_collisions; - - // Remove subducted oceanic lithosphere from plate. - // This is crucial for - // a) having correct amount of colliding crust (below) - // b) protecting subducted locations from receiving - // crust from other subductions/collisions. - plates[i]->setCrust(x_mod, y_mod, this_map[j] - - OCEANIC_BASE, this_timestamp); - - if (this_map[j] <= 0) - continue; // Nothing more to collide. - } - else if (prev_is_oceanic) - { - const float sediment = SUBDUCT_RATIO * OCEANIC_BASE * - (CONTINENTAL_BASE - hmap[k]) / - CONTINENTAL_BASE; - - plateCollision coll(imap[k], x_mod, y_mod, sediment); - subductions[i].push_back(coll); - ++oceanic_collisions; - - plates[imap[k]]->setCrust(x_mod, y_mod, hmap[k] - - OCEANIC_BASE, prev_timestamp); - hmap[k] -= OCEANIC_BASE; - - if (hmap[k] <= 0) - { - imap[k] = i; - hmap[k] = this_map[j]; - amap[k] = this_age[j]; - - continue; - } - } - - resolveJuxtapositions(i, j, k, x_mod, y_mod, - this_map, this_age, continental_collisions); - } - } - } -} - -void lithosphere::updateCollisions() -{ - for (uint32_t i = 0; i < num_plates; ++i) - { - for (uint32_t j = 0; j < collisions[i].size(); ++j) - { - const plateCollision& coll = collisions[i][j]; - uint32_t coll_count, coll_count_i, coll_count_j; - float coll_ratio, coll_ratio_i, coll_ratio_j; - - #ifdef DEBUG - if (i == coll.index) - { - puts("when colliding: SRC == DEST!"); - exit(1); - } - #endif - - // Collision causes friction. Apply it to both plates. - plates[i]->applyFriction(coll.crust); - plates[coll.index]->applyFriction(coll.crust); - - plates[i]->getCollisionInfo(coll.wx, coll.wy, - &coll_count_i, &coll_ratio_i); - plates[coll.index]->getCollisionInfo(coll.wx, - coll.wy, &coll_count_j, &coll_ratio_j); - - // Find the minimum count of collisions between two - // continents on different plates. - // It's minimum because large plate will get collisions - // from all over where as smaller plate will get just - // a few. It's those few that matter between these two - // plates, not what the big plate has with all the - // other plates around it. - coll_count = coll_count_i; - coll_count -= (coll_count - coll_count_j) & - -(coll_count > coll_count_j); - - // Find maximum amount of collided surface area between - // two continents on different plates. - // Like earlier, it's the "experience" of the smaller - // plate that matters here. - coll_ratio = coll_ratio_i; - coll_ratio += (coll_ratio_j - coll_ratio) * - (coll_ratio_j > coll_ratio); - - if ((coll_count > aggr_overlap_abs) | - (coll_ratio > aggr_overlap_rel)) - { - float amount = plates[i]->aggregateCrust( - plates[coll.index], - coll.wx, coll.wy); - - // Calculate new direction and speed for the - // merged plate system, that is, for the - // receiving plate! - plates[coll.index]->collide(*plates[i], - coll.wx, coll.wy, amount); - } - } - - collisions[i].clear(); - } -} - -// Remove empty plates from the system. -void lithosphere::removeEmptyPlates(uint32_t*& indexFound) -{ - for (uint32_t i = 0; i < num_plates; ++i) - { - if (num_plates == 1) - puts("ONLY ONE PLATE LEFT!"); - else if (indexFound[i] == 0) - { - delete plates[i]; - plates[i] = plates[num_plates - 1]; - indexFound[i] = indexFound[num_plates - 1]; - - // Life is seldom as simple as seems at first. - // Replace the moved plate's index in the index map - // to match its current position in the array! - for (uint32_t j = 0; j < _worldDimension.getArea(); ++j) - if (imap[j] == num_plates - 1) - imap[j] = i; - - --num_plates; - --i; - } - } -} - -void lithosphere::update() -{ -try { - _steps++; - float totalVelocity = 0; - float systemKineticEnergy = 0; - - for (uint32_t i = 0; i < num_plates; ++i) - { - totalVelocity += plates[i]->getVelocity(); - systemKineticEnergy += plates[i]->getMomentum(); - } - - if (systemKineticEnergy > peak_Ek) { - peak_Ek = systemKineticEnergy; - } - - // If there's no continental collisions during past iterations, - // then interesting activity has ceased and we should restart. - // Also if the simulation has been going on for too long already, - // restart, because interesting stuff has most likely ended. - if (totalVelocity < RESTART_SPEED_LIMIT || - systemKineticEnergy / peak_Ek < RESTART_ENERGY_RATIO || - last_coll_count > NO_COLLISION_TIME_LIMIT || - iter_count > 600) - { - restart(); - return; - } - - const uint32_t map_area = _worldDimension.getArea(); - const IndexMap prev_imap = imap; - imap = IndexMap(_worldDimension.getWidth(), _worldDimension.getHeight()); - - // Realize accumulated external forces to each plate. - for (uint32_t i = 0; i < num_plates; ++i) - { - plates[i]->resetSegments(); - - if (erosion_period > 0 && iter_count % erosion_period == 0) - plates[i]->erode(CONTINENTAL_BASE); - - plates[i]->move(); - } - - uint32_t oceanic_collisions = 0; - uint32_t continental_collisions = 0; - - updateHeightAndPlateIndexMaps(map_area, oceanic_collisions, continental_collisions); - - // Update the counter of iterations since last continental collision. - last_coll_count = (last_coll_count + 1) & - -(continental_collisions == 0); - - for (uint32_t i = 0; i < num_plates; ++i) - { - for (uint32_t j = 0; j < subductions[i].size(); ++j) - { - const plateCollision& coll = subductions[i][j]; - - p_assert(i != coll.index, "when subducting: SRC == DEST!"); - - // Do not apply friction to oceanic plates. - // This is a very cheap way to emulate slab pull. - // Just perform subduction and on our way we go! - plates[i]->addCrustBySubduction( - coll.wx, coll.wy, coll.crust, iter_count, - plates[coll.index]->getVelX(), - plates[coll.index]->getVelY()); - } - - subductions[i].clear(); - } - - updateCollisions(); - - uint32_t* indexFound = new uint32_t[num_plates]; - memset(indexFound, 0, sizeof(uint32_t)*num_plates); - - // Fill divergent boundaries with new crustal material, molten magma. - for (uint32_t y = 0, i = 0; y < BOOL_REGENERATE_CRUST * _worldDimension.getHeight(); ++y) { - for (uint32_t x = 0; x < _worldDimension.getWidth(); ++x, ++i) { - if (imap[i] >= num_plates) { - // The owner of this new crust is that neighbour plate - // who was located at this point before plates moved. - imap[i] = prev_imap[i]; - - // If this is oceanic crust then add buoyancy to it. - // Magma that has just crystallized into oceanic crust - // is more buoyant than that which has had a lot of - // time to cool down and become more dense. - amap[i] = iter_count; - hmap[i] = OCEANIC_BASE * BUOYANCY_BONUS_X; - - // This should probably not happen - if (imap[i] < num_plates) { - plates[imap[i]]->setCrust(x, y, OCEANIC_BASE, - iter_count); - } - - } else if (++indexFound[imap[i]] && hmap[i] <= 0){ - puts("Occupied point has no land mass!"); - exit(1); - } - } - } - - removeEmptyPlates(indexFound); - - delete[] indexFound; - - // Add some "virginity buoyancy" to all pixels for a visual boost! :) - for (uint32_t i = 0; i < (BUOYANCY_BONUS_X > 0) * map_area; ++i) - { - // Calculate the inverted age of this piece of crust. - // Force result to be minimum between inv. age and - // max buoyancy bonus age. - uint32_t crust_age = iter_count - amap[i]; - crust_age = MAX_BUOYANCY_AGE - crust_age; - crust_age &= -(crust_age <= MAX_BUOYANCY_AGE); - - hmap[i] += (hmap[i] < CONTINENTAL_BASE) * BUOYANCY_BONUS_X * - OCEANIC_BASE * crust_age * MULINV_MAX_BUOYANCY_AGE; - } - - ++iter_count; -} catch (const exception& e){ - std::string msg = "Problem during update: "; - msg = msg + e.what(); - cerr << msg << endl; - throw runtime_error(msg.c_str()); -} -} - -void lithosphere::restart() -{ -try { - - const uint32_t map_area = _worldDimension.getArea(); - - cycle_count += max_cycles > 0; // No increment if running for ever. - if (cycle_count > max_cycles) - return; - - // Update height map to include all recent changes. - hmap.set_all(0); - for (uint32_t i = 0; i < num_plates; ++i) - { - const uint32_t x0 = plates[i]->getLeftAsUint(); - const uint32_t y0 = plates[i]->getTopAsUint(); - const uint32_t x1 = x0 + plates[i]->getWidth(); - const uint32_t y1 = y0 + plates[i]->getHeight(); - - const float* this_map; - const uint32_t* this_age; - plates[i]->getMap(&this_map, &this_age); - - // Copy first part of plate onto world map. - for (uint32_t y = y0, j = 0; y < y1; ++y) - { - for (uint32_t x = x0; x < x1; ++x, ++j) - { - const uint32_t x_mod = _worldDimension.xMod(x); - const uint32_t y_mod = _worldDimension.yMod(y); - const float h0 = hmap[_worldDimension.indexOf(x_mod, y_mod)]; - const float h1 = this_map[j]; - const uint32_t a0 = amap[_worldDimension.indexOf(x_mod, y_mod)]; - const uint32_t a1 = this_age[j]; - - amap[_worldDimension.indexOf(x_mod, y_mod)] = (h0 *a0 +h1 *a1) /(h0 +h1); - hmap[_worldDimension.indexOf(x_mod, y_mod)] += this_map[j]; - } - } - } - - // Delete plates. - delete[] plates; - plates = 0; - num_plates = 0; - - // create new plates IFF there are cycles left to run! - // However, if max cycle count is "ETERNITY", then 0 < 0 + 1 always. - if (cycle_count < max_cycles + !max_cycles) - { - createPlates(num_plates = max_plates); - - // Restore the ages of plates' points of crust! - for (uint32_t i = 0; i < num_plates; ++i) - { - const uint32_t x0 = plates[i]->getLeftAsUint(); - const uint32_t y0 = plates[i]->getTopAsUint(); - const uint32_t x1 = x0 + plates[i]->getWidth(); - const uint32_t y1 = y0 + plates[i]->getHeight(); - - const float* this_map; - const uint32_t* this_age_const; - uint32_t* this_age; - - plates[i]->getMap(&this_map, &this_age_const); - this_age = (uint32_t *)this_age_const; - - for (uint32_t y = y0, j = 0; y < y1; ++y) - { - for (uint32_t x = x0; x < x1; ++x, ++j) - { - const uint32_t x_mod = _worldDimension.xMod(x); - const uint32_t y_mod = _worldDimension.yMod(y); - - this_age[j] = amap[_worldDimension.indexOf(x_mod, y_mod)]; - } - } - } - - return; - } - - // Add some "virginity buoyancy" to all pixels for a visual boost. - for (uint32_t i = 0; i < (BUOYANCY_BONUS_X > 0) * map_area; ++i) - { - uint32_t crust_age = iter_count - amap[i]; - crust_age = MAX_BUOYANCY_AGE - crust_age; - crust_age &= -(crust_age <= MAX_BUOYANCY_AGE); - - hmap[i] += (hmap[i] < CONTINENTAL_BASE) * BUOYANCY_BONUS_X * - OCEANIC_BASE * crust_age * MULINV_MAX_BUOYANCY_AGE; - } -} catch (const exception& e){ - std::string msg = "Problem during restart: "; - msg = msg + e.what(); - throw runtime_error(msg.c_str()); -} -} - -uint32_t lithosphere::getWidth() const -{ - return _worldDimension.getWidth(); -} - -uint32_t lithosphere::getHeight() const -{ - return _worldDimension.getHeight(); -} - -uint32_t* lithosphere::getPlatesMap() const throw() -{ - return imap.raw_data(); -} - -const plate* lithosphere::getPlate(uint32_t index) const -{ - if (index >= num_plates){ - throw runtime_error("(lithosphere::getPlate) invalid plate index"); - } - return plates[index]; -} diff --git a/pybindings/src/mass.cpp b/pybindings/src/mass.cpp deleted file mode 100644 index 2c5231c..0000000 --- a/pybindings/src/mass.cpp +++ /dev/null @@ -1,111 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - - #include "mass.hpp" - -// ---------------------------------------------- -// MassBuilder -// ---------------------------------------------- - -MassBuilder::MassBuilder(const float* m, const Dimension& dimension) - : mass(0), cx(0), cy(0) -{ - uint32_t k; - for (uint32_t y = k = 0; y < dimension.getHeight(); ++y) { - for (uint32_t x = 0; x < dimension.getWidth(); ++x, ++k) { - if (m[k] < 0.0f) { - throw runtime_error(string("(MassBuilder::MassBuilder) Crust should be not negative") - + ", m[k] " + Platec::to_string_f(m[k])); - } - addPoint(x, y, m[k]); - } - } -} - -MassBuilder::MassBuilder() - : mass(0), cx(0), cy(0) -{ - -} - -void MassBuilder::addPoint(uint32_t x, uint32_t y, float crust) -{ - if (crust < 0.0f) { - throw runtime_error(string("(MassBuilder::addPoint) Crust should be not negative") - + ", crust " + Platec::to_string_f(crust)); - } - mass += crust; - // Update the center coordinates weighted by mass. - cx += x * crust; - cy += y * crust; -} - -Mass MassBuilder::build() -{ - if (mass <= 0) { - return Mass(0, 0, 0); - } else { - return Mass(mass, cx/mass, cy/mass); - } -} - -// ---------------------------------------------- -// Mass -// ---------------------------------------------- - -Mass::Mass(float mass_, float cx_, float cy_) - : mass(mass_), cx(cx_), cy(cy_), _totalX(0), _totalY(0) -{ - -} - -void Mass::incMass(float delta) -{ - mass += delta; - if (mass < 0.0f && mass > -0.01f) { - mass = 0.0f; - } - if (mass < 0.0f) { - throw runtime_error(string("(Mass::incMass) a negative delta is not allowed: ") - + Platec::to_string_f(delta) + ", resulting mass: " + Platec::to_string_f(mass)); - } -} - -float Mass::getMass() const -{ - return mass; -} - -float Mass::getCx() const -{ - if (null()) throw runtime_error("(Mass::getCx)"); - return cx; -} - -float Mass::getCy() const -{ - if (null()) throw runtime_error("(Mass::getCy)"); - return cy; -} - -bool Mass::null() const -{ - return mass <= 0; -} - diff --git a/pybindings/src/movement.cpp b/pybindings/src/movement.cpp deleted file mode 100644 index e11be09..0000000 --- a/pybindings/src/movement.cpp +++ /dev/null @@ -1,211 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "movement.hpp" -#include "plate.hpp" -#include "mass.hpp" - -// Missing on Windows -#ifndef M_PI -#define M_PI 3.141592654 -#endif - -Movement::Movement(SimpleRandom randsource, const WorldDimension& worldDimension) - : _randsource(randsource), - velocity(1), - rot_dir(randsource.next() % 2 ? 1 : -1), - dx(0), dy(0), - _worldDimension(worldDimension) -{ - const double angle = 2 * M_PI * _randsource.next_double(); - vx = cos(angle) * INITIAL_SPEED_X; - vy = sin(angle) * INITIAL_SPEED_X; -} - -void Movement::applyFriction(float deformed_mass, float mass) -{ - if (deformed_mass < 0 || mass < 0) { - throw runtime_error("(Movement::applyFriction) negative masses make not sense"); - } - if (0.0f == mass) { - velocity = 0; - return; - } - float vel_dec = DEFORMATION_WEIGHT * deformed_mass / mass; - vel_dec = vel_dec < velocity ? vel_dec : velocity; - - // Altering the source variable causes the order of calls to - // this function to have difference when it shouldn't! - // However, it's a hack well worth the outcome. :) - velocity -= vel_dec; -} - -void Movement::move() -{ - float len; - - // Apply any new impulses to the plate's trajectory. - vx += dx; - vy += dy; - dx = 0; - dy = 0; - - // Force direction of plate to be unit vector. - // Update velocity so that the distance of movement doesn't change. - len = sqrt(vx*vx+vy*vy); - vx /= len; - vy /= len; - velocity += len - 1.0; - velocity *= velocity > 0; // Round negative values to zero. - - // Apply some circular motion to the plate. - // Force the radius of the circle to remain fixed by adjusting - // angular velocity (which depends on plate's velocity). - uint32_t world_avg_side = (_worldDimension.getWidth() + _worldDimension.getHeight()) / 2; - float alpha = rot_dir * velocity / (world_avg_side * 0.33); - float _cos = cos(alpha * velocity); - float _sin = sin(alpha * velocity); - float _vx = vx * _cos - vy * _sin; - float _vy = vy * _cos + vx * _sin; - vx = _vx; - vy = _vy; -} - -float Movement::velocityOnX() const -{ - return vx * velocity; -} - -float Movement::velocityOnY() const -{ - return vy * velocity; -} - -float Movement::velocityOnX(float length) const -{ - if (length < 0) { - throw runtime_error("(Movement::velocityOnX) negative length makes not sense"); - } - return vx * length; -} - -float Movement::velocityOnY(float length) const -{ - if (length < 0) { - throw runtime_error("(Movement::velocityOnY) negative length makes not sense"); - } - return vy * length; -} - -float Movement::dot(float dx_, float dy_) const -{ - return vx * dx_ + vy * dy_; -} - -float Movement::relativeUnitVelocityOnX(float otherVx) const -{ - return vx - otherVx; -} - -float Movement::relativeUnitVelocityOnY(float otherVy) const -{ - return vy - otherVy; -} - -float Movement::momentum(const Mass& mass) const throw() { - return mass.getMass() * velocity; -} - -void Movement::collide(const IMass& thisMass, - IPlate& otherPlate, - uint32_t wx, uint32_t wy, float coll_mass) -{ - const float coeff_rest = 0.0; // Coefficient of restitution. - // 1 = fully elastic, 0 = stick together. - Platec::IntVector massCentersDistance = Platec::IntVector::fromDistance( - otherPlate.massCenter().toInt(), thisMass.massCenter().toInt()); - if (massCentersDistance.length() <= 0) { - return; // Avoid division by zero! - } - - // Scaling is required at last when impulses are added to plates! - // Compute relative velocity between plates at the collision point. - // Because torque is not included, calc simplifies to v_ab = v_a - v_b. - Platec::FloatVector collisionDirection = massCentersDistance.toUnitVector(); - Platec::FloatVector relativeVelocity = velocityUnitVector() - otherPlate.velocityUnitVector(); - - // Get the dot product of relative velocity vector and collision vector. - // Then get the projection of v_ab along collision vector. - // Note that vector n must be a unit vector! - const float rel_dot_n = collisionDirection.dotProduct(relativeVelocity); - if (rel_dot_n <= 0) { - return; // Exit if objects are moving away from each other. - } - - // Calculate the denominator of impulse: n . n * (1 / m_1 + 1 / m_2). - // Use the mass of the colliding crust for the "donator" plate. - float denom = collisionDirection.length() * collisionDirection.length() * (1.0/otherPlate.getMass() + 1.0/coll_mass); - - // Calculate force of impulse. - float J = -(1 + coeff_rest) * rel_dot_n / denom; - - // Compute final change of trajectory. - // The plate that is the "giver" of the impulse should receive a - // force according to its pre-collision mass, not the current mass! - dx += (collisionDirection * (J / thisMass.getMass())).x(); - dy += (collisionDirection * (J / thisMass.getMass())).y(); - otherPlate.decImpulse( collisionDirection * (J / (coll_mass + otherPlate.getMass())) ); - - // In order to prove that the code above works correctly, here is an - // example calculation with ball A (mass 10) moving right at velocity - // 1 and ball B (mass 100) moving up at velocity 1. Collision point - // is at rightmost point of ball A and leftmost point of ball B. - // Radius of both balls is 2. - // ap_dx = 2; - // ap_dy = 0; - // bp_dx = -2; - // bp_dy = 0; - // nx = 2 - -2 = 4; - // ny = 0 - 0 = 0; - // n_len = sqrt(4 * 4 + 0) = 4; - // nx = 4 / 4 = 1; - // ny = 0 / 4 = 0; - // - // So far so good, right? Normal points into ball B like it should. - // - // rel_vx = 1 - 0 = 1; - // rel_vy = 0 - -1 = 1; - // rel_dot_n = 1 * 1 + 1 * 0 = 1; - // denom = (1 * 1 + 0 * 0) * (1/10 + 1/100) = 1 * 11/100 = 11/100; - // J = -(1 + 0) * 1 / (11/100) = -100/11; - // dx = 1 * (-100/11) / 10 = -10/11; - // dy = 0; - // p.dx = -1 * (-100/11) / 100 = 1/11; - // p.dy = -0; - // - // So finally: - // vx = 1 - 10/11 = 1/11 - // vy = 0 - // p.vx = 0 + 1/11 = 1/11 - // p.vy = -1 - // - // We see that in with restitution 0, both balls continue at same - // speed along X axis. However at the same time ball B continues its - // path upwards like it should. Seems correct right? -} diff --git a/pybindings/src/noise.cpp b/pybindings/src/noise.cpp deleted file mode 100644 index a88e96d..0000000 --- a/pybindings/src/noise.cpp +++ /dev/null @@ -1,146 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include -#include -#include "noise.hpp" -#include "sqrdmd.hpp" -#include "simplexnoise.hpp" -#include "utils.hpp" - -static const float SQRDMD_ROUGHNESS = 0.35f; -static const float SIMPLEX_PERSISTENCE = 0.25f; -#define PI 3.14159265 - -static uint32_t nearest_pow(uint32_t num) -{ - uint32_t n = 1; - - while (n < num){ - n <<= 1; - } - - return n; -} - -void createSlowNoise(float* map, const WorldDimension& tmpDim, SimpleRandom randsource) -{ - long seed = randsource.next(); - uint32_t width = tmpDim.getWidth(); - uint32_t height = tmpDim.getHeight(); - float persistence = 0.25f; - - float ka = 256/seed; - float kb = seed*567%256; - float kc = (seed*seed) % 256; - float kd = (567-seed) % 256; - for (int y = 0; y < height; y++) { - for (int x = 0; x < width; x++) { - float fNX = x/(float)width; // we let the x-offset define the circle - float fNY = y/(float)height; // we let the x-offset define the circle - float fRdx = fNX*2*PI; // a full circle is two pi radians - float fRdy = fNY*4*PI; // a full circle is two pi radians - float fYSin = sinf(fRdy); - float fRdsSin = 1.0f; - float noiseScale = 0.593; - float a = fRdsSin*sinf(fRdx); - float b = fRdsSin*cosf(fRdx); - float c = fRdsSin*sinf(fRdy); - float d = fRdsSin*cosf(fRdy); - float v = scaled_octave_noise_4d(4.0f, - persistence, - 0.25f, - 0.0f, - 1.0f, - ka+a*noiseScale, - kb+b*noiseScale, - kc+c*noiseScale, - kd+d*noiseScale); - if (map[y * width + x] == 0.0f) map[y * width + x] = v; - } - } -} - -void createNoise(float* tmp, const WorldDimension& tmpDim, SimpleRandom randsource, bool useSimplex) -{ -try { - if (useSimplex) { - simplexnoise(randsource.next(), tmp, - tmpDim.getWidth(), - tmpDim.getHeight(), - SQRDMD_ROUGHNESS); - } else { - uint32_t side = tmpDim.getMax(); - side = nearest_pow(side)+1; - float* squareTmp = new float[side*side]; - memset(squareTmp, 0, sizeof(float)*side*side); - for (int y=0; y // FT_EPSILON -#ifdef __MINGW32__ // this is to avoid a problem with the hypot function which is messed up by Python... -#undef __STRICT_ANSI__ -#endif -#include // sin, cos -#include // rand -#include -#include // std::invalid_argument -#include - -#include "plate.hpp" -#include "heightmap.hpp" -#include "rectangle.hpp" -#include "utils.hpp" -#include "plate_functions.hpp" - -using namespace std; - -plate::plate(long seed, const float* m, uint32_t w, uint32_t h, uint32_t _x, uint32_t _y, - uint32_t plate_age, WorldDimension worldDimension) : - _randsource(seed), - _mass(MassBuilder(m, Dimension(w, h)).build()), - map(w, h), - age_map(w, h), - _worldDimension(worldDimension), - _movement(_randsource, worldDimension) -{ - if (NULL == m) { - throw invalid_argument("the given heightmap should not be null"); - } - if (w <= 0 || h <= 0) { - throw invalid_argument("width and height of the plate should be greater than zero"); - } - if (_x < 0 || _y <0) { - throw invalid_argument("coordinates of the plate should be greater or equal to zero"); - } - if (plate_age < 0) { - throw invalid_argument("age of the plate should be greater or equal to zero"); - } - - const uint32_t plate_area = w * h; - - _bounds = new Bounds(worldDimension, FloatPoint(_x, _y), Dimension(w, h)); - - uint32_t k; - for (uint32_t y = k = 0; y < _bounds->height(); ++y) { - for (uint32_t x = 0; x < _bounds->width(); ++x, ++k) { - // Clone map data and count crust mass. - map[k] = m[k]; - - // Set the age of ALL points in this plate to same - // value. The right thing to do would be to simulate - // the generation of new oceanic crust as if the plate - // had been moving to its current direction until all - // plate's (oceanic) crust receive an age. - age_map.set(x, y, plate_age & -(m[k] > 0)); - } - } - Segments* segments = new Segments(plate_area); - _segments = segments; - _mySegmentCreator = new MySegmentCreator(*_bounds, _segments, map, _worldDimension); - segments->setSegmentCreator(_mySegmentCreator); - segments->setBounds(_bounds); -} - -plate::~plate() -{ - delete _mySegmentCreator; - delete _segments; - delete _bounds; -} - -uint32_t plate::addCollision(uint32_t wx, uint32_t wy) -{ - ISegmentData& seg = getContinentAt(wx, wy); - seg.incCollCount(); - return seg.area(); -} - -void plate::addCrustByCollision(uint32_t x, uint32_t y, float z, uint32_t time, ContinentId activeContinent) -{ - // Add crust. Extend plate if necessary. - setCrust(x, y, getCrust(x, y) + z, time); - - uint32_t index = _bounds->getValidMapIndex(&x, &y); - _segments->setId(index, activeContinent); - - ISegmentData& data = (*_segments)[activeContinent]; - data.incArea(); - data.enlarge_to_contain(x, y); -} - -void plate::addCrustBySubduction(uint32_t x, uint32_t y, float z, uint32_t t, - float dx, float dy) -{ - // TODO: Create an array of coordinate changes that would create - // a circle around current point. Array is static and it is - // initialized at the first call to this function. - // After all points of the circle are checked around subduction - // point the one with most land mass around it will be chosen as - // "most in land" point and subducting crust is added there. - // However to achieve a little more "natural" look normal - // distributed randomness is added around the "center" point. - // Benefits: - // NEVER adds crust outside plate. - // ALWAYS goes inland as much as possible - // Drawbacks: - // Additional logic required - // Might place crust on other continent on same plate! - uint32_t index = _bounds->getValidMapIndex(&x, &y); - - // Take vector difference only between plates that move more or less - // to same direction. This makes subduction direction behave better. - // - // Use of "this" pointer is not necessary, but it make code clearer. - // Cursed be those who use "m_" prefix in member names! >( - float dot = _movement.dot(dx, dy); - dx -= _movement.velocityOnX(dot > 0); - dy -= _movement.velocityOnY(dot > 0); - - float offset = (float)_randsource.next_double(); - float offset_sign = 2 * (int)(_randsource.next() % 2) - 1; - offset *= offset * offset * offset_sign; - float offset2 = (float)_randsource.next_double(); - float offset_sign2 = 2 * (int)(_randsource.next() % 2) - 1; - offset2 *= offset2 * offset2 * offset_sign2; - dx = 10 * dx + 3 * offset; - dy = 10 * dy + 3 * offset2; - - float fx = x + dx; - float fy = y + dy; - - if (_bounds->isInLimits(fx, fy)) - { - index = _bounds->index(fx, fy); - if (map[index] > 0) - { - t = (map[index] * age_map[index] + z * t) / (map[index] + z); - age_map[index] = t * (z > 0); - - map[index] += z; - _mass.incMass(z); - } - } -} - -float plate::aggregateCrust(plate* p, uint32_t wx, uint32_t wy) -{ - uint32_t lx = wx, ly = wy; - const uint32_t index = _bounds->getValidMapIndex(&lx, &ly); - - const ContinentId seg_id = _segments->id(index); - - // This check forces the caller to do things in proper order! - // - // Usually continents collide at several locations simultaneously. - // Thus if this segment that is being merged now is removed from - // segmentation bookkeeping, then the next point of collision that is - // processed during the same iteration step would cause the test - // below to be true and system would experience a premature abort. - // - // Therefore, segmentation bookkeeping is left intact. It doesn't - // cause significant problems because all crust is cleared and empty - // points are not processed at all. (Test on (seg_id >= seg_data.size()) removed) - - // One continent may have many points of collision. If one of them - // causes continent to aggregate then all successive collisions and - // attempts of aggregation would necessarily change nothing at all, - // because the continent was removed from this plate earlier! - if ((*_segments)[seg_id].isEmpty()) { - return 0; // Do not process empty continents. - } - - ContinentId activeContinent = p->selectCollisionSegment(wx, wy); - - // Wrap coordinates around world edges to safeguard subtractions. - wx += _worldDimension.getWidth(); - wy += _worldDimension.getHeight(); - - // Aggregating segment [%u, %u]x[%u, %u] vs. [%u, %u]@[%u, %u]\n", - // seg_data[seg_id].x0, seg_data[seg_id].y0, - // seg_data[seg_id].x1, seg_data[seg_id].y1, - // _dimension.getWidth(), _dimension.getHeight(), lx, ly); - - float old_mass = _mass.getMass(); - - // Add all of the collided continent's crust to destination plate. - for (uint32_t y = (*_segments)[seg_id].getTop(); y <= (*_segments)[seg_id].getBottom(); ++y) - { - for (uint32_t x = (*_segments)[seg_id].getLeft(); x <= (*_segments)[seg_id].getRight(); ++x) - { - const uint32_t i = y * _bounds->width() + x; - if ((_segments->id(i) == seg_id) && (map[i] > 0)) - { - p->addCrustByCollision(wx + x - lx, wy + y - ly, - map[i], age_map[i], activeContinent); - - _mass.incMass(-1.0f * map[i]); - map[i] = 0.0f; - } - } - } - - (*_segments)[seg_id].markNonExistent(); // Mark segment as non-existent - return old_mass - _mass.getMass(); -} - -void plate::applyFriction(float deformed_mass) -{ - // Remove the energy that deformation consumed from plate's kinetic - // energy: F - dF = ma - dF => a = dF/m. - if (!_mass.null()) - { - _movement.applyFriction(deformed_mass, _mass.getMass()); - } -} - -void plate::collide(plate& p, uint32_t wx, uint32_t wy, float coll_mass) -{ - _movement.collide(this->_mass, p, wx, wy, coll_mass); -} - -void plate::calculateCrust(uint32_t x, uint32_t y, uint32_t index, - float& w_crust, float& e_crust, float& n_crust, float& s_crust, - uint32_t& w, uint32_t& e, uint32_t& n, uint32_t& s) -{ - ::calculateCrust(x, y, index, w_crust, e_crust, n_crust, s_crust, - w, e, n, s, - _worldDimension, map, _bounds->width(), _bounds->height()); -} - -void plate::findRiverSources(float lower_bound, vector* sources) -{ - // Find all tops. - for (uint32_t y = 0; y < _bounds->height(); ++y) { - for (uint32_t x = 0; x < _bounds->width(); ++x) { - const uint32_t index = _bounds->index(x, y); - - if (map[index] < lower_bound) { - continue; - } - - float w_crust, e_crust, n_crust, s_crust; - uint32_t w, e, n, s; - calculateCrust(x, y, index, w_crust, e_crust, n_crust, s_crust, - w, e, n, s); - - // This location is either at the edge of the plate or it is not the - // tallest of its neightbours. Don't start a river from here. - if (w_crust * e_crust * n_crust * s_crust == 0) { - continue; - } - - sources->push_back(index); - } - } -} - -void plate::flowRivers(float lower_bound, vector* sources, HeightMap& tmp) -{ - vector sinks_data; - vector* sinks = &sinks_data; - - uint32_t* isDone = new uint32_t[_bounds->area()]; - memset(isDone, 0, _bounds->area() * sizeof(uint32_t)); - - // From each top, start flowing water along the steepest slope. - while (!sources->empty()) { - while (!sources->empty()) { - const uint32_t index = sources->back(); - const uint32_t y = index / _bounds->width(); - const uint32_t x = index - y * _bounds->width(); - - sources->pop_back(); - - if (map[index] < lower_bound) { - continue; - } - - float w_crust, e_crust, n_crust, s_crust; - uint32_t w, e, n, s; - calculateCrust(x, y, index, w_crust, e_crust, n_crust, s_crust, - w, e, n, s); - - // If this is the lowest part of its neighbourhood, stop. - if (w_crust + e_crust + n_crust + s_crust == 0) { - continue; - } - - w_crust += (w_crust == 0) * map[index]; - e_crust += (e_crust == 0) * map[index]; - n_crust += (n_crust == 0) * map[index]; - s_crust += (s_crust == 0) * map[index]; - - // Find lowest neighbour. - float lowest_crust = w_crust; - uint32_t dest = index - 1; - - if (e_crust < lowest_crust) { - lowest_crust = e_crust; - dest = index + 1; - } - - if (n_crust < lowest_crust) { - lowest_crust = n_crust; - dest = index - _bounds->width(); - } - - if (s_crust < lowest_crust) { - lowest_crust = s_crust; - dest = index + _bounds->width(); - } - - // if it's not handled yet, add it as new sink. - if (dest < _bounds->area() && !isDone[dest]) { - sinks->push_back(dest); - isDone[dest] = 1; - } - - // Erode this location with the water flow. - tmp[index] -= (tmp[index] - lower_bound) * 0.2; - } - - - vector* v_tmp = sources; - sources = sinks; - sinks = v_tmp; - sinks->clear(); - } - - delete[] isDone; -} - -void plate::erode(float lower_bound) -{ - vector sources_data; - vector* sources = &sources_data; - - HeightMap tmpHm(map); - findRiverSources(lower_bound, sources); - flowRivers(lower_bound, sources, tmpHm); - - // Add random noise (10 %) to heightmap. - for (uint32_t i = 0; i < _bounds->area(); ++i){ - float alpha = 0.2 * (float)_randsource.next_double(); - tmpHm[i] += 0.1 * tmpHm[i] - alpha * tmpHm[i]; - } - - map = tmpHm; - tmpHm.set_all(0.0f); - MassBuilder massBuilder; - - for (uint32_t y = 0; y < _bounds->height(); ++y) - { - for (uint32_t x = 0; x < _bounds->width(); ++x) - { - const uint32_t index = y * _bounds->width() + x; - massBuilder.addPoint(x, y, map[index]); - tmpHm[index] += map[index]; // Careful not to overwrite earlier amounts. - - if (map[index] < lower_bound) - continue; - - float w_crust, e_crust, n_crust, s_crust; - uint32_t w, e, n, s; - calculateCrust(x, y, index, w_crust, e_crust, n_crust, s_crust, - w, e, n, s); - - // This location has no neighbours (ARTIFACT!) or it is the lowest - // part of its area. In either case the work here is done. - if (w_crust + e_crust + n_crust + s_crust == 0) - continue; - - // The steeper the slope, the more water flows along it. - // The more downhill (sources), the more water flows to here. - // 1+1+10 = 12, avg = 4, stdev = sqrt((3*3+3*3+6*6)/3) = 4.2, var = 18, - // 1*1+1*1+10*10 = 102, 102/4.2=24 - // 1+4+7 = 12, avg = 4, stdev = sqrt((3*3+0*0+3*3)/3) = 2.4, var = 6, - // 1*1+4*4+7*7 = 66, 66/2.4 = 27 - // 4+4+4 = 12, avg = 4, stdev = sqrt((0*0+0*0+0*0)/3) = 0, var = 0, - // 4*4+4*4+4*4 = 48, 48/0 = inf -> 48 - // If there's a source slope of height X then it will always cause - // water erosion of amount Y. Then again from one spot only so much - // water can flow. - // Thus, the calculated non-linear flow value for this location is - // multiplied by the "water erosion" constant. - // The result is max(result, 1.0). New height of this location could - // be e.g. h_lowest + (1 - 1 / result) * (h_0 - h_lowest). - - // Calculate the difference in height between this point and its - // nbours that are lower than this point. - float w_diff = map[index] - w_crust; - float e_diff = map[index] - e_crust; - float n_diff = map[index] - n_crust; - float s_diff = map[index] - s_crust; - - float min_diff = w_diff; - min_diff -= (min_diff - e_diff) * (e_diff < min_diff); - min_diff -= (min_diff - n_diff) * (n_diff < min_diff); - min_diff -= (min_diff - s_diff) * (s_diff < min_diff); - - // Calculate the sum of difference between lower neighbours and - // the TALLEST lower neighbour. - float diff_sum = (w_diff - min_diff) * (w_crust > 0) + - (e_diff - min_diff) * (e_crust > 0) + - (n_diff - min_diff) * (n_crust > 0) + - (s_diff - min_diff) * (s_crust > 0); - - // Erosion difference sum is negative! - assert(diff_sum >= 0); - - if (diff_sum < min_diff) - { - // There's NOT enough room in neighbours to contain all the - // crust from this peak so that it would be as tall as its - // tallest lower neighbour. Thus first step is make ALL - // lower neighbours and this point equally tall. - tmpHm[w] += (w_diff - min_diff) * (w_crust > 0); - tmpHm[e] += (e_diff - min_diff) * (e_crust > 0); - tmpHm[n] += (n_diff - min_diff) * (n_crust > 0); - tmpHm[s] += (s_diff - min_diff) * (s_crust > 0); - tmpHm[index] -= min_diff; - - min_diff -= diff_sum; - - // Spread the remaining crust equally among all lower nbours. - min_diff /= 1 + (w_crust > 0) + (e_crust > 0) + - (n_crust > 0) + (s_crust > 0); - - tmpHm[w] += min_diff * (w_crust > 0); - tmpHm[e] += min_diff * (e_crust > 0); - tmpHm[n] += min_diff * (n_crust > 0); - tmpHm[s] += min_diff * (s_crust > 0); - tmpHm[index] += min_diff; - } - else - { - float unit = min_diff / diff_sum; - - // Remove all crust from this location making it as tall as - // its tallest lower neighbour. - tmpHm[index] -= min_diff; - - // Spread all removed crust among all other lower neighbours. - tmpHm[w] += unit * (w_diff - min_diff) * (w_crust > 0); - tmpHm[e] += unit * (e_diff - min_diff) * (e_crust > 0); - tmpHm[n] += unit * (n_diff - min_diff) * (n_crust > 0); - tmpHm[s] += unit * (s_diff - min_diff) * (s_crust > 0); - } - } - } - map = tmpHm; - _mass = massBuilder.build(); -} - -void plate::getCollisionInfo(uint32_t wx, uint32_t wy, uint32_t* count, float* ratio) const -{ - const ISegmentData& seg = getContinentAt(wx, wy); - - *count = 0; - *ratio = 0; - - *count = seg.collCount(); - *ratio = (float)seg.collCount() / - (float)(1 + seg.area()); // +1 avoids DIV with zero. -} - -uint32_t plate::getContinentArea(uint32_t wx, uint32_t wy) const -{ - const uint32_t index = _bounds->getValidMapIndex(&wx, &wy); - assert(_segments->id(index) < _segments->size()); - return (*_segments)[_segments->id(index)].area(); -} - -float plate::getCrust(uint32_t x, uint32_t y) const -{ - const uint32_t index = _bounds->getMapIndex(&x, &y); - return index != BAD_INDEX ? map[index] : 0; -} - -uint32_t plate::getCrustTimestamp(uint32_t x, uint32_t y) const -{ - const uint32_t index = _bounds->getMapIndex(&x, &y); - return index != BAD_INDEX ? age_map[index] : 0; -} - -void plate::getMap(const float** c, const uint32_t** t) const -{ - if (c) { - *c = map.raw_data(); - } - if (t) { - *t = age_map.raw_data(); - } -} - -void plate::move() -{ - _movement.move(); - - // Location modulations into range [0..world width/height[ are a have to! - // If left undone SOMETHING WILL BREAK DOWN SOMEWHERE in the code! - - _bounds->shift(_movement.velocityOnX(), _movement.velocityOnY()); -} - -void plate::resetSegments() -{ - p_assert(_bounds->area()==_segments->area(), "Segments has not the expected area"); - _segments->reset(); -} - -void plate::setCrust(uint32_t x, uint32_t y, float z, uint32_t t) -{ - if (z < 0) { // Do not accept negative values. - z = 0; - } - - uint32_t _x = x; - uint32_t _y = y; - uint32_t index = _bounds->getMapIndex(&_x, &_y); - - if (index == BAD_INDEX) - { - // Extending plate for nothing! - assert(z>0); - - const uint32_t ilft = _bounds->leftAsUint(); - const uint32_t itop = _bounds->topAsUint(); - const uint32_t irgt = _bounds->rightAsUintNonInclusive(); - const uint32_t ibtm = _bounds->bottomAsUintNonInclusive(); - - _worldDimension.normalize(x, y); - - // Calculate distance of new point from plate edges. - const uint32_t _lft = ilft - x; - const uint32_t _rgt = (_worldDimension.getWidth() & -(x < ilft)) + x - irgt; - const uint32_t _top = itop - y; - const uint32_t _btm = (_worldDimension.getHeight() & -(y < itop)) + y - ibtm; - - // Set larger of horizontal/vertical distance to zero. - // A valid distance is NEVER larger than world's side's length! - uint32_t d_lft = _lft & -(_lft < _rgt) & -(_lft < _worldDimension.getWidth()); - uint32_t d_rgt = _rgt & -(_rgt <= _lft) & -(_rgt < _worldDimension.getWidth()); - uint32_t d_top = _top & -(_top < _btm) & -(_top < _worldDimension.getHeight()); - uint32_t d_btm = _btm & -(_btm <= _top) & -(_btm < _worldDimension.getHeight()); - - // Scale all changes to multiple of 8. - d_lft = ((d_lft > 0) + (d_lft >> 3)) << 3; - d_rgt = ((d_rgt > 0) + (d_rgt >> 3)) << 3; - d_top = ((d_top > 0) + (d_top >> 3)) << 3; - d_btm = ((d_btm > 0) + (d_btm >> 3)) << 3; - - // Make sure plate doesn't grow bigger than the system it's in! - if (_bounds->width() + d_lft + d_rgt > _worldDimension.getWidth()) - { - d_lft = 0; - d_rgt = _worldDimension.getWidth() - _bounds->width(); - } - - if (_bounds->height() + d_top + d_btm > _worldDimension.getHeight()) - { - d_top = 0; - d_btm = _worldDimension.getHeight() - _bounds->height(); - } - - // Index out of bounds, but nowhere to grow! - assert(d_lft + d_rgt + d_top + d_btm != 0); - - const uint32_t old_width = _bounds->width(); - const uint32_t old_height = _bounds->height(); - - _bounds->shift(-1.0*d_lft, -1.0*d_top); - _bounds->grow(d_lft + d_rgt, d_top + d_btm); - - HeightMap tmph = HeightMap(_bounds->width(), _bounds->height()); - AgeMap tmpa = AgeMap(_bounds->width(), _bounds->height()); - uint32_t* tmps = new uint32_t[_bounds->area()]; - tmph.set_all(0); - tmpa.set_all(0); - memset(tmps, 255, _bounds->area()*sizeof(uint32_t)); - - // copy old plate into new. - for (uint32_t j = 0; j < old_height; ++j) - { - const uint32_t dest_i = (d_top + j) * _bounds->width() + d_lft; - const uint32_t src_i = j * old_width; - memcpy(&tmph[dest_i], &map[src_i], old_width * - sizeof(float)); - memcpy(&tmpa[dest_i], &age_map[src_i], old_width * - sizeof(uint32_t)); - memcpy(&tmps[dest_i], &_segments->id(src_i), old_width * - sizeof(uint32_t)); - } - - map = tmph; - age_map = tmpa; - _segments->reassign(_bounds->area(), tmps); - - // Shift all segment data to match new coordinates. - _segments->shift(d_lft, d_top); - - _x = x, _y = y; - index = _bounds->getValidMapIndex(&_x, &_y); - - assert(index < _bounds->area()); - } - - // Update crust's age. - // If old crust exists, new age is mean of original and supplied ages. - // If no new crust is added, original time remains intact. - const uint32_t old_crust = -(map[index] > 0); - const uint32_t new_crust = -(z > 0); - t = (t & ~old_crust) | ((uint32_t)((map[index] * age_map[index] + z * t) / - (map[index] + z)) & old_crust); - age_map[index] = (t & new_crust) | (age_map[index] & ~new_crust); - - _mass.incMass(-1.0f * map[index]); - _mass.incMass(z); // Update mass counter. - map[index] = z; // Set new crust height to desired location. -} - -ContinentId plate::selectCollisionSegment(uint32_t coll_x, uint32_t coll_y) -{ - uint32_t index = _bounds->getValidMapIndex(&coll_x, &coll_y); - ContinentId activeContinent = _segments->id(index); - return activeContinent; -} - -/////////////////////////////////////////////////////////////////////////////// -/// Private methods /////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// - -uint32_t plate::createSegment(uint32_t x, uint32_t y) throw() -{ - return _mySegmentCreator->createSegment(x, y); -} - -ISegmentData& plate::getContinentAt(int x, int y) -{ - return (*_segments)[_segments->getContinentAt(x, y)]; -} - -const ISegmentData& plate::getContinentAt(int x, int y) const -{ - return (*_segments)[_segments->getContinentAt(x, y)]; -} diff --git a/pybindings/src/plate_functions.cpp b/pybindings/src/plate_functions.cpp deleted file mode 100644 index 85a580a..0000000 --- a/pybindings/src/plate_functions.cpp +++ /dev/null @@ -1,67 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "plate_functions.hpp" - -void calculateCrust( - uint32_t x, uint32_t y, - uint32_t index, - float& w_crust, float& e_crust, float& n_crust, float& s_crust, - uint32_t& w, uint32_t& e, uint32_t& n, uint32_t& s, - const WorldDimension& worldDimension, - HeightMap& map, - const uint32_t width, const uint32_t height) -{ -try { - // Build masks for accessible directions (4-way). - // Allow wrapping around map edges if plate has world wide dimensions. - uint32_t w_mask = -((x > 0) | (width == worldDimension.getWidth())); - uint32_t e_mask = -((x < width - 1) | (width == worldDimension.getWidth())); - uint32_t n_mask = -((y > 0) | (height == worldDimension.getHeight())); - uint32_t s_mask = -((y < height - 1) | (height == worldDimension.getHeight())); - - // Calculate the x and y offset of neighbour directions. - // If neighbour is out of plate edges, set it to zero. This protects - // map memory reads from segment faulting. - w = w_mask==-1 ? worldDimension.xMod(x-1) : 0; - e = e_mask==-1 ? worldDimension.xMod(x+1) : 0; - n = n_mask==-1 ? worldDimension.yMod(y-1) : 0; - s = s_mask==-1 ? worldDimension.yMod(y+1) : 0; - - // Calculate offsets within map memory. - w = y * width + w; - e = y * width + e; - n = n * width + x; - s = s * width + x; - - // Extract neighbours heights. Apply validity filtering: 0 is invalid. - w_crust = map[w] * (w_mask & (map[w] < map[index])); - e_crust = map[e] * (e_mask & (map[e] < map[index])); - n_crust = map[n] * (n_mask & (map[n] < map[index])); - s_crust = map[s] * (s_mask & (map[s] < map[index])); -} catch (const exception& e){ - std::string msg = "Problem during plate::calculateCrust (width: "; - msg = msg + Platec::to_string(width) - + ", height: " + Platec::to_string(height) - + ", x: " + Platec::to_string(x) - + ", y:" + Platec::to_string(y) + ") :" - + e.what(); - throw runtime_error(msg.c_str()); -} -} diff --git a/pybindings/src/platecapi.cpp b/pybindings/src/platecapi.cpp deleted file mode 100644 index 0a98a76..0000000 --- a/pybindings/src/platecapi.cpp +++ /dev/null @@ -1,140 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "lithosphere.hpp" -#include "plate.hpp" -#include "platecapi.hpp" -#include -#include - -#include - -class platec_api_list_elem -{ - public: - platec_api_list_elem(uint32_t _id, lithosphere* _data) : - data(_data), id(_id) { } - - lithosphere* data; - uint32_t id; -}; - -extern lithosphere* platec_api_get_lithosphere(uint32_t); - -static std::vector lithospheres; -static uint32_t last_id = 1; - - -void* platec_api_create(long seed, uint32_t width, uint32_t height, float sea_level, - uint32_t erosion_period, float folding_ratio, - uint32_t aggr_overlap_abs, float aggr_overlap_rel, - uint32_t cycle_count, uint32_t num_plates) -{ - /* Miten nykyisen opengl-mainin koodit refaktoroidaan tänne? - * parametrien tarkistus, kommentit eli dokumentointi, muuta? */ - - lithosphere* litho = new lithosphere(seed, width, height, sea_level, - erosion_period, folding_ratio, aggr_overlap_abs, - aggr_overlap_rel, cycle_count); - litho->createPlates(num_plates); - - platec_api_list_elem elem(++last_id, litho); - lithospheres.push_back(elem); - - return litho; -} - -void platec_api_destroy(void* litho) -{ - for (uint32_t i = 0; i < lithospheres.size(); ++i) - if (lithospheres[i].data == litho) { - lithospheres.erase(lithospheres.begin()+i); - break; - } -} - -const uint32_t* platec_api_get_agemap(uint32_t id) -{ - lithosphere* litho = platec_api_get_lithosphere(id); - if (!litho) - return NULL; - - return litho->getAgemap(); -} - -float* platec_api_get_heightmap(void *pointer) -{ - lithosphere* litho = (lithosphere*)pointer; - float *res = litho->getTopography(); - return res; -} - -uint32_t* platec_api_get_platesmap(void *pointer) -{ - lithosphere* litho = (lithosphere*)pointer; - uint32_t *res = litho->getPlatesMap(); - return res; -} - -lithosphere* platec_api_get_lithosphere(uint32_t id) -{ - for (uint32_t i = 0; i < lithospheres.size(); ++i) - if (lithospheres[i].id == id) - return lithospheres[i].data; - - return NULL; -} - -uint32_t platec_api_is_finished(void *pointer) -{ - lithosphere* litho = (lithosphere*)pointer; - if (litho->isFinished()) { - return 1; - } else { - return 0; - } -} - -void platec_api_step(void *pointer) -{ - lithosphere* litho = (lithosphere*)pointer; - litho->update(); -} - -uint32_t lithosphere_getMapWidth ( void* object) -{ - return static_cast( object)->getWidth(); -} - -uint32_t lithosphere_getMapHeight ( void* object) -{ - return static_cast( object)->getHeight(); -} - -float platec_api_velocity_unity_vector_x(void* pointer, uint32_t plate_index) -{ - lithosphere* litho = (lithosphere*)pointer; - return litho->getPlate(plate_index)->velocityUnitVector().x(); -} - -float platec_api_velocity_unity_vector_y(void* pointer, uint32_t plate_index) -{ - lithosphere* litho = (lithosphere*)pointer; - return litho->getPlate(plate_index)->velocityUnitVector().y(); -} diff --git a/pybindings/src/rectangle.cpp b/pybindings/src/rectangle.cpp deleted file mode 100644 index 4bcfca5..0000000 --- a/pybindings/src/rectangle.cpp +++ /dev/null @@ -1,92 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include -#include "rectangle.hpp" -#include "utils.hpp" - -namespace Platec { - -/// Return a valid index or BAD_INDEX -uint32_t Rectangle::getMapIndex(uint32_t* px, uint32_t* py) const -{ - uint32_t x = *px % _worldDimension.getWidth(); - uint32_t y = *py % _worldDimension.getHeight(); - - const uint32_t ilft = (uint32_t)(int)_left; - const uint32_t itop = (uint32_t)(int)_top; - const uint32_t irgt = (uint32_t)(int)_right + (((uint32_t)(int)_right < ilft) ? (uint32_t)(int)_worldDimension.getWidth() : 0); - const uint32_t ibtm = (uint32_t)(int)_bottom + (((uint32_t)(int)_bottom < itop) ? (uint32_t)(int)_worldDimension.getHeight() : 0); - const int width = irgt - ilft; - if (width < 0) { - throw std::invalid_argument("(Rectangle::getMapIndex) negative width"); - } - - /////////////////////////////////////////////////////////////////////// - // If you think you're smart enough to optimize this then PREPARE to be - // smart as HELL to debug it! - /////////////////////////////////////////////////////////////////////// - - const uint32_t xOkA = (x >= ilft) && (x < irgt); - const uint32_t xOkB = (x + _worldDimension.getWidth() >= ilft) && (x + _worldDimension.getWidth() < irgt); - const uint32_t xOk = xOkA || xOkB; - - const uint32_t yOkA = (y >= itop) && (y < ibtm); - const uint32_t yOkB = (y + _worldDimension.getHeight() >= itop) && (y + _worldDimension.getHeight() < ibtm); - const uint32_t yOk = yOkA || yOkB; - - x += (x < ilft) ? _worldDimension.getWidth() : 0; // Point is within plate's map: wrap - y += (y < itop) ? _worldDimension.getHeight() : 0; // it around world edges if necessary. - - x -= ilft; // Calculate offset within local map. - y -= itop; - - if (x < 0) { - throw std::invalid_argument("failure x"); - } - if (y < 0) { - throw std::invalid_argument("failure y"); - } - - if (xOk && yOk) { - *px = x; - *py = y; - return (y * width + x); - } else { - return BAD_INDEX; - } -} - -void Rectangle::enlarge_to_contain(uint32_t x, uint32_t y) -{ - if (y < _top) { - _top = y; - } - if (y > _bottom) { - _bottom = y; - } - if (x < _left) { - _left = x; - } - if (x > _right) { - _right = x; - } -} - -} \ No newline at end of file diff --git a/pybindings/src/segment_creator.cpp b/pybindings/src/segment_creator.cpp deleted file mode 100644 index cf6f174..0000000 --- a/pybindings/src/segment_creator.cpp +++ /dev/null @@ -1,259 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "segment_creator.hpp" -#include "movement.hpp" -#include "segments.hpp" -#include "bounds.hpp" - -uint32_t MySegmentCreator::calcDirection(uint32_t x, uint32_t y, const uint32_t origin_index, const uint32_t ID) const -{ - uint32_t canGoLeft = x > 0 && map[origin_index - 1] >= CONT_BASE; - uint32_t canGoRight = x < _bounds.width() - 1 && map[origin_index+1] >= CONT_BASE; - uint32_t canGoUp = y > 0 && map[origin_index - _bounds.width()] >= CONT_BASE; - uint32_t canGoDown = y < _bounds.height() - 1 && map[origin_index + _bounds.width()] >= CONT_BASE; - uint32_t nbour_id = ID; - - // This point belongs to no segment yet. - // However it might be a neighbour to some segment created earlier. - // If such neighbour is found, associate this point with it. - if (canGoLeft && _segments->id(origin_index - 1) < ID) { - nbour_id = _segments->id(origin_index - 1); - } else if (canGoRight && _segments->id(origin_index + 1) < ID) { - nbour_id = _segments->id(origin_index + 1); - } else if (canGoUp && _segments->id(origin_index - _bounds.width()) < ID) { - nbour_id = _segments->id(origin_index - _bounds.width()); - } else if (canGoDown && _segments->id(origin_index + _bounds.width()) < ID) { - nbour_id = _segments->id(origin_index + _bounds.width()); - } - - return nbour_id; -} - -void MySegmentCreator::scanSpans(const uint32_t line, uint32_t& start, uint32_t& end, - std::vector* spans_todo, std::vector* spans_done) const -{ - do // Find an unscanned span on this line. - { - end = spans_todo[line].back(); - spans_todo[line].pop_back(); - - start = spans_todo[line].back(); - spans_todo[line].pop_back(); - - // Reduce any done spans from this span. - for (uint32_t j = 0; j < spans_done[line].size(); - j += 2) - { - // Saved coordinates are AT the point - // that was included last to the span. - // That's why equalities matter. - - if (start >= spans_done[line][j] && - start <= spans_done[line][j+1]) - start = spans_done[line][j+1] + 1; - - if (end >= spans_done[line][j] && - end <= spans_done[line][j+1]) - end = spans_done[line][j] - 1; - } - - // Unsigned-ness hacking! - // Required to fix the underflow of end - 1. - start |= -(end >= _bounds.width()); - end -= (end >= _bounds.width()); - - } while (start > end && spans_todo[line].size()); -} - -ContinentId MySegmentCreator::createSegment(uint32_t x, uint32_t y) const throw() -{ - const uint32_t origin_index = _bounds.index(x, y); - const uint32_t ID = _segments->size(); - - if (_segments->id(origin_index) < ID) { - return _segments->id(origin_index); - } - - uint32_t nbour_id = calcDirection(x, y, origin_index, ID); - - if (nbour_id < ID) - { - _segments->setId(origin_index, nbour_id); - (*_segments)[nbour_id].incArea(); - - (*_segments)[nbour_id].enlarge_to_contain(x, y); - - return nbour_id; - } - - uint32_t lines_processed; - Platec::Rectangle r = Platec::Rectangle(_worldDimension, x, x, y, y); - - SegmentData* pData = new SegmentData(r, 0); - - std::vector* spans_todo = new std::vector[_bounds.height()]; - std::vector* spans_done = new std::vector[_bounds.height()]; - - _segments->setId(origin_index, ID); - spans_todo[y].push_back(x); - spans_todo[y].push_back(x); - - do - { - lines_processed = 0; - for (uint32_t line = 0; line < _bounds.height(); ++line) - { - uint32_t start, end; - - if (spans_todo[line].size() == 0) - continue; - - scanSpans(line, start, end, spans_todo, spans_done); - - if (start > end) // Nothing to do here anymore... - continue; - - // Calculate line indices. Allow wrapping around map edges. - const uint32_t row_above = ((line - 1) & -(line > 0)) | - ((_bounds.height() - 1) & -(line == 0)); - const uint32_t row_below = (line + 1) & -(line < _bounds.height() - 1); - const uint32_t line_here = line * _bounds.width(); - const uint32_t line_above = row_above * _bounds.width(); - const uint32_t line_below = row_below * _bounds.width(); - - // Extend the beginning of line. - while (start > 0 && _segments->id(line_here+start-1) > ID && - map[line_here+start-1] >= CONT_BASE) - { - --start; - _segments->setId(line_here + start, ID); - - // Count volume of pixel... - } - - // Extend the end of line. - while (end < _bounds.width() - 1 && - _segments->id(line_here + end + 1) > ID && - map[line_here + end + 1] >= CONT_BASE) - { - ++end; - _segments->setId(line_here + end, ID); - - // Count volume of pixel... - } - - // Check if should wrap around left edge. - if (_bounds.width() == _worldDimension.getWidth() && start == 0 && - _segments->id(line_here+_bounds.width()-1) > ID && - map[line_here+_bounds.width()-1] >= CONT_BASE) - { - _segments->setId(line_here + _bounds.width() - 1, ID); - spans_todo[line].push_back(_bounds.width() - 1); - spans_todo[line].push_back(_bounds.width() - 1); - - // Count volume of pixel... - } - - // Check if should wrap around right edge. - if (_bounds.width() == _worldDimension.getWidth() && end == _bounds.width() - 1 && - _segments->id(line_here+0) > ID && - map[line_here+0] >= CONT_BASE) - { - _segments->setId(line_here + 0, ID); - spans_todo[line].push_back(0); - spans_todo[line].push_back(0); - - // Count volume of pixel... - } - - pData->incArea(1 + end - start); // Update segment area counter. - - // Record any changes in extreme dimensions. - if (line < pData->getTop()) pData->setTop(line); - if (line > pData->getBottom()) pData->setBottom(line); - if (start < pData->getLeft()) pData->setLeft(start); - if (end > pData->getRight()) pData->setRight(end); - - if (line > 0 || _bounds.height() == _worldDimension.getHeight()) { - for (uint32_t j = start; j <= end; ++j) - if (_segments->id(line_above + j) > ID && - map[line_above + j] >= CONT_BASE) - { - uint32_t a = j; - _segments->setId(line_above + a, ID); - - // Count volume of pixel... - - while (++j < _bounds.width() && - _segments->id(line_above + j) > ID && - map[line_above + j] >= CONT_BASE) - { - _segments->setId(line_above + j, ID); - - // Count volume of pixel... - } - - uint32_t b = --j; // Last point is invalid. - - spans_todo[row_above].push_back(a); - spans_todo[row_above].push_back(b); - ++j; // Skip the last scanned point. - } - } - - if (line < _bounds.height() - 1 || _bounds.height() == _worldDimension.getHeight()) { - for (uint32_t j = start; j <= end; ++j) - if (_segments->id(line_below + j) > ID && - map[line_below + j] >= CONT_BASE) - { - uint32_t a = j; - _segments->setId(line_below + a, ID); - - // Count volume of pixel... - - while (++j < _bounds.width() && - _segments->id(line_below + j) > ID && - map[line_below + j] >= CONT_BASE) - { - _segments->setId(line_below + j, ID); - - // Count volume of pixel... - } - - uint32_t b = --j; // Last point is invalid. - - spans_todo[row_below].push_back(a); - spans_todo[row_below].push_back(b); - ++j; // Skip the last scanned point. - } - } - - spans_done[line].push_back(start); - spans_done[line].push_back(end); - ++lines_processed; - } - } while (lines_processed > 0); - - delete[] spans_todo; - delete[] spans_done; - _segments->add(pData); - - return ID; -} diff --git a/pybindings/src/segment_data.cpp b/pybindings/src/segment_data.cpp deleted file mode 100644 index 2737dc0..0000000 --- a/pybindings/src/segment_data.cpp +++ /dev/null @@ -1,109 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "segment_data.hpp" - -SegmentData::SegmentData(Platec::Rectangle& rectangle, - uint32_t area) : _rectangle(rectangle), - _area(area), _coll_count(0) {}; - -void SegmentData::enlarge_to_contain(uint32_t x, uint32_t y) -{ - _rectangle.enlarge_to_contain(x, y); -}; - -uint32_t SegmentData::getLeft() const -{ - return _rectangle.getLeft(); -}; - -uint32_t SegmentData::getRight() const -{ - return _rectangle.getRight(); -}; - -uint32_t SegmentData::getTop() const -{ - return _rectangle.getTop(); -}; - -uint32_t SegmentData::getBottom() const -{ - return _rectangle.getBottom(); -}; - -void SegmentData::shift(uint32_t dx, uint32_t dy) -{ - _rectangle.shift(dx, dy); -}; - -void SegmentData::setLeft(uint32_t v) -{ - _rectangle.setLeft(v); -}; - -void SegmentData::setRight(uint32_t v) -{ - _rectangle.setRight(v); -}; - -void SegmentData::setTop(uint32_t v) -{ - _rectangle.setTop(v); -}; - -void SegmentData::setBottom(uint32_t v) -{ - _rectangle.setBottom(v); -}; - -bool SegmentData::isEmpty() const -{ - return _area == 0; -}; - -void SegmentData::incCollCount() -{ - _coll_count++; -}; - -void SegmentData::incArea() -{ - _area++; -}; - -void SegmentData::incArea(uint32_t amount) -{ - _area += amount; -}; - -uint32_t SegmentData::area() const -{ - return _area; -}; - -uint32_t SegmentData::collCount() const -{ - return _coll_count; -} - -void SegmentData::markNonExistent() -{ - _area = 0; -} diff --git a/pybindings/src/segments.cpp b/pybindings/src/segments.cpp deleted file mode 100644 index 0545fcf..0000000 --- a/pybindings/src/segments.cpp +++ /dev/null @@ -1,128 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "segments.hpp" - -Segments::Segments(uint32_t plate_area) -{ - _area = plate_area; - segment = new uint32_t[plate_area]; - memset(segment, 255, plate_area * sizeof(uint32_t)); -} - -Segments::~Segments() -{ - delete[] segment; - segment = NULL; - _area = 0; - for (int i=0; ishift(d_lft, d_top); - } -} - -uint32_t Segments::size() const -{ - return seg_data.size(); -} - -const ISegmentData& Segments::operator[](uint32_t index) const -{ - if (index >= seg_data.size()) throw runtime_error("(Segments::operator[]) unvalid index"); - return *seg_data[index]; -} - -ISegmentData& Segments::operator[](uint32_t index) -{ - if (index >= seg_data.size()) throw runtime_error("(Segments::operator[]) unvalid index"); - return *seg_data[index]; -} - -void Segments::add(ISegmentData* data){ - seg_data.push_back(data); -} - -const ContinentId& Segments::id(uint32_t index) const { - if (index>=_area) { - throw runtime_error("unvalid index"); - } - return segment[index]; -} - -ContinentId& Segments::id(uint32_t index) { - if (index>=_area) { - throw runtime_error("unvalid index"); - } - return segment[index]; -} - -void Segments::setId(uint32_t index, ContinentId id) { - if (index>=_area) { - throw runtime_error("unvalid index"); - } - segment[index] = id; -} - -ContinentId Segments::getContinentAt(int x, int y) const -{ - if (_bounds == NULL) throw runtime_error("(Segments::getContinentAt) bounds not set"); - if (_segmentCreator == NULL) throw runtime_error("(Segments::getContinentAt) segmentCreator not set"); - uint32_t lx = x, ly = y; - uint32_t index = _bounds->getValidMapIndex(&lx, &ly); - ContinentId seg = id(index); - - if (seg >= size()) { - // in this case, we consider as const this call because we calculate - // something that we would calculate anyway, so the segments are - // a sort of cache - //seg = const_cast(this)->createSegment(lx, ly); - seg = _segmentCreator->createSegment(lx, ly); - } - - if (seg >= size()) - { - throw invalid_argument("Could not create segment"); - } - return seg; -} diff --git a/pybindings/src/simplerandom.cpp b/pybindings/src/simplerandom.cpp deleted file mode 100644 index f66f6e0..0000000 --- a/pybindings/src/simplerandom.cpp +++ /dev/null @@ -1,145 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2010 Craig McQueen (http://craig.mcqueen.id.au) - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This is code from the Simple Pseudo-random Number Generators - * Available on GitHub https://github.com/cmcqueen/simplerandom - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "simplerandom.hpp" -#include -#include "utils.hpp" - -void simplerandom_cong_seed(SimpleRandomCong_t * p_cong, uint32_t seed); -void simplerandom_cong_mix(SimpleRandomCong_t * p_cong, const uint32_t * p_data, uint32_t num_data); -uint32_t simplerandom_cong_next(SimpleRandomCong_t * p_cong); - -SimpleRandom::SimpleRandom(uint32_t seed) -{ - internal = new SimpleRandomCong_t(); - simplerandom_cong_seed(internal, seed); -} - -SimpleRandom::SimpleRandom(SimpleRandom& other) -{ - internal = new SimpleRandomCong_t(); - internal->cong = other.internal->cong; -} - -SimpleRandom::SimpleRandom(const SimpleRandom& other) -{ - internal = new SimpleRandomCong_t(); - internal->cong = other.internal->cong; -} - -SimpleRandom::~SimpleRandom() -{ - delete internal; -} - -uint32_t SimpleRandom::next() -{ - uint32_t res = simplerandom_cong_next(internal); - - return res; -} - -double SimpleRandom::next_double() -{ - return ((double)next() / (double)maximum()); -} - -float SimpleRandom::next_float_signed() -{ - float value = next_double(); - p_assert(value >= 0.0f && value <= 1.0f, "(SimpleRandom::next_float_signed)"); - return next_double() - 0.5f; -} - -int32_t SimpleRandom::next_signed() -{ - int32_t value = (int32_t)next(); - return value; -} - -uint32_t SimpleRandom::maximum() -{ - return 4294967295; -} - -uint32_t simplerandom_cong_num_seeds(const SimpleRandomCong_t * p_cong) -{ - (const void *)p_cong; /* We only use this parameter for type checking. */ - - return 1u; -} - -uint32_t simplerandom_cong_seed_array(SimpleRandomCong_t * p_cong, const uint32_t * p_seeds, uint32_t num_seeds, bool mix_extras) -{ - uint32_t seed = 0; - uint32_t num_seeds_used = 0; - - if (num_seeds >= 1u && p_seeds != NULL) - { - seed = p_seeds[0]; - num_seeds_used = 1u; - } - simplerandom_cong_seed(p_cong, seed); - - if (mix_extras && p_seeds != NULL) - { - simplerandom_cong_mix(p_cong, p_seeds + num_seeds_used, num_seeds - num_seeds_used); - num_seeds_used = num_seeds; - } - return num_seeds_used; -} - -void simplerandom_cong_seed(SimpleRandomCong_t * p_cong, uint32_t seed) -{ - p_cong->cong = seed; - /* No sanitize is needed because for Cong, all state values are valid. */ -} - -void simplerandom_cong_sanitize(SimpleRandomCong_t * p_cong) -{ - /* All state values are valid for Cong. No sanitizing needed. */ - (const void *) p_cong; -} - -uint32_t simplerandom_cong_next(SimpleRandomCong_t * p_cong) -{ - uint32_t cong; - - cong = UINT32_C(69069) * p_cong->cong + 12345u; - p_cong->cong = cong; - - return cong; -} - -void simplerandom_cong_mix(SimpleRandomCong_t * p_cong, const uint32_t * p_data, uint32_t num_data) -{ - if (p_data != NULL) - { - while (num_data) - { - --num_data; - p_cong->cong ^= *p_data++; - simplerandom_cong_sanitize(p_cong); - simplerandom_cong_next(p_cong); - } - } -} diff --git a/pybindings/src/simplexnoise.cpp b/pybindings/src/simplexnoise.cpp deleted file mode 100644 index 98412f4..0000000 --- a/pybindings/src/simplexnoise.cpp +++ /dev/null @@ -1,521 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2007-2012 Eliot Eshelman - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This is code from the Simple Pseudo-random Number Generators - * Available on GitHub https://github.com/cmcqueen/simplerandom - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#ifdef __MINGW32__ // this is to avoid a problem with the hypot function which is messed up by Python... -#undef __STRICT_ANSI__ -#endif -#include -#include - -#include "simplexnoise.hpp" - - -/* 2D, 3D and 4D Simplex Noise functions return 'random' values in (-1, 1). - -This algorithm was originally designed by Ken Perlin, but my code has been -adapted from the implementation written by Stefan Gustavson (stegu@itn.liu.se) - -Raw Simplex noise functions return the value generated by Ken's algorithm. - -Scaled Raw Simplex noise functions adjust the range of values returned from the -traditional (-1, 1) to whichever bounds are passed to the function. - -Multi-Octave Simplex noise functions compine multiple noise values to create a -more complex result. Each successive layer of noise is adjusted and scaled. - -Scaled Multi-Octave Simplex noise functions scale the values returned from the -traditional (-1,1) range to whichever range is passed to the function. - -In many cases, you may think you only need a 1D noise function, but in practice -2D is almost always better. For instance, if you're using the current frame -number as the parameter for the noise, all objects will end up with the same -noise value at each frame. By adding a second parameter on the second -dimension, you can ensure that each gets a unique noise value and they don't -all look identical. -*/ - - -// 2D Multi-octave Simplex noise. -// -// For each octave, a higher frequency/lower amplitude function will be added to the original. -// The higher the persistence [0-1], the more of each succeeding octave will be added. -float octave_noise_2d( const float octaves, const float persistence, const float scale, const float x, const float y ) { - float total = 0; - float frequency = scale; - float amplitude = 1; - - // We have to keep track of the largest possible amplitude, - // because each octave adds more, and we need a value in [-1, 1]. - float maxAmplitude = 0; - - for( int i=0; i < octaves; i++ ) { - total += raw_noise_2d( x * frequency, y * frequency ) * amplitude; - - frequency *= 2; - maxAmplitude += amplitude; - amplitude *= persistence; - } - - return total / maxAmplitude; -} - - -// 3D Multi-octave Simplex noise. -// -// For each octave, a higher frequency/lower amplitude function will be added to the original. -// The higher the persistence [0-1], the more of each succeeding octave will be added. -float octave_noise_3d( const float octaves, const float persistence, const float scale, const float x, const float y, const float z ) { - float total = 0; - float frequency = scale; - float amplitude = 1; - - // We have to keep track of the largest possible amplitude, - // because each octave adds more, and we need a value in [-1, 1]. - float maxAmplitude = 0; - - for( int i=0; i < octaves; i++ ) { - total += raw_noise_3d( x * frequency, y * frequency, z * frequency ) * amplitude; - - frequency *= 2; - maxAmplitude += amplitude; - amplitude *= persistence; - } - - return total / maxAmplitude; -} - - -// 4D Multi-octave Simplex noise. -// -// For each octave, a higher frequency/lower amplitude function will be added to the original. -// The higher the persistence [0-1], the more of each succeeding octave will be added. -float octave_noise_4d( const float octaves, const float persistence, const float scale, const float x, const float y, const float z, const float w ) { - float total = 0; - float frequency = scale; - float amplitude = 1; - - // We have to keep track of the largest possible amplitude, - // because each octave adds more, and we need a value in [-1, 1]. - float maxAmplitude = 0; - - for( int i=0; i < octaves; i++ ) { - total += raw_noise_4d( x * frequency, y * frequency, z * frequency, w * frequency ) * amplitude; - - frequency *= 2; - maxAmplitude += amplitude; - amplitude *= persistence; - } - - return total / maxAmplitude; -} - - - -// 2D Scaled Multi-octave Simplex noise. -// -// Returned value will be between loBound and hiBound. -float scaled_octave_noise_2d( const float octaves, const float persistence, const float scale, const float loBound, const float hiBound, const float x, const float y ) { - return octave_noise_2d(octaves, persistence, scale, x, y) * (hiBound - loBound) / 2 + (hiBound + loBound) / 2; -} - - -// 3D Scaled Multi-octave Simplex noise. -// -// Returned value will be between loBound and hiBound. -float scaled_octave_noise_3d( const float octaves, const float persistence, const float scale, const float loBound, const float hiBound, const float x, const float y, const float z ) { - return octave_noise_3d(octaves, persistence, scale, x, y, z) * (hiBound - loBound) / 2 + (hiBound + loBound) / 2; -} - -// 4D Scaled Multi-octave Simplex noise. -// -// Returned value will be between loBound and hiBound. -float scaled_octave_noise_4d( const float octaves, const float persistence, const float scale, const float loBound, const float hiBound, const float x, const float y, const float z, const float w ) { - return octave_noise_4d(octaves, persistence, scale, x, y, z, w) * (hiBound - loBound) / 2 + (hiBound + loBound) / 2; -} - - - -// 2D Scaled Simplex raw noise. -// -// Returned value will be between loBound and hiBound. -float scaled_raw_noise_2d( const float loBound, const float hiBound, const float x, const float y ) { - return raw_noise_2d(x, y) * (hiBound - loBound) / 2 + (hiBound + loBound) / 2; -} - - -// 3D Scaled Simplex raw noise. -// -// Returned value will be between loBound and hiBound. -float scaled_raw_noise_3d( const float loBound, const float hiBound, const float x, const float y, const float z ) { - return raw_noise_3d(x, y, z) * (hiBound - loBound) / 2 + (hiBound + loBound) / 2; -} - -// 4D Scaled Simplex raw noise. -// -// Returned value will be between loBound and hiBound. -float scaled_raw_noise_4d( const float loBound, const float hiBound, const float x, const float y, const float z, const float w ) { - return raw_noise_4d(x, y, z, w) * (hiBound - loBound) / 2 + (hiBound + loBound) / 2; -} - - - -// 2D raw Simplex noise -float raw_noise_2d( const float x, const float y ) { - // Noise contributions from the three corners - float n0, n1, n2; - - // Skew the input space to determine which simplex cell we're in - float F2 = 0.5 * (sqrtf(3.0) - 1.0); - // Hairy factor for 2D - float s = (x + y) * F2; - int i = fastfloor( x + s ); - int j = fastfloor( y + s ); - - float G2 = (3.0 - sqrtf(3.0)) / 6.0; - float t = (i + j) * G2; - // Unskew the cell origin back to (x,y) space - float X0 = i-t; - float Y0 = j-t; - // The x,y distances from the cell origin - float x0 = x-X0; - float y0 = y-Y0; - - // For the 2D case, the simplex shape is an equilateral triangle. - // Determine which simplex we are in. - int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords - if(x0>y0) {i1=1; j1=0;} // lower triangle, XY order: (0,0)->(1,0)->(1,1) - else {i1=0; j1=1;} // upper triangle, YX order: (0,0)->(0,1)->(1,1) - - // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and - // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where - // c = (3-sqrt(3))/6 - float x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords - float y1 = y0 - j1 + G2; - float x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords - float y2 = y0 - 1.0 + 2.0 * G2; - - // Work out the hashed gradient indices of the three simplex corners - int ii = i & 255; - int jj = j & 255; - int gi0 = perm[ii+perm[jj]] % 12; - int gi1 = perm[ii+i1+perm[jj+j1]] % 12; - int gi2 = perm[ii+1+perm[jj+1]] % 12; - - // Calculate the contribution from the three corners - float t0 = 0.5 - x0*x0-y0*y0; - if(t0<0) n0 = 0.0; - else { - t0 *= t0; - n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient - } - - float t1 = 0.5 - x1*x1-y1*y1; - if(t1<0) n1 = 0.0; - else { - t1 *= t1; - n1 = t1 * t1 * dot(grad3[gi1], x1, y1); - } - - float t2 = 0.5 - x2*x2-y2*y2; - if(t2<0) n2 = 0.0; - else { - t2 *= t2; - n2 = t2 * t2 * dot(grad3[gi2], x2, y2); - } - - // Add contributions from each corner to get the final noise value. - // The result is scaled to return values in the interval [-1,1]. - return 70.0 * (n0 + n1 + n2); -} - - -// 3D raw Simplex noise -float raw_noise_3d( const float x, const float y, const float z ) { - float n0, n1, n2, n3; // Noise contributions from the four corners - - // Skew the input space to determine which simplex cell we're in - float F3 = 1.0/3.0; - float s = (x+y+z)*F3; // Very nice and simple skew factor for 3D - int i = fastfloor(x+s); - int j = fastfloor(y+s); - int k = fastfloor(z+s); - - float G3 = 1.0/6.0; // Very nice and simple unskew factor, too - float t = (i+j+k)*G3; - float X0 = i-t; // Unskew the cell origin back to (x,y,z) space - float Y0 = j-t; - float Z0 = k-t; - float x0 = x-X0; // The x,y,z distances from the cell origin - float y0 = y-Y0; - float z0 = z-Z0; - - // For the 3D case, the simplex shape is a slightly irregular tetrahedron. - // Determine which simplex we are in. - int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords - int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords - - if(x0>=y0) { - if(y0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } // X Y Z order - else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } // X Z Y order - else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } // Z X Y order - } - else { // x0 y0) ? 32 : 0; - int c2 = (x0 > z0) ? 16 : 0; - int c3 = (y0 > z0) ? 8 : 0; - int c4 = (x0 > w0) ? 4 : 0; - int c5 = (y0 > w0) ? 2 : 0; - int c6 = (z0 > w0) ? 1 : 0; - int c = c1 + c2 + c3 + c4 + c5 + c6; - - int i1, j1, k1, l1; // The integer offsets for the second simplex corner - int i2, j2, k2, l2; // The integer offsets for the third simplex corner - int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner - - // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order. - // Many values of c will never occur, since e.g. x>y>z>w makes x=3 ? 1 : 0; - j1 = simplex[c][1]>=3 ? 1 : 0; - k1 = simplex[c][2]>=3 ? 1 : 0; - l1 = simplex[c][3]>=3 ? 1 : 0; - // The number 2 in the "simplex" array is at the second largest coordinate. - i2 = simplex[c][0]>=2 ? 1 : 0; - j2 = simplex[c][1]>=2 ? 1 : 0; - k2 = simplex[c][2]>=2 ? 1 : 0; - l2 = simplex[c][3]>=2 ? 1 : 0; - // The number 1 in the "simplex" array is at the second smallest coordinate. - i3 = simplex[c][0]>=1 ? 1 : 0; - j3 = simplex[c][1]>=1 ? 1 : 0; - k3 = simplex[c][2]>=1 ? 1 : 0; - l3 = simplex[c][3]>=1 ? 1 : 0; - // The fifth corner has all coordinate offsets = 1, so no need to look that up. - - float x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords - float y1 = y0 - j1 + G4; - float z1 = z0 - k1 + G4; - float w1 = w0 - l1 + G4; - float x2 = x0 - i2 + 2.0*G4; // Offsets for third corner in (x,y,z,w) coords - float y2 = y0 - j2 + 2.0*G4; - float z2 = z0 - k2 + 2.0*G4; - float w2 = w0 - l2 + 2.0*G4; - float x3 = x0 - i3 + 3.0*G4; // Offsets for fourth corner in (x,y,z,w) coords - float y3 = y0 - j3 + 3.0*G4; - float z3 = z0 - k3 + 3.0*G4; - float w3 = w0 - l3 + 3.0*G4; - float x4 = x0 - 1.0 + 4.0*G4; // Offsets for last corner in (x,y,z,w) coords - float y4 = y0 - 1.0 + 4.0*G4; - float z4 = z0 - 1.0 + 4.0*G4; - float w4 = w0 - 1.0 + 4.0*G4; - - // Work out the hashed gradient indices of the five simplex corners - int ii = i & 255; - int jj = j & 255; - int kk = k & 255; - int ll = l & 255; - int gi0 = perm[ii+perm[jj+perm[kk+perm[ll]]]] % 32; - int gi1 = perm[ii+i1+perm[jj+j1+perm[kk+k1+perm[ll+l1]]]] % 32; - int gi2 = perm[ii+i2+perm[jj+j2+perm[kk+k2+perm[ll+l2]]]] % 32; - int gi3 = perm[ii+i3+perm[jj+j3+perm[kk+k3+perm[ll+l3]]]] % 32; - int gi4 = perm[ii+1+perm[jj+1+perm[kk+1+perm[ll+1]]]] % 32; - - // Calculate the contribution from the five corners - float t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0; - if(t0<0) n0 = 0.0; - else { - t0 *= t0; - n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); - } - - float t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1; - if(t1<0) n1 = 0.0; - else { - t1 *= t1; - n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); - } - - float t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2; - if(t2<0) n2 = 0.0; - else { - t2 *= t2; - n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); - } - - float t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3; - if(t3<0) n3 = 0.0; - else { - t3 *= t3; - n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); - } - - float t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4; - if(t4<0) n4 = 0.0; - else { - t4 *= t4; - n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); - } - - // Sum up and scale the result to cover the range [-1,1] - return 27.0 * (n0 + n1 + n2 + n3 + n4); -} - - -int fastfloor( const float x ) { return x > 0 ? (int) x : (int) x - 1; } - -float dot( const int* g, const float x, const float y ) { return g[0]*x + g[1]*y; } -float dot( const int* g, const float x, const float y, const float z ) { return g[0]*x + g[1]*y + g[2]*z; } -float dot( const int* g, const float x, const float y, const float z, const float w ) { return g[0]*x + g[1]*y + g[2]*z + g[3]*w; } - -#define PI 3.14159265 - -int simplexnoise(int32_t seed, float* map, int width, int height, float persistence) -{ - float ka = 256/seed; - float kb = seed*567%256; - float kc = (seed*seed) % 256; - float kd = (567-seed) % 256; - for (int y = 0; y < height; y++) { - for (int x = 0; x < width; x++) { - float fNX = x/(float)width; // we let the x-offset define the circle - float fNY = y/(float)height; // we let the x-offset define the circle - float fRdx = fNX*2*PI; // a full circle is two pi radians - float fRdy = fNY*4*PI; // a full circle is two pi radians - float fYSin = sinf(fRdy); - float fRdsSin = 1.0f; - float noiseScale = 0.593; - float a = fRdsSin*sinf(fRdx); - float b = fRdsSin*cosf(fRdx); - float c = fRdsSin*sinf(fRdy); - float d = fRdsSin*cosf(fRdy); - float v = scaled_octave_noise_4d(16.0f, - persistence, - 0.5f, - 0.0f, - 1.0f, - ka+a*noiseScale, - kb+b*noiseScale, - kc+c*noiseScale, - kd+d*noiseScale); - if (map[y * width + x] == 0.0f) map[y * width + x] = v; - } - } - - return 0; -} - diff --git a/pybindings/src/sqrdmd.cpp b/pybindings/src/sqrdmd.cpp deleted file mode 100644 index ca4fd79..0000000 --- a/pybindings/src/sqrdmd.cpp +++ /dev/null @@ -1,227 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -/** @file hmapgen_sqrdmd.c - * @brief Contains functions to generate fractal height maps. - * - * @author Lauri Viitanen - * @date 2011-08-09 - */ -#include -#include -#include -#ifdef __MINGW32__ // this is to avoid a problem with the hypot function which is messed up by Python... -#undef __STRICT_ANSI__ -#endif -#include "simplerandom.hpp" - -#include "sqrdmd.hpp" - -using namespace std; - -#define CALC_SUM(a, b, c, d, rnd)\ -{\ - sum = ((a) + (b) + (c) + (d)) * 0.25f;\ - sum = sum + slope * rnd;\ -} - -#define SAVE_SUM(a)\ -{\ - bool isZero = (int)map[a] == 0; \ - if (isZero) { \ - map[a] = sum; \ - } \ -} - -void normalize(float* arr, int size) -{ - float min = arr[0], max = arr[0], diff; - - for (int i = 1; i < size; ++i) - { - min = min < arr[i] ? min : arr[i]; - max = max > arr[i] ? max : arr[i]; - } - - diff = max - min; - - if (diff > 0) - for (int i = 0; i < size; ++i) - arr[i] = (arr[i] - min) / diff; -} - -class Coord -{ -public: - Coord(int width, int height) : _width(width), _height(height) - {}; - int indexOf(int x, int y) const - { - if (x < 0 || x >= _width) throw invalid_argument("x is not valid"); - if (y < 0 || y >= _height) throw invalid_argument("y is not valid"); - return y * _width + x; - } -private: - int _width, _height; -}; - -int sqrdmd(long seed, float* map, int size, float rgh) -{ - SimpleRandom _randsource(seed); - - const int full_size = size * size; - - int i, temp; - int x, y, dx, dy; - int x0, x1, y0, y1; - int p0, p1, p2, p3; - int step, line_jump, masked; - float slope, sum, center_sum; - i = 0; - temp = size - 1; - // MUST EQUAL TO 2^x + 1! - if (temp & (temp - 1) || temp & 3) { - throw invalid_argument("Side should be 2**n +1"); - } - temp = size; - slope = rgh; - step = size & ~1; - - /* Calculate midpoint ("diamond step"). */ - dy = step * size; - CALC_SUM(map[0], map[step], map[dy], map[dy + step], _randsource.next_float_signed()); - SAVE_SUM(i); - center_sum = sum; - - /* Calculate each sub diamonds' center points ("square step"). */ - /* Top row. */ - p0 = step >> 1; - CALC_SUM(map[0], map[step], center_sum, center_sum, _randsource.next_float_signed()); - SAVE_SUM(p0); - /* Left column. */ - p1 = p0 * size; - CALC_SUM(map[0], map[dy], center_sum, center_sum, _randsource.next_float_signed()); - SAVE_SUM(p1); - map[full_size + p0 - size] = map[p0]; /* Copy top val into btm row. */ - map[p1 + size - 1] = map[p1]; /* Copy left value into right column. */ - slope *= rgh; - step >>= 1; - - while (step > 1) /* Enter the main loop. */ - { - /************************************************************* - * Calc midpoint of sub squares on the map ("diamond step"). * - *************************************************************/ - dx = step; - dy = step * size; - i = (step >> 1) * (size + 1); - line_jump = step * size + 1 + step - size; - for (y0 = 0, y1 = dy; y1 < size * size; y0 += dy, y1 += dy) - { - for (x0 = 0, x1 = dx; x1 < size; x0 += dx, x1 += dx, i += step) - { - sum = (map[y0+x0] + map[y0+x1] + - map[y1+x0] + map[y1+x1]) * 0.25f; - sum = sum + slope * _randsource.next_float_signed(); - masked = !((int)map[i]); - map[i] = map[i] * !masked + sum * masked; - } - /* There's additional step taken at the end of last - * valid loop. That step actually isn't valid because - * the row ends right then. Thus we are forced to - * manually remove it after the loop so that 'i' - * points again to the index accessed last. - */ - i += line_jump - step; - } - - /************************************************************** - * Calculate each sub diamonds' center point ("square step"). - * Diamond gets its left and right vertices from the square - * corners of last iteration and its top and bottom vertices - * from the "diamond step" we just performed. - *************************************************************/ - i = step >> 1; - p0 = step; /* right */ - p1 = i * size + i; /* bottom */ - p2 = 0; /* left */ - p3 = full_size + i - (i + 1) * size; /* top (wrapping edges) */ - - /* Calculate "diamond" values for top row in map. */ - while (p0 < size) - { - sum = (map[p0] + map[p1] + map[p2] + map[p3]) * 0.25f; - sum = sum + slope * _randsource.next_float_signed(); - masked = !((int)map[i]); - map[i] = map[i] * !masked + sum * masked; - /* Copy it into bottom row. */ - map[full_size + i - size] = map[i]; - p0 += step; p1 += step; p2 += step; - p3 += step; i += step; - } - - /* Now that top row's values are calculated starting from - * 'y = step >> 1' both saves us from recalculating same things - * twice and guarantees that data will not be read beyond top - * row of map. 'size - (step >> 1)' guarantees that data will - * not be read beyond bottom row of map. - */ - for (y = step >> 1, temp = 0; y < size - (step >> 1); y += step >> 1, temp = !temp) - { - p0 = step >> 1; /* right */ - p1 = p0 * size; /* bottom */ - p2 = -p0; /* left */ - p3 = -p1; /* top */ - /* For even rows add step/2. Otherwise add nothing. */ - x = i = p0 * temp; /* Init 'x' while it's easy. */ - i += y * size; /* Move 'i' into correct row. */ - p0 += i; - p1 += i; - /* For odd rows p2 (left) wraps around map edges. */ - p2 += i + (size - 1) * !temp; - p3 += i; - /* size - (step >> 1) guarantees that data will not be - * read beyond rightmost column of map. */ - for (; x < size - (step >> 1); x += step) - { - sum = (map[p0] + map[p1] + - map[p2] + map[p3]) * 0.25f; - sum = sum + slope * _randsource.next_float_signed(); - masked = !((int)map[i]); - map[i] = map[i] * !masked + sum * masked; - p0 += step; - p1 += step; - p2 += step; - p3 += step; - i += step; - /* if we start from leftmost column -> left - * point (p2) is going over the right border -> - * wrap it around into the beginning of - * previous rows left line. */ - p2 -= (size - 1) * !x; - } - /* copy rows first element into its last */ - i = y * size; - map[i + size - 1] = map[i]; - } - slope *= rgh; /* reduce amount of randomness for next round */ - step >>= 1; /* split squares and diamonds in half */ - } - return (0); -} diff --git a/pybindings/src/utils.cpp b/pybindings/src/utils.cpp deleted file mode 100644 index c757629..0000000 --- a/pybindings/src/utils.cpp +++ /dev/null @@ -1,51 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "utils.hpp" -#include -#include -#include - -namespace Platec { - -std::string to_string(uint32_t value) -{ - std::stringstream ss; - ss << value; - std::string str = ss.str(); - return str; -} - -std::string to_string_f(float value) -{ - std::stringstream ss; - ss << value; - std::string str = ss.str(); - return str; -} - -} - -void p_assert(bool condition, const std::string& message) -{ - if (!condition) { - std::cerr << "Assertion failed: " << message << std::endl; - exit(1); - } -} \ No newline at end of file diff --git a/pybindings/src/world_point.cpp b/pybindings/src/world_point.cpp deleted file mode 100644 index 8c73077..0000000 --- a/pybindings/src/world_point.cpp +++ /dev/null @@ -1,58 +0,0 @@ -/****************************************************************************** - * plate-tectonics, a plate tectonics simulation library - * Copyright (C) 2012-2013 Lauri Viitanen - * Copyright (C) 2014-2015 Federico Tomassetti, Bret Curtis - * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, see http://www.gnu.org/licenses/ - *****************************************************************************/ - -#include "world_point.hpp" -#include "rectangle.hpp" - -WorldPoint::WorldPoint(uint32_t x, uint32_t y, const WorldDimension& worldDimension) - : _x(x), _y(y) -{ - if (_x >= worldDimension.getWidth()) { - throw runtime_error("WorldPoint::WorldPoint"); - } - if (_y >= worldDimension.getHeight()) { - throw runtime_error("WorldPoint::WorldPoint"); - } -} - -WorldPoint::WorldPoint(const WorldPoint& other) - : _x(other.x()), - _y(other.y()) -{ } - -uint32_t WorldPoint::x() const -{ - return _x; -} - -uint32_t WorldPoint::y() const -{ - return _y; -} - -uint32_t WorldPoint::toIndex(const WorldDimension& dim) const -{ - if (_x >= dim.getWidth()) { - throw runtime_error("WorldPoint::WorldPoint"); - } - if (_y >= dim.getHeight()) { - throw runtime_error("WorldPoint::WorldPoint"); - } - return _y * dim.getWidth() + _x; -}