diff --git a/src/ComplexNets/GraphGenerator.cpp b/src/ComplexNets/GraphGenerator.cpp index 758160e..d5ff283 100644 --- a/src/ComplexNets/GraphGenerator.cpp +++ b/src/ComplexNets/GraphGenerator.cpp @@ -13,12 +13,6 @@ using namespace graphpp; -typedef struct Position -{ - float x; - float y; -} Position; - static const double PI = atan(1) * 4; std::vector vertexesPositions; @@ -175,7 +169,7 @@ Graph* GraphGenerator::generateBarabasiAlbertGraph(unsigned int m_0, unsigned in * minimizes the euclidean distance from U to V and the hoops distance from U to root. This process * is repeated m times for each new vertex added. * Step 3) Create q new edges on the graph for the vertex added - * Step 4) A new root is chosen with probability dependant on the node degree every t rounds. + * Step 4) A new root is chosen with probability dependant on the node degree every t rounds. (@TODO) * */ Graph* GraphGenerator::generateHotExtendedGraph( @@ -183,218 +177,162 @@ Graph* GraphGenerator::generateHotExtendedGraph( { auto graph = new Graph(); - // Every time an edge is added, two entries are added in vertexIndexes, the two indexes of the - // nodes joined. - std::vector vertexIndexes; - unsigned int root = 1; + // add root vertex + Vertex* root = new Vertex(0); + graph->addVertex(root); - vertexesPositions.clear(); - // Firts vertex + addPosition(root); + + // for each of the nodes we will be adding to the graph + for (unsigned i = 1; i < n; i++) + { + Vertex* v = new Vertex(i); + addPosition(v); - // Step 1 + // setup minWeight, minVertex, minHop + auto minWeight = distance(v, root); + auto minVertex = root; + auto minHop = 0; - addOriginalVertex(graph); + // TODO: where is maxHop updated? + double maxHop = 0; - // For each of the nodes we will be adding to the graph - for (unsigned int i = 2; i <= n; i++) - { - // Step 2 + // for hop j from 1 to maxHop + for (unsigned hop = 1; hop <= maxHop; hop++) + { + auto nn = singleHopNN(hop); + insertNN(nn, v); - addFKPNode(i, graph, root, xi, &vertexIndexes, m); + auto neighbor = getNN(graph, nn, v); + auto weight = distance(v, neighbor) + xi * hop; - // Step 3 + if (weight < minWeight) + { + minWeight = weight; + minVertex = neighbor; + minHop = hop; + } - addExtendedEdges(q, i, graph, root, r, &vertexIndexes); + removeNN(nn, v); + } - // Step 4 + auto nn = singleHopNN(minHop + 1); + insertNN(nn, v); - root = chooseNewRoot(i, t, root, vertexIndexes); - } - return graph; -} + // TODO: update Ti to Tr? + graph->addEdge(v, minVertex); -/** - The main Vertex is created, added to the graph and a position is generated -*/ -void GraphGenerator::addOriginalVertex(Graph* graph) -{ - Vertex* newVertex = new Vertex(1); - graph->addVertex(newVertex); - addVertexPosition(); -} + nn = doubleHopNN(minHop + 1, minHop + 3); + insertNN(nn, v); -/** - A new vertex is created and added. Then m edges are added according to the original FKP - function - - TODO: Optimize this function using Voronoi Diagrams -*/ -void GraphGenerator::addFKPNode( - unsigned int vertexIndex, - Graph* graph, - unsigned int root, - float xi, - std::vector* vertexIndexes, - unsigned int m) -{ - std::map distance; - - // Creation of vertex - Vertex* newVertex = new Vertex(vertexIndex); - graph->addVertex(newVertex); - addVertexPosition(); - - for (unsigned int j = 1; j < vertexIndex; j++) - { // this for evaluated "w" for each vertex already in the graph - unsigned hopsDistance = graph->hops( - graph->getVertexById(j), - graph->getVertexById(root) - ); // distance between vertex evaluated and root vertex - float euclidianDistance = distanceBetweenVertex(j, vertexIndex); // Distance between vertex evaluated and new vertex - - // Original FKP Algorithm chooses a new connection between the new vertex and the one with - // minimum W - - float w = euclidianDistance + xi * hopsDistance; - distance[w] = j; // stores 'w' as key of a map - } - addEdges( - graph, newVertex, distance, m, - vertexIndexes); // m edges are added acording to the minimum distances - distance.clear(); -} + nn = doubleHopNN(minHop - 1, minHop + 1); + insertNN(nn, v); -/** - q edges are added according to the function in the paper - - TODO: Optimize the complexity of the algorithm using Voronoi diagrams. -*/ -void GraphGenerator::addExtendedEdges( - unsigned int q, - unsigned int vertexIndex, - Graph* graph, - unsigned int root, - float r, - std::vector* vertexIndexes) -{ - // We will go through this function q times, adding q edges - for (unsigned int qfinal = 0; qfinal < q; qfinal++) - { - // At the end of the loop we only add one edge, so we save the minimum distance and the - // indexes of the nodes - float minDist = 0; - unsigned int finalJ = 0; - unsigned int finalK = 0; - - // We loop through every combination of nodes - for (unsigned int j = 1; j <= vertexIndex; j++) + std::pair minEdge; + minWeight = INFINITY; + minVertex = nullptr; + + v = nullptr; + + for (unsigned hop = 0; hop <= maxHop - 2; hop++) { - for (unsigned int k = 1; k <= vertexIndex; k++) + for (auto const& k : singleHopNN(hop)) { - // For each combination of nodes we calculate the function - float euclidianDistance = 0; - unsigned int HopsWithEdge = 0; - unsigned int HopsWithOutEdge = 0; - if (k != j && !graph->getVertexById(j)->isNeighbourOf(graph->getVertexById(k))) + nn = doubleHopNN(hop, hop + 2); + v = getNN(graph, nn, k); + + auto weight = distance(k, v) - (r / n) * childCount(v); + if (weight < minWeight) { - for (unsigned int l = 1; l <= vertexIndex; l++) - { - // We calculate the hops to the root with and without the edge - HopsWithOutEdge = - graph->hops(graph->getVertexById(l), graph->getVertexById(root)) + - HopsWithOutEdge; // Hops between evaluated vertex and root vertex - // without new edge - graph->addEdge( - graph->getVertexById(j), - graph->getVertexById( - k)); // Add new edges only for evaluation, later will be removed - HopsWithEdge = - graph->hops(graph->getVertexById(l), graph->getVertexById(root)) + - HopsWithEdge; // Hops between evaluated vertex and root vertex with new - // edge - graph->removeEdge( - graph->getVertexById(j), - graph->getVertexById(k)); // remove the edge previously added - } - - euclidianDistance = distanceBetweenVertex(j, k); - float w = - euclidianDistance + (r / vertexIndex) * (HopsWithEdge - HopsWithOutEdge); - - // We save only the minimum - - if (minDist == 0 || w < minDist) - { - minDist = w; - finalJ = j; - finalK = k; - } + minWeight = weight; + minEdge = std::make_pair(k, v); + minVertex = v; } } } - // Finally, we add a new edge + graph->addEdge(minEdge.first, minEdge.second); - if (minDist != 0) + auto hp = 0; + + // TODO: see whiteboard + for (auto const& j : BFS(root)) { - vertexIndexes->push_back(finalJ); - vertexIndexes->push_back(finalK); - graph->addEdge(graph->getVertexById(finalJ), graph->getVertexById(finalK)); + auto hj = weightChange(j); + if (hj > hp) + { + nn = doubleHopNN(hj - 2, hj); + removeNN(nn, j); + + nn = doubleHopNN(hj, hj + 2); + removeNN(nn, j); + + nn = singleHopNN(hj); + removeNN(nn, j); + + nn = doubleHopNN(hp - 2, hp); + insertNN(nn, j); + + nn = doubleHopNN(hp, hp + 2); + insertNN(nn, j); + + nn = singleHopNN(hp); + insertNN(nn, j); + } } } + return graph; } -/** - A new root is choosen according to the parameter t and the indexes distribution -*/ -int GraphGenerator::chooseNewRoot( - unsigned int vertexIndex, - unsigned int t, - unsigned int root, - std::vector vertexIndexes) +void GraphGenerator::addPosition(Vertex* v) { - if ((vertexIndex - 1) % t == 0) - { - return vertexIndexes[rand() % vertexIndexes.size()]; - } - return root; + // assert(vertex.id == vertexesPositions.size()) + Position p = { + (float)rand() / RAND_MAX, + (float)rand() / RAND_MAX + }; + vertexesPositions.push_back(p); } -void GraphGenerator::addEdges( - Graph* graph, - Vertex* vertex, - std::map distance, - unsigned int q, - std::vector* vertexIndexes) +Position GraphGenerator::position(Vertex* v) { - for (unsigned int k = 0; k < q && !distance.empty(); k++) - { // Adding "q" new edges. The processes is similar to added vertex. - if (!graph->getVertexById(distance.begin()->second)->isNeighbourOf(vertex)) - { - vertexIndexes->push_back(distance.begin()->second); - vertexIndexes->push_back(vertex->getVertexId()); - graph->addEdge(graph->getVertexById(distance.begin()->second), vertex); - distance.erase(distance.begin()); // Remove the lowest value of w from the list, - // because was used previously. This process repeat - // "q" times. - if (distance.empty()) - break; - } - } + return vertexesPositions[v->getVertexId()]; +} + +double GraphGenerator::distance(Vertex* v1, Vertex* v2) +{ + Position p1 = position(v1); + Position p2 = position(v2); + + return sqrt(pow(p2.x - p1.x, 2) + pow(p2.y - p1.y, 2)); +} + +void GraphGenerator::insertNN(NearestNeighbor* nn, Vertex* v) +{ + nn->insert(v->getVertexId(), position(v)); +} + +void GraphGenerator::removeNN(NearestNeighbor* nn, Vertex* v) +{ + nn->remove(v->getVertexId()); +} + +Vertex* GraphGenerator::getNN(Graph *g, NearestNeighbor* nn, Vertex* v) +{ + auto id = nn->get(v->getVertexId()); + return g->getVertexById(id); } -float GraphGenerator::distanceBetweenVertex(unsigned int vertex1Id, unsigned int vertex2Id) +double GraphGenerator::weightChange(Vertex* v) { - return sqrt( - pow(vertexesPositions[vertex1Id - 1].x - vertexesPositions[vertex2Id - 1].x, 2) + - pow(vertexesPositions[vertex1Id - 1].y - vertexesPositions[vertex2Id - 1].y, 2)); + // TODO: stub + return 0; } -void GraphGenerator::addVertexPosition() +unsigned GraphGenerator::childCount(Vertex* v) { - vertexesPositions.push_back(Position()); - vertexesPositions.back().x = (float)rand() / RAND_MAX; - vertexesPositions.back().y = (float)rand() / RAND_MAX; + // TODO: stub + return 0; } Graph* GraphGenerator::generateMolloyReedGraph(std::string path) @@ -413,7 +351,7 @@ Graph* GraphGenerator::generateMolloyReedGraph(std::string path) * Computes the distance in a hiperbolic space between two points */ -inline double GraphGenerator::hiperbolicDistance(PolarPosition p1, PolarPosition p2) +double GraphGenerator::hiperbolicDistance(PolarPosition p1, PolarPosition p2) { return acosh( cosh(p1.r) * cosh(p2.r) - @@ -423,7 +361,7 @@ inline double GraphGenerator::hiperbolicDistance(PolarPosition p1, PolarPosition /* * Computes random polar hyperbolic coordinates */ -inline GraphGenerator::PolarPosition GraphGenerator::getRandomHyperbolicCoordinates( +PolarPosition GraphGenerator::getRandomHyperbolicCoordinates( float a, double maxr) { PolarPosition pos; diff --git a/src/ComplexNets/GraphGenerator.h b/src/ComplexNets/GraphGenerator.h index 0a3986c..ff5d14a 100644 --- a/src/ComplexNets/GraphGenerator.h +++ b/src/ComplexNets/GraphGenerator.h @@ -5,54 +5,38 @@ #pragma once #include "ComplexNets/MolloyReedGraphReader.h" +#include "ComplexNets/NearestNeighbor.h" +#include "ComplexNets/Position.h" #include "ComplexNets/typedefs.h" class GraphGenerator { private: - typedef struct PolarPosition - { - double r; - double theta; - } PolarPosition; - GraphGenerator(); static GraphGenerator* instance; - float distanceBetweenVertex(unsigned int vertex1Id, unsigned int vertex2Id); - void addVertexPosition(); - void addEdges( - Graph* graph, - Vertex* vertex, - std::map distance, - unsigned int quant, - std::vector* vertexIndexes); - inline double hiperbolicDistance(PolarPosition p1, PolarPosition p2); - inline double getMaxRadius(int i, float a, float c); - inline PolarPosition getRandomHyperbolicCoordinates(float a, double maxr); + + // members related to hyperbolic graphs + double hiperbolicDistance(PolarPosition p1, PolarPosition p2); + double getMaxRadius(int i, float a, float c); + PolarPosition getRandomHyperbolicCoordinates(float a, double maxr); + + // members related to HOT Extended graphs + void addPosition(Vertex* v); + NearestNeighbor* singleHopNN(unsigned hop); + NearestNeighbor* doubleHopNN(unsigned hop1, unsigned hop2); + + void insertNN(NearestNeighbor* nn, Vertex *v); + void removeNN(NearestNeighbor* nn, Vertex *v); + Vertex* getNN(Graph* g, NearestNeighbor* nn, Vertex *v); + + double distance(Vertex* v1, Vertex* v2); + double weightChange(Vertex* v); + unsigned childCount(Vertex* v); + Position position(Vertex* v); public: static GraphGenerator* getInstance(); - void addOriginalVertex(Graph* graph); - void addFKPNode( - unsigned int vertexIndex, - Graph* graph, - unsigned int root, - float xi, - std::vector* vertexIndexes, - unsigned int m); - void addExtendedEdges( - unsigned int q, - unsigned int vertexIndex, - Graph* graph, - unsigned int root, - float r, - std::vector* vertexIndexes); - int chooseNewRoot( - unsigned int vertexIndex, - unsigned int t, - unsigned int root, - std::vector vertexIndexes); Graph* generateGraphFromFile(std::string path, bool directed, bool multigraph); DirectedGraph* generateDirectedGraphFromFile(std::string path, bool multigraph); diff --git a/src/ComplexNets/NearestNeighbor.cpp b/src/ComplexNets/NearestNeighbor.cpp new file mode 100644 index 0000000..c086ff9 --- /dev/null +++ b/src/ComplexNets/NearestNeighbor.cpp @@ -0,0 +1,21 @@ +// This toolbox is licensed under the Academic Free License 3.0. +// Instituto Tecnológico de Buenos Aires (ITBA). + +#include "NearestNeighbor.h" + + +void NearestNeighbor::insert(unsigned id, Position pos) +{ + // TODO: stub +} + +void NearestNeighbor::remove(unsigned id) +{ + // TODO: stub +} + +unsigned NearestNeighbor::insert(unsigned id, Position pos) +{ + // TODO: stub + return 0; +} diff --git a/src/ComplexNets/NearestNeighbor.h b/src/ComplexNets/NearestNeighbor.h new file mode 100644 index 0000000..1413ffa --- /dev/null +++ b/src/ComplexNets/NearestNeighbor.h @@ -0,0 +1,14 @@ +// This toolbox is licensed under the Academic Free License 3.0. +// Instituto Tecnológico de Buenos Aires (ITBA). + +#pragma once + +#include "ComplexNets/Position.h" + +class NearestNeighbor +{ +public: + void insert(unsigned id, Position pos); + void remove(unsigned id); + unsigned get(unsigned id); +}; diff --git a/src/ComplexNets/Position.h b/src/ComplexNets/Position.h new file mode 100644 index 0000000..e037be3 --- /dev/null +++ b/src/ComplexNets/Position.h @@ -0,0 +1,17 @@ +// This toolbox is licensed under the Academic Free License 3.0. +// Instituto Tecnológico de Buenos Aires (ITBA). + +#pragma once + +typedef struct Position +{ + // TODO: change to double? + float x; + float y; +} Position; + +typedef struct PolarPosition +{ + double r; + double theta; +} PolarPosition; diff --git a/src/ComplexNetsGui.pro b/src/ComplexNetsGui.pro index 71a0ae1..59285b7 100644 --- a/src/ComplexNetsGui.pro +++ b/src/ComplexNetsGui.pro @@ -8,6 +8,7 @@ SOURCES += ComplexNetsGui/src/*.cpp \ ComplexNets/*.cpp \ ComplexNetsCmd/*.cpp \ ComplexNetsCmd/*.c + vendor/ HEADERS += ComplexNetsGui/inc/*.h \ ComplexNets/*.h \ ComplexNetsCmd/*.h diff --git a/src/ComplexNetsGui/inc/ui_mainwindow.h b/src/ComplexNetsGui/inc/ui_mainwindow.h index 07c826f..e139fb4 100644 --- a/src/ComplexNetsGui/inc/ui_mainwindow.h +++ b/src/ComplexNetsGui/inc/ui_mainwindow.h @@ -1,7 +1,7 @@ /******************************************************************************** ** Form generated from reading UI file 'mainwindow.ui' ** -** Created by: Qt User Interface Compiler version 5.9.0 +** Created by: Qt User Interface Compiler version 5.9.1 ** ** WARNING! All changes made in this file will be lost when recompiling UI file! ********************************************************************************/ diff --git a/src/__init__.py b/src/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/optimized_hot.py b/src/optimized_hot.py new file mode 100644 index 0000000..a07380c --- /dev/null +++ b/src/optimized_hot.py @@ -0,0 +1,222 @@ +# coding: utf-8 +from collections import defaultdict +from random import random +import math +import sys + +import logging as l +l.basicConfig(format="%(message)s", level=l.DEBUG) + + + +class Vertex(object): + + def __init__(self, id): + self.id = id + self.position = (random(), random()) + + def __repr__(self): + return "Vertex<%s>(%s, %s)" % (self.id, self.position[0], self.position[1]) + + +class Graph(object): + + def __init__(self): + self.vertices = {} + + def __repr__(self): + return "Graph(%s)" % self.vertices + + def add_vertex(self, v): + self.vertices[v] = [] + + def add_edge(self, v1, v2): + self.vertices[v1].append(v2) + self.vertices[v2].append(v1) + + def breadth_first_search(self, root): + visited = set() + q = [root] + + while q: + v = q.pop(0) + visited.add(v) + for j in self.vertices[v]: + if j not in visited: + q.append(j) + yield v + + +def export(graph, file_obj): + for v in graph.vertices: + for j in graph.vertices[v]: + file_obj.write("%s\t%s\n" % (v.id, j.id)) + + +# HOT Extended Model +# The Paper can be checked in: +# http://cnet.fi.uba.ar/ignacio.alvarez-hamelin/pdf/model_internet_jiah_ns.pdf +# +# m: +# number of edges in each new vertex (k in the paper) +# n: +# total number of nodes +# xi: +# parameter used to select the neighbors for new vertex. (theta in the paper) +# q: +# number of edges added in the graph after of connect a vertex. (q in the paper) +# r: +# parameter used to select the edges in the graph after connecting a vertex (γ in the paper) +# t: +# parameter that represents how many rounds to do until a new root selection (τ in the paper) + +def hot_extended(m, n, xi, q, r, t): + + graph = Graph() + + # first we add the root vertex + root = Vertex(0) + graph.add_vertex(root) + + # we then proceed to add the n-1 remaining nodes + for i in range(1, n): + + # we add the new vertex with a random position + v = Vertex(i) + graph.add_vertex(v) + + w_min = distance(v, root) + v_min = root + h_min = 0 + + h_max = 0; + + # for hop j from 1 to h_max + for j in range(1, h_max): + + # we insert our new vertex in the delaunay for each hop + d_j = delaunay(j) + d_j.insert(v) + + # we get the nearest neighbor of v in this hop + neighbor = d_j.nearest_neighbor(v) + + # we compute the weighted distance to this neighbor + w = distance(v, neighbor) + xi * j + + # if this is the nearest neighbor up to this hop, update the min neighbor + if w < w_min: + w_min = w + v_min = neighbor + h_min = j + + # leave delaunay for hop as we without changes + d_j.remove(v) + + # we will connect v to a neighbor that is at hop h_min, so add it to the next hop delaunay + delaunay(h_min + 1).insert(v) + + # connect v with its nearest weighted neighbor + # @question: why do we add a single edge here? + graph.add_edge(v, v_min) + + # this should end the step 1 described in the paper + # @question: shouldn't we add m new edges (q in the paper)? + + # now v should be available in the delaunay for hops pairs 1+ and 1- + delaunay(h_min + 1, h_min + 3).insert(v) + delaunay(h_min - 1, h_min + 1).insert(v) + + # this should be the start of step 2 described in the paper + e_min = None + w_min = float('inf') + v_min = None + + for j in range(0, h_max - 2): + + for k in delaunay(j): + + d_ = delaunay(hop, hop + 2) + v = d_.nearest_neighbor(k) + + # @question: how do we compute child count here? we need to mark the parent edge somehow to skip it + w = distance(k, v) - (r / n) * child_count(v); + + if w < w_min: + w_min = w + v_min = v + e_min = (k, v) + + # graph.add_edge(*e_min); + # @question: why do we add a single edge here instead of q? + + hp = 0 + + for j in graph.breadth_first_search(root): + # @question: is this where we compute T(j)? + + hj = weight_change(j); + + if (hj > hp): + delaunay(hj - 2, hj).remove(j) + delaunay(hj, hj + 2).remove(j) + delaunay(hj).remove(j) + + delaunay(hp - 2, hp).insert(j) + delaunay(hp, hp + 2).insert(j) + delaunay(hp).insert(j) + + # @question: shouldn't this be hp = hj? + hj = hp + + # @question: original algorithm should change root from time to time? this changes hop calculation + + return graph + +def distance(v1, v2): + return math.hypot(v2.position[0] - v1.position[0], v2.position[1] - v1.position[1]) + +def weight_change(v): + return 0 # stub + +def child_count(v): + return float('inf') # stub + + +class Delaunay(object): # stub for now + + def __init__(self): + self.vertices = [] + + def insert(self, v): + self.vertices.append(v) + + def remove(self, v): + self.vertices.remove(v) + + def nearest_neighbor(self, v): + min_d = float('inf') + min_w = None + + for w in self.vertices: + if v == w: + continue + + d = distance(v, w) + if d < min_d: + min_d = d + min_w = w + + return min_w + + +delaunay_map = defaultdict(Delaunay) + +def delaunay(hop, hop_s=None): + key = (hop, hop_s) + return delaunay_map[key] + + +if __name__ == "__main__": + g = hot_extended(2, 10, 0.3, 2, 0.2, 0.1) + # export(g, sys.stdout) diff --git a/vendor/delaunay/.gitignore b/vendor/delaunay/.gitignore new file mode 100644 index 0000000..7563863 --- /dev/null +++ b/vendor/delaunay/.gitignore @@ -0,0 +1,7 @@ +test +delaunay +*.swp +build/ +peda-session-delaunay.txt +*.lua +.gdb_history diff --git a/vendor/delaunay/LICENSE b/vendor/delaunay/LICENSE new file mode 100644 index 0000000..71568cf --- /dev/null +++ b/vendor/delaunay/LICENSE @@ -0,0 +1,24 @@ +Copyright (c) 2015 Simon Zeni (simonzeni@gmail.com) + + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. + + diff --git a/vendor/delaunay/README.md b/vendor/delaunay/README.md new file mode 100644 index 0000000..c0c6c85 --- /dev/null +++ b/vendor/delaunay/README.md @@ -0,0 +1,62 @@ +# delaunay-triangulation + +## Pseudo-code algorithm + +``` +function BowyerWatson (pointList) + // pointList is a set of coordinates defining the points to be triangulated + triangulation := empty triangle mesh data structure + add super-triangle to triangulation // must be large enough to completely contain all the points in pointList + for each point in pointList do // add all the points one at a time to the triangulation + badTriangles := empty set + for each triangle in triangulation do // first find all the triangles that are no longer valid due to the insertion + if point is inside circumcircle of triangle + add triangle to badTriangles + polygon := empty set + for each triangle in badTriangles do // find the boundary of the polygonal hole + for each edge in triangle do + if edge is not shared by any other triangles in badTriangles + add edge to polygon + for each triangle in badTriangles do // remove them from the data structure + remove triangle from triangulation + for each edge in polygon do // re-triangulate the polygonal hole + newTri := form a triangle from edge to point + add newTri to triangulation + for each triangle in triangulation // done inserting points, now clean up + if triangle contains a vertex from original super-triangle + remove triangle from triangulation + return triangulation +``` + +## Sample + +![alt text](https://github.com/Bl4ckb0ne/delaunay-triangulation/blob/master/sample.png "Sample image (if you see this, then the image can't load or hasn't loaded yet)") + + +From the [Wikipedia page of the algorithm](https://en.wikipedia.org/wiki/Bowyer%E2%80%93Watson_algorithm) + +## Requirement + +You will need [SFML 2+](http://www.sfml-dev.org/download/sfml/2.3.2/) to run the example, and C++11 to compile it. + +## Usage + +To build it, you can type in : +```sh +make +``` +You may change the compiler on the makefile (using the CXX var) +```sh +make CXX=g++ # to use the GCC compiler +make CXX=clang++ # default compiler +``` + +The executable name is ``` delaunay ```, without arguments +```sh +./delaunay +``` + +You also can clear the executable and the build folder. +```sh +make clean +``` diff --git a/vendor/delaunay/delaunay.h b/vendor/delaunay/delaunay.h new file mode 100644 index 0000000..e04ce23 --- /dev/null +++ b/vendor/delaunay/delaunay.h @@ -0,0 +1,145 @@ +#ifndef H_DELAUNAY +#define H_DELAUNAY + +#include "vector2.h" +#include "edge.h" +#include "triangle.h" + +#include +#include + +template +class Delaunay +{ + public: + using TriangleType = Triangle; + using EdgeType = Edge; + using VertexType = Vector2; + + const std::vector& triangulate(std::vector &vertices) + { + // Store the vertices localy + _vertices = vertices; + + // Determinate the super triangle + float minX = vertices[0].x; + float minY = vertices[0].y; + float maxX = minX; + float maxY = minY; + + for(std::size_t i = 0; i < vertices.size(); ++i) + { + if (vertices[i].x < minX) minX = vertices[i].x; + if (vertices[i].y < minY) minY = vertices[i].y; + if (vertices[i].x > maxX) maxX = vertices[i].x; + if (vertices[i].y > maxY) maxY = vertices[i].y; + } + + float dx = maxX - minX; + float dy = maxY - minY; + float deltaMax = std::max(dx, dy); + float midx = (minX + maxX) / 2.f; + float midy = (minY + maxY) / 2.f; + + VertexType p1(midx - 20 * deltaMax, midy - deltaMax); + VertexType p2(midx, midy + 20 * deltaMax); + VertexType p3(midx + 20 * deltaMax, midy - deltaMax); + + //std::cout << "Super triangle " << std::endl << Triangle(p1, p2, p3) << std::endl; + + // Create a list of triangles, and add the supertriangle in it + _triangles.push_back(TriangleType(p1, p2, p3)); + + for(auto p = begin(vertices); p != end(vertices); p++) + { + //std::cout << "Traitement du point " << *p << std::endl; + //std::cout << "_triangles contains " << _triangles.size() << " elements" << std::endl; + + std::vector badTriangles; + std::vector polygon; + + for(auto t = begin(_triangles); t != end(_triangles); t++) + { + //std::cout << "Processing " << std::endl << *t << std::endl; + + if(t->circumCircleContains(*p)) + { + //std::cout << "Pushing bad triangle " << *t << std::endl; + badTriangles.push_back(*t); + polygon.push_back(t->e1); + polygon.push_back(t->e2); + polygon.push_back(t->e3); + } + else + { + //std::cout << " does not contains " << *p << " in his circum center" << std::endl; + } + } + + _triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [badTriangles](TriangleType &t){ + for(auto bt = begin(badTriangles); bt != end(badTriangles); bt++) + { + if(*bt == t) + { + //std::cout << "Removing bad triangle " << std::endl << *bt << " from _triangles" << std::endl; + return true; + } + } + return false; + }), end(_triangles)); + + std::vector badEdges; + for(auto e1 = begin(polygon); e1 != end(polygon); e1++) + { + for(auto e2 = begin(polygon); e2 != end(polygon); e2++) + { + if(e1 == e2) + continue; + + if(*e1 == *e2) + { + badEdges.push_back(*e1); + badEdges.push_back(*e2); + } + } + } + + polygon.erase(std::remove_if(begin(polygon), end(polygon), [badEdges](EdgeType &e){ + for(auto it = begin(badEdges); it != end(badEdges); it++) + { + if(*it == e) + return true; + } + return false; + }), end(polygon)); + + for(auto e = begin(polygon); e != end(polygon); e++) + _triangles.push_back(TriangleType(e->p1, e->p2, *p)); + + } + + _triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [p1, p2, p3](TriangleType &t){ + return t.containsVertex(p1) || t.containsVertex(p2) || t.containsVertex(p3); + }), end(_triangles)); + + for(auto t = begin(_triangles); t != end(_triangles); t++) + { + _edges.push_back(t->e1); + _edges.push_back(t->e2); + _edges.push_back(t->e3); + } + + return _triangles; + } + + const std::vector& getTriangles() const { return _triangles; }; + const std::vector& getEdges() const { return _edges; }; + const std::vector& getVertices() const { return _vertices; }; + + private: + std::vector _triangles; + std::vector _edges; + std::vector _vertices; +}; + +#endif diff --git a/vendor/delaunay/edge.h b/vendor/delaunay/edge.h new file mode 100644 index 0000000..cc9e8ac --- /dev/null +++ b/vendor/delaunay/edge.h @@ -0,0 +1,33 @@ +#ifndef H_EDGE +#define H_EDGE + +#include "vector2.h" + +template +class Edge +{ + public: + using VertexType = Vector2; + + Edge(const VertexType &p1, const VertexType &p2) : p1(p1), p2(p2) {}; + Edge(const Edge &e) : p1(e.p1), p2(e.p2) {}; + + VertexType p1; + VertexType p2; +}; + +template +inline std::ostream &operator << (std::ostream &str, Edge const &e) +{ + return str << "Edge " << e.p1 << ", " << e.p2; +} + +template +inline bool operator == (const Edge & e1, const Edge & e2) +{ + return (e1.p1 == e2.p1 && e1.p2 == e2.p2) || + (e1.p1 == e2.p2 && e1.p2 == e2.p1); +} + +#endif + diff --git a/vendor/delaunay/main.cpp b/vendor/delaunay/main.cpp new file mode 100644 index 0000000..0f1f649 --- /dev/null +++ b/vendor/delaunay/main.cpp @@ -0,0 +1,99 @@ +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "vector2.h" +#include "triangle.h" +#include "delaunay.h" + +float RandomFloat(float a, float b) { + float random = ((float) rand()) / (float) RAND_MAX; + float diff = b - a; + float r = random * diff; + return a + r; +} + +int main() +{ + srand (time(NULL)); + float numberPoints = roundf(RandomFloat(4, 40)); + + std::cout << "Generating " << numberPoints << " random points" << std::endl; + + std::vector> points; + for(int i = 0; i < numberPoints; i++) { + points.push_back(Vector2(RandomFloat(0, 800), RandomFloat(0, 600))); + } + + Delaunay triangulation; + std::vector> triangles = triangulation.triangulate(points); + std::cout << triangles.size() << " triangles generated\n"; + std::vector> edges = triangulation.getEdges(); + + std::cout << " ========= "; + + std::cout << "\nPoints : " << points.size() << std::endl; + for(auto &p : points) + std::cout << p << std::endl; + + std::cout << "\nTriangles : " << triangles.size() << std::endl; + for(auto &t : triangles) + std::cout << t << std::endl; + + std::cout << "\nEdges : " << edges.size() << std::endl; + for(auto &e : edges) + std::cout << e << std::endl; + + // SFML window + sf::RenderWindow window(sf::VideoMode(800, 600), "Delaunay triangulation"); + + // Transform each points of each vector as a rectangle + std::vector squares; + + for(auto p = begin(points); p != end(points); p++) { + sf::RectangleShape *c1 = new sf::RectangleShape(sf::Vector2f(4, 4)); + c1->setPosition(p->x, p->y); + squares.push_back(c1); + } + + // Make the lines + std::vector > lines; + for(auto e = begin(edges); e != end(edges); e++) { + lines.push_back({{ + sf::Vertex(sf::Vector2f((*e).p1.x + 2, (*e).p1.y + 2)), + sf::Vertex(sf::Vector2f((*e).p2.x + 2, (*e).p2.y + 2)) + }}); + } + + while (window.isOpen()) + { + sf::Event event; + while (window.pollEvent(event)) + { + if (event.type == sf::Event::Closed) + window.close(); + } + + window.clear(); + + // Draw the squares + for(auto s = begin(squares); s != end(squares); s++) { + window.draw(**s); + } + + // Draw the lines + for(auto l = begin(lines); l != end(lines); l++) { + window.draw((*l).data(), 2, sf::Lines); + } + + window.display(); + } + + return 0; +} diff --git a/vendor/delaunay/sample.png b/vendor/delaunay/sample.png new file mode 100644 index 0000000..fef352e Binary files /dev/null and b/vendor/delaunay/sample.png differ diff --git a/vendor/delaunay/triangle.h b/vendor/delaunay/triangle.h new file mode 100644 index 0000000..e068a5c --- /dev/null +++ b/vendor/delaunay/triangle.h @@ -0,0 +1,65 @@ +#ifndef H_TRIANGLE +#define H_TRIANGLE + +#include "vector2.h" +#include "edge.h" + +#include +#include + +template +class Triangle +{ + public: + using EdgeType = Edge; + using VertexType = Vector2; + + Triangle(const VertexType &_p1, const VertexType &_p2, const VertexType &_p3) + : p1(_p1), p2(_p2), p3(_p3), + e1(_p1, _p2), e2(_p2, _p3), e3(_p3, _p1) + {} + + bool containsVertex(const VertexType &v) + { + return p1 == v || p2 == v || p3 == v; + } + + bool circumCircleContains(const VertexType &v) + { + float ab = (p1.x * p1.x) + (p1.y * p1.y); + float cd = (p2.x * p2.x) + (p2.y * p2.y); + float ef = (p3.x * p3.x) + (p3.y * p3.y); + + float circum_x = (ab * (p3.y - p2.y) + cd * (p1.y - p3.y) + ef * (p2.y - p1.y)) / (p1.x * (p3.y - p2.y) + p2.x * (p1.y - p3.y) + p3.x * (p2.y - p1.y)) / 2.f; + float circum_y = (ab * (p3.x - p2.x) + cd * (p1.x - p3.x) + ef * (p2.x - p1.x)) / (p1.y * (p3.x - p2.x) + p2.y * (p1.x - p3.x) + p3.y * (p2.x - p1.x)) / 2.f; + float circum_radius = sqrtf(((p1.x - circum_x) * (p1.x - circum_x)) + ((p1.y - circum_y) * (p1.y - circum_y))); + + float dist = sqrtf(((v.x - circum_x) * (v.x - circum_x)) + ((v.y - circum_y) * (v.y - circum_y))); + return dist <= circum_radius; + } + + VertexType p1; + VertexType p2; + VertexType p3; + EdgeType e1; + EdgeType e2; + EdgeType e3; +}; + +template +inline std::ostream &operator << (std::ostream &str, const Triangle & t) +{ + return str << "Triangle:" << std::endl << "\t" << t.p1 << std::endl << "\t" << t.p2 << std::endl << "\t" << t.p3 << std::endl << "\t" << t.e1 << std::endl << "\t" << t.e2 << std::endl << "\t" << t.e3 << std::endl; + +} + +template +inline bool operator == (const Triangle &t1, const Triangle &t2) +{ + return (t1.p1 == t2.p1 || t1.p1 == t2.p2 || t1.p1 == t2.p3) && + (t1.p2 == t2.p1 || t1.p2 == t2.p2 || t1.p2 == t2.p3) && + (t1.p3 == t2.p1 || t1.p3 == t2.p2 || t1.p3 == t2.p3); +} + + +#endif diff --git a/vendor/delaunay/vector2.h b/vendor/delaunay/vector2.h new file mode 100644 index 0000000..d187fe3 --- /dev/null +++ b/vendor/delaunay/vector2.h @@ -0,0 +1,71 @@ +#ifndef H_VECTOR2 +#define H_VECTOR2 + +#include +#include + +template +class Vector2 +{ + public: + // + // Constructors + // + + Vector2() + { + x = 0; + y = 0; + } + + Vector2(T _x, T _y) + { + x = _x; + y = _y; + } + + Vector2(const Vector2 &v) + { + x = v.x; + y = v.y; + } + + void set(const Vector2 &v) + { + x = v.x; + y = v.y; + } + + // + // Operations + // + T dist2(const Vector2 &v) + { + T dx = x - v.x; + T dy = y - v.y; + return dx * dx + dy * dy; + } + + float dist(const Vector2 &v) + { + return sqrtf(dist2(v)); + } + + T x; + T y; + +}; + +template +std::ostream &operator << (std::ostream &str, Vector2 const &point) +{ + return str << "Point x: " << point.x << " y: " << point.y; +} + +template +bool operator == (Vector2 v1, Vector2 v2) +{ + return (v1.x == v2.x) && (v1.y == v2.y); +} + +#endif