diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000000000000000000000000000000000000..332dc4ff12da94374f2284efa0b8ac01a8d0b6b7
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,5 @@
+include README.md
+include requirements.txt
+include images/*
+include examples/*
+exclude */.gitignore
\ No newline at end of file
diff --git a/README.md b/README.md
index dce75ca78b14911f3a9231ed3d4e8ae743e39eb4..54d88c0cabbc58b293291021ddd82fe2f869ac8d 100644
--- a/README.md
+++ b/README.md
@@ -7,9 +7,9 @@
 * [OSMnx](https://osmnx.readthedocs.io/) that includes [NetworkX](https://networkx.org/) and [pandas](https://osmnx.readthedocs.io/)
 * [argparse](https://docs.python.org/3/library/argparse.html)
 
-## Usage
+## Examples
 
-The main script to use is ```get-crossroad-description.py```. You will find a complete description of the possible parameters using the following command:
+The main example to use is ```get-crossroad-description.py```. You will find a complete description of the possible parameters using the following command:
 
 * ```get-crossroad-description.py --help```
 
diff --git a/build/lib/crseg/__init__.py b/build/lib/crseg/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c12f34c653fa81bbe28973516082e785c4fac3e9
--- /dev/null
+++ b/build/lib/crseg/__init__.py
@@ -0,0 +1 @@
+__version__ = '0.1.1'
\ No newline at end of file
diff --git a/src/crseg/crossroad.py b/build/lib/crseg/crossroad.py
similarity index 100%
rename from src/crseg/crossroad.py
rename to build/lib/crseg/crossroad.py
diff --git a/src/crseg/crossroad_connections.py b/build/lib/crseg/crossroad_connections.py
similarity index 100%
rename from src/crseg/crossroad_connections.py
rename to build/lib/crseg/crossroad_connections.py
diff --git a/src/crseg/lane_description.py b/build/lib/crseg/lane_description.py
similarity index 100%
rename from src/crseg/lane_description.py
rename to build/lib/crseg/lane_description.py
diff --git a/src/crseg/link.py b/build/lib/crseg/link.py
similarity index 100%
rename from src/crseg/link.py
rename to build/lib/crseg/link.py
diff --git a/src/crseg/region.py b/build/lib/crseg/region.py
similarity index 100%
rename from src/crseg/region.py
rename to build/lib/crseg/region.py
diff --git a/src/crseg/regionfactory.py b/build/lib/crseg/regionfactory.py
similarity index 100%
rename from src/crseg/regionfactory.py
rename to build/lib/crseg/regionfactory.py
diff --git a/src/crseg/reliability.py b/build/lib/crseg/reliability.py
similarity index 100%
rename from src/crseg/reliability.py
rename to build/lib/crseg/reliability.py
diff --git a/src/crseg/segmentation.py b/build/lib/crseg/segmentation.py
similarity index 100%
rename from src/crseg/segmentation.py
rename to build/lib/crseg/segmentation.py
diff --git a/src/crseg/utils.py b/build/lib/crseg/utils.py
similarity index 100%
rename from src/crseg/utils.py
rename to build/lib/crseg/utils.py
diff --git a/crossroads_segmentation.egg-info/PKG-INFO b/crossroads_segmentation.egg-info/PKG-INFO
new file mode 100644
index 0000000000000000000000000000000000000000..09f1022ab9697445f385969c310a4398354e74cf
--- /dev/null
+++ b/crossroads_segmentation.egg-info/PKG-INFO
@@ -0,0 +1,72 @@
+Metadata-Version: 2.1
+Name: crossroads-segmentation
+Version: 0.1.1
+Summary: Crossroads segmentation is a python tool that produces automatic segmentations of data from OpenStreetMap.
+Home-page: https://gitlab.limos.fr/jmafavre/crossroads-segmentation/
+Author: Jean-Marie Favreau
+Author-email: j-marie.favreau@uca.fr
+License: AGPL-3.0
+Platform: UNKNOWN
+Classifier: License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)
+Classifier: Programming Language :: Python :: 3
+Classifier: Programming Language :: Python :: 3.7
+Description-Content-Type: text/markdown
+
+# crossroads segmentation
+
+**Crossroads segmentation** is a python tool that produces automatic segmentations of data from OpenStreetMap.
+
+## Dependancies
+
+* [OSMnx](https://osmnx.readthedocs.io/) that includes [NetworkX](https://networkx.org/) and [pandas](https://osmnx.readthedocs.io/)
+* [argparse](https://docs.python.org/3/library/argparse.html)
+
+## Examples
+
+The main example to use is ```get-crossroad-description.py```. You will find a complete description of the possible parameters using the following command:
+
+* ```get-crossroad-description.py --help```
+
+This tool is using OSMnx to download OpenStreetMap data from the selected region. It uses a cache, stored in ```cache/``` directory. If a region has already been asked, it will use the cached data and not download it again. You can of course delete the cache directory to download again the data.
+
+The location of the region can be choosen using coordinates (```--by-coordinates LAT LNG```) or using an predefined coordinate defined by a name (```--by-name NAME```). A radius (```-r VALUE```) with a default value of 150 meters can be adjusted to choose the size of the region to consider.
+
+Several outputs are possible:
+
+* to display the segmentation with all the crossings in the region (```--display-segmentation```), or only focussing on the main crossroad (```--display-main-crossroad```) closest to the input coordinate. This second display gives also the branches of the crossroad.
+* to produce a text version of the selection (```--to-text```, ```--to-text-all```) in the standard output
+* to produce a ```json``` file that contains all the detected crossroads (```--to-json-all FILENAME```) or only the main one (```--to-json FILENAME```). Branches are also contained in this output.
+
+
+* 3 parameters (C0, C1 and C2) to drive the creation and merge of the crossroads
+* a maximum number of crossroads in a ring (```--max-cycle-elements NB```), with default value of 10 for the last step of the segmentation.
+
+
+Several of these outputs (```--to-json```, ```--to-json-all```, ```--display-main-crossroad```, ```--to-geopackage```) can be adjusted using the parameter ```--multiscale``` to describe the small crossroad that has been merged to produce the large ones.
+
+## Non regression tests
+
+A very basic non regression test is provided in ```test``` directory. Usage:
+
+* run first ```./regenerate_references.sh```
+* run ```./test.sh``` each time you want to check for regressions
+
+## Visual evaluation
+
+A separated project has been split from this one in order to evaluate segmentation quality. See [crossroads-evaluation](https://gitlab.limos.fr/jmafavre/crossroads-evaluation).
+
+## Examples
+
+
+```./get-crossroad-description.py --by-name POC1  --display-main-crossroad --multiscale```
+
+![POC1](images/POC1.png)
+
+
+```./get-crossroad-description.py --by-name obélisque  --display-segmentation -r 1000```
+
+![R1000](images/R1000.png)
+
+
+
+
diff --git a/crossroads_segmentation.egg-info/SOURCES.txt b/crossroads_segmentation.egg-info/SOURCES.txt
new file mode 100644
index 0000000000000000000000000000000000000000..71e896e9df08f26575afcc93a1fbd68dc6c51878
--- /dev/null
+++ b/crossroads_segmentation.egg-info/SOURCES.txt
@@ -0,0 +1,24 @@
+MANIFEST.in
+README.md
+requirements.txt
+setup.py
+crossroads_segmentation.egg-info/PKG-INFO
+crossroads_segmentation.egg-info/SOURCES.txt
+crossroads_segmentation.egg-info/dependency_links.txt
+crossroads_segmentation.egg-info/top_level.txt
+crseg/__init__.py
+crseg/crossroad.py
+crseg/crossroad_connections.py
+crseg/lane_description.py
+crseg/link.py
+crseg/region.py
+crseg/regionfactory.py
+crseg/reliability.py
+crseg/segmentation.py
+crseg/utils.py
+examples/crossroads-by-name.json
+examples/get-crossroad-description.py
+examples/get-paris-streets.py
+examples/stats-intersection-info.py
+images/POC1.png
+images/R1000.png
\ No newline at end of file
diff --git a/crossroads_segmentation.egg-info/dependency_links.txt b/crossroads_segmentation.egg-info/dependency_links.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8b137891791fe96927ad78e64b0aad7bded08bdc
--- /dev/null
+++ b/crossroads_segmentation.egg-info/dependency_links.txt
@@ -0,0 +1 @@
+
diff --git a/crossroads_segmentation.egg-info/top_level.txt b/crossroads_segmentation.egg-info/top_level.txt
new file mode 100644
index 0000000000000000000000000000000000000000..056069e03c4f6168142565b7fccd6f78978f0f27
--- /dev/null
+++ b/crossroads_segmentation.egg-info/top_level.txt
@@ -0,0 +1 @@
+crseg
diff --git a/src/crseg/.gitignore b/crseg/.gitignore
similarity index 100%
rename from src/crseg/.gitignore
rename to crseg/.gitignore
diff --git a/crseg/__init__.py b/crseg/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c12f34c653fa81bbe28973516082e785c4fac3e9
--- /dev/null
+++ b/crseg/__init__.py
@@ -0,0 +1 @@
+__version__ = '0.1.1'
\ No newline at end of file
diff --git a/crseg/crossroad.py b/crseg/crossroad.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d202a4c2d5473e275d3a35274249f8d3a452988
--- /dev/null
+++ b/crseg/crossroad.py
@@ -0,0 +1,588 @@
+
+import osmnx as ox
+import itertools
+
+
+from . import reliability as rl
+from . import region as r
+from . import utils as u
+from . import lane_description as ld
+
+
+class Crossroad(r.Region):
+    
+
+    def __init__(self, G, node = None, target_id = -1, scale = 2):
+        r.Region.__init__(self, G, target_id)
+
+        # multiplicative coefficient applied to street width 
+        # to obtain distance between the center of the crossroad
+        # and the boundary
+        self.ratio_boundary = scale
+
+        self.max_distance_boundary_polyline = { "motorway": 100, 
+                                                "trunk": 100,
+                                                "primary": 50, 
+                                                "secondary": 40, 
+                                                "tertiary": 30, 
+                                                "unclassified": 25, 
+                                                "residential": 20,
+                                                "living_street": 15,
+                                                "service": 10,
+                                                "default": 10
+                                                }
+
+        if node != None:
+            self.propagate(node)
+            self.build_lanes_description()
+
+    def __str__(self):
+        if hasattr(self, "branches"):
+            return "* id: %s\n* center: %s\n* lanes: %s\n* branches: %s" % (self.id, self.center, len(self.lanes), len(self.branches))
+        else:
+            return "* id: %s\n* center: %s\n* lanes: %s\n" % (self.id, self.center, len(self.lanes))
+
+    def __repr__(self):
+        return "id: %s, center: %s, #nodes: %s" % (self.id, self.center, len(self.nodes))
+
+    def get_center(self):
+        return self.center
+
+    def is_crossroad(self):
+        return True
+
+    def to_text(self):
+        text = "General description:\n" + self.__str__()
+        text += "\nDetails:\n"
+        text += str(self.to_json_data())
+        return text
+
+    def to_json_array(tp, innerNodes, borderNodes, edges, G):
+        crdata = {}
+        crdata["type"] = tp
+        crdata["nodes"] = {}
+        crdata["nodes"]["inner"] = innerNodes
+        crdata["nodes"]["border"] = borderNodes
+        crdata["edges_by_nodes"] = []
+        for e in edges:
+            crdata["edges_by_nodes"].append(e)
+        crdata["coordinates"] = {n: {"x": G.nodes[n]["x"], "y": G.nodes[n]["y"]} for n in innerNodes + borderNodes}
+        return crdata
+
+    def to_json_data(self):
+        data = []
+
+        innerNodes = []
+        borderNodes = []
+        for n in self.nodes:
+            if self.is_boundary_node(n):
+                borderNodes.append(n)
+            else:
+                innerNodes.append(n)            
+
+        data.append(Crossroad.to_json_array("crossroad", innerNodes, borderNodes, self.edges, self.G))
+
+        # for each branch
+        for branch in self.branches:
+            nodes = set()
+            for lane in branch:
+                nodes.add(lane.edge[0])
+                nodes.add(lane.edge[1])
+            edges = [lane.edge for lane in branch]
+            data.append(Crossroad.to_json_array("branch", [], list(nodes), edges, self.G))
+
+
+        return data
+
+    def set_graph_attributes(self, crossroad_attr, branch_attr = None):
+        rid = self.id
+        # set crossroad attribute
+        for n in self.nodes:
+            if len(self.G.nodes[n][crossroad_attr]) == 0:
+                self.G.nodes[n][crossroad_attr] = str(rid)
+            else:
+                self.G.nodes[n][crossroad_attr] += ";" + str(rid)
+        
+        for e in self.edges:
+            if len(self.G[e[0]][e[1]][0][crossroad_attr]) == 0:
+                self.G[e[0]][e[1]][0][crossroad_attr] = str(rid)
+            else:
+                self.G[e[0]][e[1]][0][crossroad_attr] += ";" + str(rid)
+
+        if branch_attr != None:
+            for bid, branch in enumerate(self.branches):
+                cid = str(rid) + "-" + str(bid)
+                for lane in branch:
+                    if len(self.G[lane.edge[0]][lane.edge[1]][0][branch_attr]) == 0:
+                        self.G[lane.edge[0]][lane.edge[1]][0][branch_attr] = cid
+                    else:
+                        self.G[lane.edge[0]][lane.edge[1]][0][branch_attr] += ";" + cid
+
+ 
+
+
+    def get_lane_description_from_edge(self, edge):
+        e = self.G[edge[0]][edge[1]][0]
+        angle = u.Util.bearing(self.G, self.get_center(), edge[0 if self.get_center() == edge[1] else 1])
+        name = e["name"] if "name" in e else None
+        if name == None:
+            # build the path starting from this edge
+            path = u.Util.get_path_to_biffurcation(self.G, edge[0], edge[1])
+
+            # if one of the edges has a name, it's the name of the lane
+            for p1, p2 in zip(path, path[1:]):
+                 e = self.G[p1][p2][0]
+                 if "name" in e:
+                     name = e["name"]
+                     break
+
+            if name == None:
+                # if not found,
+                # consider the last node of this path
+                end = path[-1]
+                # and check if it exists other paths between this end and the crossroad
+                other_paths = []
+                for nb in self.G.neighbors(end):
+                    op = u.Util.get_path_to_biffurcation(self.G, end, nb)
+                    if self.has_node(op[-1]):
+                        # TODO: check if they are parallel
+                        other_paths.append(op)
+                # if only one path exists
+                if len(other_paths) == 1:
+                    o_e = self.G[other_paths[0][0]][other_paths[0][1]][0]
+                    # if yes, the current edge has probably the same name
+                    name = o_e["name"] if "name" in o_e else None
+        return ld.LaneDescription(angle, name, edge)
+
+    def get_lanes_description_from_node(self, border):
+        edges = [(border, nb) for nb in self.G.neighbors(border) if not self.has_edge((nb, border))]
+        return [self.get_lane_description_from_edge(e) for e in edges]
+
+    # estimate the width of the given edge, and deduce the maximum
+    # distance from the center of a crossroad to the boundary of the crossroad
+    def estimate_max_distance_to_boundary(self, edge):
+        w = u.Util.estimate_edge_width(self.G, edge)
+        return w * self.ratio_boundary
+
+    def get_max_lane_width_around_node(self, n):
+        m = 0
+        for nb in self.G.neighbors(n):
+                v = self.estimate_max_distance_to_boundary((n, nb))
+                if v > m:
+                    m = v
+        return m
+
+    def get_max_lane_width(self):
+        m = 0
+        for n in self.nodes:
+            v = self.get_max_lane_width_around_node(n)
+            if v > m:
+                m = v
+        return m
+
+
+    def get_open_paths(self, point, radius):
+        result = []
+
+        for nb in self.G.neighbors(point):
+            if not self.has_edge((nb, point)):    
+                result.append(u.Util.get_path_to_biffurcation(self.G, point, nb, radius))
+
+        return result
+
+    def build_lanes_description(self):
+        self.lanes = []
+
+        center = self.get_center()
+        radius = self.get_max_lane_width() * self.ratio_boundary
+
+        borders = [n for n in self.nodes if self.is_boundary_node(n)]
+
+        for b in borders:
+            if b != center:
+                self.lanes = self.lanes + self.get_lanes_description_from_node(b)
+            else:
+                # go trough all possible paths starting from the center
+                # and add the corresponding lanes
+                open_lanes = self.get_open_paths(center, radius)
+                for ol in open_lanes:
+                    self.lanes.append(self.get_lane_description_from_edge((ol[1], ol[0])))
+        
+
+    def build_crossroads(G, scale):
+        crossroads = []
+        for n in G.nodes:
+            if r.Region.unknown_region_node_in_graph(G, n):
+                if Crossroad.is_reliable_crossroad_node(G, n):
+                    c = Crossroad(G, n, scale = scale)
+
+                    if c.is_straight_crossing():
+                        c.clear_region()
+                    else:
+                        crossroads.append(c)
+        return crossroads
+
+    def is_straight_crossing(self):
+        for n in self.nodes:
+            if len(list(self.G.neighbors(n))) > 2:
+                return False
+        
+        return True
+
+
+    def is_reliable_crossroad_node(G, n):
+        if rl.Reliability.is_weakly_in_crossroad(G, n):
+            return True
+        
+        return False
+
+    def propagate(self, n):
+
+        self.add_node(n)
+        self.center = n
+
+        for nb in self.G.neighbors(n):
+            if self.unknown_region_edge((n, nb)):
+                paths = self.get_possible_paths(n, nb)
+                for path in paths[::-1]:
+                    if path != None and self.is_correct_inner_path(path):
+                        self.add_path(path)
+                        break
+
+
+    def is_correct_inner_node(self, node):
+        return not rl.Reliability.is_weakly_boundary(self.G, node)
+
+            
+
+    def get_closest_possible_biffurcation(self, point):
+        result = -1
+        length = -1
+        for nb in self.G.neighbors(point):
+            path = u.Util.get_path_to_biffurcation(self.G, point, nb)
+            l = u.Util.length(self.G, path)
+            if length < 0 or l < length:
+                length = l
+                result = path[len(path) - 1]
+
+        return result
+
+
+    def is_inner_path_by_osmdata(self, path):
+        for p1, p2 in zip(path, path[1:]):
+            if not "junction" in self.G[p1][p2][0]:
+                return False
+        return True
+
+    def is_correct_inner_path(self, path):
+        if len(path) < 2:
+            return False
+        # loops are not correct inner path in a crossing
+        if path[0] == path[len(path) - 1]:
+            return False
+
+        # use "junction" OSM tag as a good clue
+        if self.is_inner_path_by_osmdata(path):
+            return True
+
+        first = path[0]
+        last = path[len(path) - 1]
+        if rl.Reliability.is_weakly_in_crossroad(self.G, first) and rl.Reliability.is_weakly_boundary(self.G, last):
+            d =  u.Util.length(self.G, path)
+            r = 1
+            if len(list(self.G.neighbors(first))) > 4: # a crossroad with many lanes is larger, thus 
+                r = 2
+            dmax = self.get_max_lane_width_around_node(path[0])
+            if d < dmax * r:
+                return True
+            else:
+                return False
+
+
+    def get_possible_paths(self, n1, n2):
+        results = []
+
+        path = [n1, n2]
+
+        # check first for a boundary
+        while self.is_middle_path_node(path[len(path) - 1]):
+            next = self.get_next_node_along_polyline(path[len(path) - 1], path[len(path) - 2])                
+
+            if next == None:
+                print("ERROR: cannot follow a path")
+                return results
+            path.append(next)
+
+            # if we reach a known region, we stop the expension process
+            if not self.unknown_region_node(next):
+                break
+
+        # if we reach a point with cardinality > 2, we do not consider it
+        if len(list(self.G.neighbors(path[len(path) - 1]))) > 2:
+            return results
+
+        results.append(path)
+
+        # if we reach a strong boundary, we find our path
+        if not self.is_middle_path_node(path[len(path) - 1], True):
+            return results
+
+        path = path.copy()
+
+        # if it's a weak border, we continue until we reach a strong one
+        while self.is_middle_path_node(path[len(path) - 1], True):
+            next = self.get_next_node_along_polyline(path[len(path) - 1], path[len(path) - 2])                
+
+            if next == None:
+                print("ERROR: cannot follow a path")
+                return results
+            path.append(next)
+
+            # if we reach a known region, we stop the expension process
+            if not self.unknown_region_node(next):
+                break
+
+        results.append(path)
+
+        return results
+    
+    def is_middle_path_node(self, node, strong = False):
+        if len(list(self.G.neighbors(node))) != 2:
+            return False
+
+        if strong:
+            return not (rl.Reliability.is_strong_boundary(self.G, node) \
+                    or rl.Reliability.is_strong_in_crossroad(self.G, node))
+        else:
+            return not (rl.Reliability.is_weakly_boundary(self.G, node) \
+                    or rl.Reliability.is_weakly_in_crossroad(self.G, node))
+
+    def get_next_node_along_polyline(self, current, pred):
+        for n in self.G.neighbors(current):
+            if n != pred:
+                return n
+        # cannot append
+        return None
+
+
+    def get_crossroads_in_neighborhood(self, crossroads, scale = 3):
+        result = []
+
+        center = self.get_center()
+        radius = self.get_max_lane_width() * scale
+
+        for c in crossroads:
+            if c.id != self.id and u.Util.distance(self.G, center, c.get_center()) < radius:
+                result.append(c)
+
+        return result
+
+
+    def find_direct_path_to_possible_adjacent_biffurcation(self, point):
+        center = self.get_center()
+        for nb in self.G.neighbors(center):
+            path = u.Util.get_path_to_biffurcation(self.G, center, nb)
+            if path[len(path) - 1] == point:
+                return path
+        return None
+
+    def in_same_cluster(self, crossroad, scale):
+
+        if self.id == crossroad.id:
+            return False  
+
+        angle = u.Util.bearing(self.G, self.get_center(), crossroad.get_center())
+
+
+        # if their is no direct path between centers, they are not
+        # in the same cluster
+        path = self.find_direct_path_to_possible_adjacent_biffurcation(crossroad.get_center())
+        if path == None:
+            return False
+
+        center = self.get_center()
+        d = u.Util.distance(self.G, center, crossroad.get_center())
+        width = self.get_max_lane_width()
+        threshold = width * scale
+        # if it does not exists a strong border between the two crossings, 
+        # reduce the thresold distance by two
+        if not rl.Reliability.has_weakly_boundary_in_path(self.G, path):
+            threshold = threshold / 2
+
+        # if the distance is up to the threshold (considering the possible reduction of the threshold)
+        # the merge cannot be applied
+        if d >= threshold:
+                return False
+        
+        # if the distance is smaller than half the threshold
+        # we merge (without checking for similarity in the name)
+        if d <= threshold / 2:
+            return True
+
+
+        # consider similar branches orthogonal to the junction
+        for b1 in self.lanes:
+            for b2 in crossroad.lanes:
+                if b1.is_similar(b2) and (b1.is_orthogonal(angle) or b2.is_orthogonal(angle)):
+                    return True
+
+
+        return False
+
+    # merge crossroads if they are in a neigborhood defined by scale times the radius of the crossroad
+    # and if they are considered as "in same cluster" (using branch similarities)
+    def get_clusters(crossroads, scale = 2):
+        result = []
+
+        visited = []
+
+        for crossroad in crossroads:
+            if not crossroad.id in visited:
+                visited.append(crossroad.id)
+                cluster = [crossroad]
+            else:
+                cluster = [c for c in result if crossroad in c]
+                if len(cluster) != 1:
+                    cluster = [crossroad]
+                else:
+                    cluster = cluster[0]
+                    result = [c for c in result if not crossroad in c]
+
+            cr_in_neigborhood = crossroad.get_crossroads_in_neighborhood(crossroads, scale)
+            for cr in cr_in_neigborhood:
+                if crossroad.in_same_cluster(cr, scale):
+                    if not cr.id in visited:
+                        visited.append(cr.id)
+                        cluster.append(cr)
+                    else:
+                        # merge clusters
+                        other_cluster = [c for c in result if cr in c]
+                        if len(other_cluster) != 1:
+                            # we check if the merge wasn't processed before
+                            if not cr in cluster:
+                                print("Error while merging two clusters:", crossroad, cr)
+                                print("Other cluster", other_cluster)
+                                print("result", result)
+                        else:
+                            other_cluster = other_cluster[0]
+                            cluster = cluster + other_cluster
+                            result = [c for c in result if not cr in c]
+                        
+            if len(cluster) >= 1:
+                result.append(cluster)
+
+        # finally remove single clusters
+        result = [r for r in result if len(r) > 1]
+        return result
+
+    # add to the current crossing the direct paths that connect the given points
+    def add_direct_paths_between_nodes(self, points):
+        # TODO: avoid too long paths
+        for p1 in points:
+            for n in self.G.neighbors(p1):
+                if not self.has_edge((p1, n)):
+                    path = u.Util.get_path_to_biffurcation(self.G, p1, n)
+                    if path[len(path) - 1] in points and u.Util.length_with_shortcut(self.G, path) < self.diameter():
+                        self.add_path(path)
+
+    # merge all given regions with the current one
+    def merge(self, regions):
+        # add nodes and edges from the other regions
+        for region in regions:
+            for n in region.nodes:
+                self.add_node(n)
+            for e in region.edges:
+                self.add_edge(e)
+
+        old_centers = [r.get_center() for r in regions]
+        old_centers.append(self.get_center())
+
+        # add missing paths between old centers
+        self.add_direct_paths_between_nodes(old_centers)
+
+        # set a new center
+        center = u.Util.centroid(self.G, old_centers)
+        distance = -1
+        new_center = None
+        for n in self.nodes:
+            d = u.Util.distance_to(self.G, n, center)
+            if distance < 0 or d < distance:
+                distance = d
+                new_center = n
+        if new_center != None:
+            self.center = new_center
+        
+        # finally rebuild the branch descriptions
+        self.build_lanes_description()
+
+    # add missing paths (inner paths and paths to boundaries)
+    def add_missing_paths(self, scale = 2, boundaries = True):
+        # add inner paths
+        self.add_direct_paths_between_nodes(self.nodes)
+
+        if boundaries:
+            # add paths to missing boundaries within the given scale
+            max_length = scale * self.get_max_lane_width()
+            for p1 in self.nodes:
+                if self.G.nodes[p1][rl.Reliability.boundary_reliability] <= rl.Reliability.uncertain:
+                    for n in self.G.neighbors(p1):
+                        if not self.has_edge((p1, n)) and self.G[p1][n][0][r.Region.label_region] == -1:
+                            path = rl.Reliability.get_path_to_boundary(self.G, p1, n)
+                            while len(path) > 2 and self.G[path[-2]][path[-1]][0][r.Region.label_region] != -1:
+                                path.pop()
+                            # find a boundary node inside the path and cut it
+                            if len(path) > 0 and u.Util.length(self.G, path) < max_length:
+                                self.add_path(path)
+
+        # finally rebuild the branch descriptions
+        self.build_lanes_description()
+
+    def compute_branches(self):
+
+        self.branches = []
+        
+        # for each lane
+        for lane in self.lanes:
+            mbranches = []
+
+            # check if it's similar to a lane already in a built branch
+            for i, branch in enumerate(self.branches):
+                nb = len([l for l in branch if l.is_similar(lane)])
+                if nb > 0:
+                    mbranches.append(i)
+            # if not, create a new branch
+            if len(mbranches) == 0:
+                self.branches.append([lane])
+            else:
+                # merge the similar branches
+                self.branches[mbranches[0]].append(lane)
+                for idb in mbranches[1:]:
+                    self.branches[mbranches[0]] += self.branches[idb]
+                    self.branches[idb] = []
+            
+            # remove empty branches
+            self.branches = [ b for b in self.branches if len(b) > 0]
+
+    # if the given edge is not part of a branch, it returns -1
+    # otherwise, it returns the index of the branch (in this crossing) that contains the given edge
+    def get_branch_id(self, e):
+
+        for i, branch in enumerate(self.branches):
+            for lane in branch:
+                if lane.equals(e):
+                    return i
+
+        return -1
+
+    # get the maximum estimated width of a branch
+    def max_branch_width(self):
+        if not hasattr(self, "branches"):
+            self.compute_branches()
+        return max([self.estimate_branch_width(b) for b in self.branches])
+
+
+    def estimate_branch_width(self, branch):
+        return sum([self.estimate_lane_width(l) for l in branch])
+            
+    def estimate_lane_width(self, lane):
+        return u.Util.estimate_edge_width(self.G, lane.edge)
diff --git a/crseg/crossroad_connections.py b/crseg/crossroad_connections.py
new file mode 100644
index 0000000000000000000000000000000000000000..e64f92d2f46c89077df22234bb22c4810fe51ac6
--- /dev/null
+++ b/crseg/crossroad_connections.py
@@ -0,0 +1,235 @@
+from . import reliability as rel
+from . import utils as u
+import math
+
+class CrossroadConnections:
+
+    typology_crossroad = 1
+    typology_link = 2
+    typology_unknown = 0
+
+    ratio_single_path = 5
+
+    max_distance_connection = 50
+    max_loop_distance = max_distance_connection * math.pi
+
+    # the connection_threshold corresponds to a coefficient used for connection
+    # between crossroads, multiplied by an estimation of the size of a crossroad
+    # defined by an estimation of the branch width
+    def __init__(self, regions, connection_threshold = 4):
+        self.regions = regions
+        self.connection_threshold = connection_threshold
+
+        self.init_structure()
+
+    def init_structure(self):
+
+        # build the list of crossroad and links
+        self.crossroads = []
+        self.links = []
+        self.crossroads_max_branch_width = {}
+
+
+
+        for rid in self.regions:
+            if not hasattr(self, "G"):
+                self.G = self.regions[rid].G
+            if self.regions[rid].is_crossroad():
+                self.crossroads.append(rid)
+                self.crossroads_max_branch_width[rid] = self.regions[rid].max_branch_width()
+            elif self.regions[rid].is_link():
+                self.links.append(rid)
+
+        # for each node, build the list of intersecting regions
+        self.regionsByNode = {}
+        for rid in self.regions:
+            for n in self.regions[rid].nodes:
+                self.add_node_region(rid, n)
+
+        # for each region, the list of its adjacent regions (only links for crossroads, and only crossroads for links)
+        self.adjacencies = {}
+        # for each junction node
+        for n in self.regionsByNode:
+            # for reach region associated to this node
+            for r1 in self.regionsByNode[n]:
+                # add an adjacency between this region and all regions connected via the current node
+                self.add_adjacencies(r1, n, self.regionsByNode[n])
+
+        # compute a list of connections between crossroads
+        self.compute_connected_crossroads()
+
+    def compute_connected_crossroads(self):
+        self.compute_initial_connections()
+
+        merged = {}
+        
+
+        # merge multiple instances of the same pair of connected crossroads
+        for c in self.connected_crossroads:
+            if (c[0], c[1]) in merged:
+                merged[(c[0], c[1])].append(c[2])
+            else:
+                merged[(c[0], c[1])] = [c[2]]
+        
+        new_list = []
+        for c in merged:
+            new_list.append((c[0], c[1], merged[c]))
+        
+        self.connected_crossroads = new_list
+
+    def get_max_distance_connection(self, cr, cr2):
+        result = max([self.crossroads_max_branch_width[cr], self.crossroads_max_branch_width[cr2]]) * self.connection_threshold
+        if result > self.max_distance_connection:
+            result = self.max_distance_connection
+        return result
+
+    def get_max_loop_distance(self, cr):
+        result = self.crossroads_max_branch_width[cr] * self.connection_threshold * math.pi
+        if result > self.max_loop_distance:
+            result = self.max_loop_distance
+        return result
+
+    def compute_initial_connections(self):
+        self.connected_crossroads = []
+        # for each crossroad region
+        for cr in self.crossroads:
+            # for each adjacent links
+            for l in self.adjacencies[cr]:
+                # then find the reachable crossings from this link
+                for cr2 in self.adjacencies[l]:
+                    # only considering the ones with an ID higher to the ID of the initial crossroad region
+                    if self.regions[cr].id < self.regions[cr2].id:
+                        path, distance = self.get_path_in_link(l, cr, cr2)
+
+                        if path != None:
+                            # add them as a pair with the corresponding path only if the path is not too long
+                            maxD = self.get_max_distance_connection(cr, cr2)
+                            if distance < maxD:
+                                distanceCC = distance + u.Util.distance(self.G, self.regions[cr].get_center(), path[0]) + u.Util.distance(self.G, self.regions[cr2].get_center(), path[-1])
+                                close = distanceCC < maxD / self.ratio_single_path
+                                self.connected_crossroads.append((cr, cr2, (path, l, close)))
+
+    # return a path (defined by a list of nodes) contained in the given link l that connects
+    # the two given crossroad regions (cr1 and cr2)
+    def get_path_in_link(self, l, cr1, cr2):
+        cr1n = [n for n in self.regions[l].nodes if cr1 in self.regionsByNode[n]]
+        cr2n = [n for n in self.regions[l].nodes if cr2 in self.regionsByNode[n]]
+
+
+        path = self.regions[l].get_path(cr1n, cr2n, u.Util.distance_with_shortcut)
+
+        if path == None:
+            return (None, None)
+        # if we identify a node which was classified as a possible crossroad, the
+        # probability that this path is probably an inner path of a crossroad increases.
+        # We thus reduce the path length
+        lInside = [self.regions[l].G.nodes[p][rel.Reliability.crossroad_reliability] for p in path[0][1:-1]]
+        nbPossible = len([r for r in lInside if r <= rel.Reliability.strongly_yes and r > rel.Reliability.strongly_no])
+        if nbPossible > 0:
+            path = (path[0], path[1] / math.log(math.e * (nbPossible + 1)))
+
+        return path
+
+    def add_adjacencies(self, r1, node, regions):
+        if not r1 in self.adjacencies:
+            self.adjacencies[r1] = {}
+        for r2 in regions:
+            if r2 != r1 and self.get_typology(r1) != self.get_typology(r2):
+                if not r2 in self.adjacencies[r1]:
+                    self.adjacencies[r1][r2] = []
+                self.adjacencies[r1][r2].append(node)
+
+    def get_typology(self, rid):
+        if rid in self.crossroads:
+            return CrossroadConnections.typology_crossroad
+        elif rid in self.links:
+            return CrossroadConnections.typology_link
+        else:
+            return CrossroadConnections.typology_unknown
+            
+
+    def add_node_region(self, rid, nid):
+        if not nid in self.regionsByNode:
+            self.regionsByNode[nid] = []
+        self.regionsByNode[nid].append(rid)
+
+    def get_pairs(self):
+        return [connected for connected in self.connected_crossroads if len(connected[2]) >= 2 or (len(connected[2]) == 1 and connected[2][0][2])]
+
+    def get_cycles(self, max_length = 5):
+        results = []
+
+        for c in self.crossroads:
+            results += self.get_cycles_from_crossroad(c, max_length)
+
+        results = self.get_unique_cycles(results)
+
+        return results
+
+    def get_unique_cycles(self, cycles):
+        result = []
+        seen = []
+
+        for c in cycles:
+            celems = set([e[0] for e in c])
+            if not celems in seen:
+                result.append(c)
+                seen.append(celems)
+
+        return result
+
+    def get_connected_crossroads(self, cr):
+        result = []
+        for connected in self.connected_crossroads:
+            if connected[0] == cr:
+                result.append((connected[1], connected[2]))
+            elif connected[1] == cr:
+                result.append((connected[0], connected[2]))
+        return result
+
+    def get_cycle_length(self, cycle, direct):
+        result = 0
+
+        for p in cycle:
+            if len(p[1]) > 0:
+                if direct:
+                    result += u.Util.distance(self.G, p[1][0][0][0], p[1][0][0][-1])
+                else:
+                    result += u.Util.length(self.G, p[1][0][0])
+
+        return result
+
+    def get_cycles_from_crossroad(self, cr, max_length):
+        paths = [ [(cr, [])] ]
+        results = []
+
+        max_perimeter = self.get_max_loop_distance(cr)
+        # increase step by step the possible paths
+        for l in range(0, max_length):
+            new_paths = []
+            # for each existing path, compute all possible extensions (without backward)
+            for p in paths:
+                # check all possible next steps
+                for next in self.get_connected_crossroads(p[-1][0]):
+                    if len(p) == 1 or p[-2][0] != next[0] and not self.intersects_path_link(next[1], p):
+                        if next[0] == p[0][0]:
+                            # loop detection
+                            new = p + [next]
+                            if self.get_cycle_length(new, True) < max_perimeter:
+                                results.append(new)
+                        else:
+                            # ongoing loop
+                            new_paths.append(p + [next])
+            paths = new_paths
+
+        return results
+
+    def intersects_path_link(self, nextpathlinks, path):
+        n_links = set([n for c in nextpathlinks for n in c[0]])
+
+        for p in path:
+            for l in p[1]:
+                if len(set.intersection(set(l[0]), n_links)) > 1:
+                    return True
+        return False
+
diff --git a/crseg/lane_description.py b/crseg/lane_description.py
new file mode 100644
index 0000000000000000000000000000000000000000..20486abab13e93a922f8aa9c95577a2456bee5e6
--- /dev/null
+++ b/crseg/lane_description.py
@@ -0,0 +1,36 @@
+
+from . import utils as u
+
+class LaneDescription:
+
+    def __init__(self, angle, name, edge):
+        self.angle = angle
+        self.name = name
+        self.edge = edge
+
+    def is_similar(self, bd, angle_similarity = 90):
+        if self.name == None or bd.name == None:
+            return False
+        if self.name != bd.name:
+            return False
+
+        if u.Util.angular_distance(self.angle, bd.angle) < angle_similarity:
+            return True
+
+        return False
+
+    def is_orthogonal(self, angle, epsilon = 45):
+        diff = u.Util.angular_distance(self.angle, angle)
+        if diff >= 90 - epsilon and diff <= 90 + epsilon:
+            return True
+        return False
+
+    def equals(self, edge):
+        return (self.edge[0] == edge[0] and self.edge[1] == edge[1]) or \
+                (self.edge[1] == edge[0] and self.edge[0] == edge[1])
+
+    def __str__(self):
+        return "%s : %s" % (self.name, self.angle)
+
+    def __repr__(self):
+        return self.__str__()
\ No newline at end of file
diff --git a/crseg/link.py b/crseg/link.py
new file mode 100644
index 0000000000000000000000000000000000000000..9857c6c2b2dbd98abc5a16de6878de17c39b96b2
--- /dev/null
+++ b/crseg/link.py
@@ -0,0 +1,39 @@
+
+import osmnx as ox
+
+from . import reliability as rl
+from . import region as r
+from . import utils as u
+from . import lane_description as ld
+
+
+class Link(r.Region):
+    
+    def __init__(self, G, node1 = None, node2 = None, target_id = -1):
+        r.Region.__init__(self, G, target_id)
+        if node1 != None:
+            self.add_node(node1)
+            self.filled = False
+            if node2 != None:
+                if G.nodes[node2][r.Region.label_region] != -1:
+                    self.filled = True
+                self.add_node(node2)
+                self.add_edge((node1, node2))
+
+    def is_link(self):
+        return True
+
+    def propagate(self):
+        # only propagate if this link is not filled
+        if not self.filled and len(self.nodes) > 0:
+            start = self.nodes[-1]
+            self.propagate_from_node(start)
+    
+    def propagate_from_node(self, start):
+        for nb in self.G.neighbors(start):
+            if self.G[start][nb][0][r.Region.label_region] == -1:
+                open = self.G.nodes[nb][r.Region.label_region] == -1
+                self.add_node(nb)
+                self.add_edge((start, nb))
+                if open:
+                    self.propagate_from_node(nb)
diff --git a/crseg/region.py b/crseg/region.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ee7eb73e01bdd20f3ee38bd43b96e8a52542f49
--- /dev/null
+++ b/crseg/region.py
@@ -0,0 +1,151 @@
+import networkx as nx
+import osmnx as ox
+import pandas as pd
+import random
+import math
+
+
+from . import utils as u
+from . import reliability as r
+
+
+
+class Region:
+
+    id_region = 0
+
+    label_region = "region"
+    regiontag_prefix = "cr.region-"
+
+    def __init__(self, G, target_id = -1):
+        self.G = G
+        self.edges = []
+        self.nodes = []
+        if target_id == -1:
+            self.id = Region.id_region
+            Region.id_region += 1
+        else:
+            self.id = target_id
+            if Region.id_region <= target_id:
+                Region.id_region = target_id + 1
+
+    def is_crossroad(self):
+        return False
+
+    def is_link(self):
+        return False
+
+    def clear_region(self):
+        # remove edges
+        for e in self.edges:
+            self.G[e[0]][e[1]][0][Region.label_region] = -1
+        
+        # then remove nodes
+        for n in self.nodes:
+            self.G.nodes[n][Region.label_region] = -1
+    
+    # return true if all the nodes of the given region are part of the current region
+    def contains(self, region):
+        for n in region.nodes:
+            if not n in self.nodes:
+                return False
+        return True
+
+    def init_attr(G):
+        nx.set_edge_attributes(G, values=-1, name=Region.label_region)
+        nx.set_node_attributes(G, values=-1, name=Region.label_region)
+
+    def unknown_region_node_in_graph(G, n):
+        return G.nodes[n][Region.label_region] == -1
+
+    def unknown_region_edge_in_graph(G, e):
+        return G[e[0]][e[1]][0][Region.label_region] == -1
+
+    def unknown_region_node(self, n):
+        return Region.unknown_region_node_in_graph(self.G, n)
+
+    def unknown_region_edge(self, e):
+        return Region.unknown_region_edge_in_graph(self.G, e)
+
+    def clear_node_region_in_grah(G, n):
+        G.nodes[n][Region.label_region] = -1
+
+    def add_path(self, path):
+        for n in path:
+            self.add_node(n)
+        for n1, n2 in zip(path, path[1:]):
+            self.add_edge((n1, n2))
+
+    def add_paths(self, paths):
+        for path in paths:
+            self.add_path(path)
+
+    def add_node(self, n):
+        if n not in self.nodes:
+            self.nodes.append(n)
+        self.G.nodes[n][Region.label_region] = self.id
+
+    def add_edge(self, e):
+        if not self.has_edge(e):
+            self.edges.append(e)
+        self.G[e[0]][e[1]][0][Region.label_region] = self.id
+
+    def add_path(self, path):
+        for p in path:
+            self.add_node(p)
+        for p1, p2 in zip(path, path[1:]):
+            self.add_edge((p1, p2))
+
+    def has_edge(self, e):
+        return (e[0], e[1]) in self.edges or (e[1], e[0]) in self.edges
+
+    def has_node(self, n):
+        return n in self.nodes
+
+    def edges_with_node(self, n):
+        return [ e for e in self.edges if e[0] == n or e[1] == n]
+
+    def is_boundary_node(self, n):
+        nbnb = len(list(self.G.neighbors(n)))
+        nbEdgesInside = len([e for e in self.edges if e[0] == n or e[1] == n])
+        return nbnb != nbEdgesInside
+
+    def centroid(self):
+        return u.Util.centroid(self.G, self.nodes)
+
+    def boundary_nodes(self):
+        result = []
+        for n in self.nodes:
+            if self.is_boundary_node(n):
+                result.append(n)
+        return result
+
+    def diameter(self):
+        # TODO: not optimized
+        result = 0
+        for n1 in self.nodes:
+            for n2 in self.nodes:
+                d = u.Util.distance(self.G, n1, n2)
+                if d > result:
+                    result = d
+        return result
+
+    # return a shortest path inside the current region that connects a node from nodes1 and a node from nodes2
+    # if no such path exists, it returns an empty path
+    def get_path(self, nodes1, nodes2, weight_function = None):
+        if len(nodes1) == 0 or len(nodes2) == 0:
+            return []
+
+        cutoff = 3 * self.diameter() # large number in case of non straight paths
+        # get all possible paths in the current region from the input nodes
+        distances, paths = nx.multi_source_dijkstra(self.G, nodes1, 
+            weight = lambda n1, n2, d: (weight_function(self.G, n1, n2) if weight_function != None else u.Util.distance(self.G, n1, n2)) if self.has_edge((n1, n2)) else None, 
+            cutoff = cutoff)
+
+        # keep paths that reach one of the given nodes
+        distances = {k: v for k, v in distances.items() if k in nodes2}
+        if len(distances) == 0:
+            return None
+        # keep the best one
+        best_target = min(distances, key=distances.get)
+        return paths[best_target], distances[best_target]
diff --git a/crseg/regionfactory.py b/crseg/regionfactory.py
new file mode 100644
index 0000000000000000000000000000000000000000..0c26b62d26b5feaaba97ecac516f58d94567d2d3
--- /dev/null
+++ b/crseg/regionfactory.py
@@ -0,0 +1,70 @@
+from . import crossroad as cr
+from . import region as rg
+from . import utils as u
+from . import link as lk
+
+
+class RegionFactory:
+
+    def clone(region):
+        if region.is_crossroad():
+            result = cr.Crossroad(region.G)
+        elif region.is_link():
+            result = lk.Link(region.G)
+            result.filled = region.filled
+        else:
+            result = rg.Region(region.G)
+        result.nodes = region.nodes.copy()
+        result.edges = region.edges.copy()
+
+        result.lanes = region.lanes.copy()
+        result.center = region.center
+
+        return result
+
+    def rebuild_regions_from_tags(G):
+        # rebuild regions from metadata
+        # ie for each region label associated to an edge or a node, create the corresponding region if it does not exists
+        # considering the metadata to create a branch or a crossroad
+
+        regions = {}
+        for n in G.nodes:
+            if G.nodes[n][rg.Region.label_region] != -1:
+                id = G.nodes[n][rg.Region.label_region]
+                if not id in regions:
+                    regions[id] = cr.Crossroad(G, target_id = int(id)) if G.graph[rg.Region.regiontag_prefix + str(id)] == "crossroad" else rg.Region(G, target_id = id)
+                regions[id].add_node(n)
+
+        for e in G.edges:
+            if G[e[0]][e[1]][0][rg.Region.label_region] != -1:
+                id = G[e[0]][e[1]][0][rg.Region.label_region]
+                if not id in regions:
+                    regions[id] = cr.Crossroad(G, target_id = id) if G.graph[rg.Region.regiontag_prefix + str(id)] == "crossroad" else rg.Region(G, target_id = id)
+                regions[id].add_edge(e)
+                regions[id].add_node(e[0])
+                regions[id].add_node(e[1])
+
+        return regions
+
+    def build_links_between_crossings(G, crossings):
+        links = {}
+
+        for rid in crossings:
+            for b in crossings[rid].boundary_nodes():
+                # if some edges are not in a region                    
+                if u.Util.has_non_labeled_adjacent_edge(G, b):
+                    # for each edge outside of a region, create a link region
+                    for nb in G.neighbors(b):
+                        if G[b][nb][0][rg.Region.label_region] == -1:
+                            l = lk.Link(G, b, nb)
+                            l.propagate()
+                            links[l.id] = l
+                else:
+                    # create a link region with a single node (if not yet created)
+                    exists = len([idl for idl in links if links[idl].has_node(b)])
+                    if exists == 0:
+                        l = lk.Link(G, b)
+                        links[l.id] = l
+                    
+
+        return links
\ No newline at end of file
diff --git a/crseg/reliability.py b/crseg/reliability.py
new file mode 100644
index 0000000000000000000000000000000000000000..bed27424dc6ddb323cb0d7f19e934a8f7593e0a7
--- /dev/null
+++ b/crseg/reliability.py
@@ -0,0 +1,167 @@
+
+
+
+import networkx as nx
+import osmnx as ox
+import pandas as pd
+import random
+import math
+
+
+from . import region as r
+from . import utils as u
+
+
+
+class Reliability:
+
+    # distance for a path to be considered as a branch
+    distance_inner_branch = 50
+ 
+
+    boundary_reliability = "reliability boundary"
+    crossroad_reliability = "reliability.crossroad"
+
+
+    moderate_boundary = [ "stop", "traffic_signals", "motorway_junction", "give_way" ]
+    possible_boundary = [ "crossing"]
+    strongly_no_boundary_attr = [ "bus_stop", "milestone", "steps", "elevator" ]
+
+    strongly_yes = 1000.0
+    strongly_no = 0.0
+
+    uncertain = (strongly_yes + strongly_no) / 2
+
+    weakly_yes = (strongly_yes + uncertain) / 2
+    weakly_no = (strongly_no + uncertain) / 2
+
+    moderate_yes = (weakly_yes + strongly_yes) / 2
+    moderate_no = (weakly_no + strongly_no) / 2
+
+
+    def init_attr(G):
+        nx.set_node_attributes(G, values=Reliability.uncertain, name=Reliability.boundary_reliability)
+        nx.set_node_attributes(G, values=Reliability.uncertain, name=Reliability.crossroad_reliability)
+        Reliability.compute_nodes_reliability(G)
+
+        nx.set_edge_attributes(G, values=Reliability.uncertain, name=Reliability.crossroad_reliability)
+        Reliability.compute_edges_reliability(G)
+
+    def compute_edges_reliability(G):
+
+        for e in G.edges():
+            length = u.Util.distance(G, e[0], e[1])
+            if "junction" in G[e[0]][e[1]][0]:
+                G[e[0]][e[1]][0][Reliability.crossroad_reliability] = Reliability.strongly_yes
+
+    def compute_nodes_reliability(G):
+
+        for n in G.nodes:
+            nb_neighbors = len(list(G.neighbors(n)))
+
+            if "highway" in G.nodes[n]:
+                if nb_neighbors == 2:
+                    G.nodes[n][Reliability.crossroad_reliability] = Reliability.strongly_no
+                
+                
+
+
+                if G.nodes[n]["highway"] in Reliability.strongly_no_boundary_attr:
+                    G.nodes[n][Reliability.boundary_reliability] = Reliability.moderate_no
+                elif G.nodes[n]["highway"] in Reliability.possible_boundary and nb_neighbors <= 3:
+                    G.nodes[n][Reliability.boundary_reliability] = Reliability.strongly_yes
+                elif G.nodes[n]["highway"] in Reliability.moderate_boundary and nb_neighbors <= 3:
+                    G.nodes[n][Reliability.boundary_reliability] = Reliability.moderate_yes
+                    G.nodes[n][Reliability.crossroad_reliability] = Reliability.moderate_yes
+                
+                if nb_neighbors >= 3:
+                    G.nodes[n][Reliability.crossroad_reliability] = Reliability.strongly_yes
+            else:
+                if nb_neighbors == 2:
+                    G.nodes[n][Reliability.boundary_reliability] = Reliability.strongly_no
+                    G.nodes[n][Reliability.crossroad_reliability] = Reliability.strongly_no
+                elif nb_neighbors >= 4:
+                        G.nodes[n][Reliability.crossroad_reliability] = Reliability.strongly_yes
+                elif nb_neighbors == 3:
+                        adj_streetnames = u.Util.get_adjacent_streetnames(G, n)
+
+                        if len(adj_streetnames) > 1:
+                            # more than one street name, it is probably part of a crossroad
+                            G.nodes[n][Reliability.crossroad_reliability] = Reliability.moderate_yes
+                        else:
+                            # only one name
+                            if u.Util.is_part_of_local_triangle(G, n) or u.Util.is_street_separation(G, n):
+                                G.nodes[n][Reliability.crossroad_reliability] = Reliability.moderate_no
+                            else:
+                                G.nodes[n][Reliability.crossroad_reliability] = Reliability.moderate_yes
+
+            if nb_neighbors > 2:
+                # if all adjacent edges are service=parking_aisle, then it is not an intersection
+                if u.Util.is_inside_parking(G, n):
+                    G.nodes[n][Reliability.crossroad_reliability] = Reliability.strongly_no
+
+
+    def get_best_reliability_node(G, n):
+        
+        if G.nodes[n][Reliability.crossroad_reliability] > G.nodes[n][Reliability.boundary_reliability]:
+            return Reliability.crossroad_reliability
+        else:
+            return Reliability.boundary_reliability
+
+    def has_strong_boundary_in_path(G, path):
+        for p in path:
+            if Reliability.is_strong_boundary(G, p):
+                return True
+        return False
+
+    def has_weakly_boundary_in_path(G, path):
+        for p in path:
+            if Reliability.is_weakly_boundary(G, p):
+                return True
+        return False
+
+    def is_strong_boundary(G, n):
+        return G.nodes[n][Reliability.boundary_reliability] == Reliability.strongly_yes
+
+    def is_weakly_boundary(G, n):
+        return G.nodes[n][Reliability.boundary_reliability] >= Reliability.weakly_yes
+
+    def is_weakly_no_boundary(G, n):
+        return G.nodes[n][Reliability.boundary_reliability] <= Reliability.weakly_no
+
+    def is_strong_no_boundary(G, n):
+        return G.nodes[n][Reliability.boundary_reliability] == Reliability.strongly_no
+
+    def is_strong_in_crossroad(G, n):
+        return G.nodes[n][Reliability.crossroad_reliability] == Reliability.strongly_yes
+
+    def is_weakly_in_crossroad(G, n):
+        return G.nodes[n][Reliability.crossroad_reliability] >= Reliability.weakly_yes
+
+    def is_weakly_not_in_crossroad(G, n):
+        return G.nodes[n][Reliability.crossroad_reliability] <= Reliability.weakly_no
+
+    def is_strong_not_in_crossroad(G, n):
+        return G.nodes[n][Reliability.crossroad_reliability] == Reliability.strongly_no
+
+    def is_weakly_in_crossroad_edge(G, e):
+        return G[e[0]][e[1]][0][Reliability.crossroad_reliability] >= Reliability.weakly_yes
+
+    def is_strong_in_crossroad_edge(G, e):
+        return G[e[0]][e[1]][0][Reliability.crossroad_reliability] == Reliability.strongly_yes
+
+
+    def get_path_to_boundary(G, n1, n2, max = -1):
+        path = [n1, n2]
+        length = u.Util.distance(G, n1, n2)
+
+
+        while (max < 0 or length < max) and u.Util.is_middle_polyline(G, path[len(path) - 1]):
+            last = path[len(path) - 1]
+            if Reliability.is_weakly_boundary(G, last):
+                return path
+            path.append(u.Util.get_opposite_node(G, path[len(path) - 1], path[len(path) - 2]))
+            length += u.Util.distance(G, path[len(path) - 2], path[len(path) - 1])
+        # we reach a split node without reaching a boundary node
+        return []
+        
diff --git a/crseg/segmentation.py b/crseg/segmentation.py
new file mode 100644
index 0000000000000000000000000000000000000000..22132921b995ce074e53da5e7b7921a28ca090d2
--- /dev/null
+++ b/crseg/segmentation.py
@@ -0,0 +1,513 @@
+
+import networkx as nx
+import osmnx as ox
+import pandas as pd
+import random
+import math
+import json
+
+
+from . import crossroad as cr
+from . import region as rg
+from . import regionfactory as rf
+from . import reliability as rel
+from . import crossroad_connections as cc
+
+class Segmentation:
+
+    def __init__(self, G, init=True, C0 = 2, C1 = 2.5, C2 = 4, max_cycle_elements = 10):
+        self.G = G
+        self.regions = {}
+        self.C0 = C0
+        self.C1 = C1
+        self.C2 = C2
+        self.max_cycle_elements = max_cycle_elements
+        random.seed()
+        if init:
+            rel.Reliability.init_attr(self.G)
+        else:
+            self.regions = rf.RegionFactory.rebuild_regions_from_tags(self.G)
+
+
+    def set_tags_only_regions(self):
+        # clear tags
+        for n in self.G.nodes:
+            self.G.nodes[n][rg.Region.label_region] = -1
+        for u, v, a in self.G.edges(data = True):
+            self.G[u][v][0][rg.Region.label_region] = -1
+
+        # set tags wrt crossroad regions
+        for rid in self.regions:
+            region = self.regions[rid]
+            if region.is_crossroad():
+                for n in region.nodes:
+                    self.G.nodes[n][rg.Region.label_region] = rid
+                for e in region.edges:
+                    self.G[e[0]][e[1]][0][rg.Region.label_region] = rid
+                
+
+
+    def process(self):
+
+        # init flags
+        rg.Region.init_attr(self.G)
+
+        self.regions = {}
+
+        # first build crossroads
+        crossroads = cr.Crossroad.build_crossroads(self.G, self.C0)
+        for c in crossroads:
+            self.regions[c.id] = c
+
+        # group subparts of crossroads together if they are part of the same crossing (using street names)
+        clusters = cr.Crossroad.get_clusters(crossroads, self.C1)
+
+        # for each cluster
+        for cluster in clusters:
+            # merge them
+            if len(cluster) > 1:
+                cluster[0].merge(cluster[1:])
+            for o in cluster[1:]:
+                del self.regions[o.id]
+
+        # maximum length for missing maths
+        scale_missing = self.C0
+        # add inner paths and missing boundaries
+        self.add_missing_paths(scale = scale_missing)
+        
+        # build links between regions
+        links = rf.RegionFactory.build_links_between_crossings(self.G, self.regions)
+        self.regions.update(links)
+        self.set_tags_only_regions()
+
+        # merge crossings
+        self.merge_linked_crossroads()
+
+        # add inner paths and missing boundaries (again)
+        self.add_missing_paths(scale = scale_missing, boundaries = False)
+
+        # create branch regions
+        for rid in self.regions:
+            if self.regions[rid].is_crossroad():
+                self.regions[rid].compute_branches()
+        
+        for rid in self.inner_regions:
+            self.inner_regions[rid].compute_branches()
+            
+
+    def merge_linked_crossroads(self):
+        self.inner_regions = {}
+        newIDs = {}
+
+        cconnections = cc.CrossroadConnections(self.regions, self.C2)
+
+        # merge multi crossings (triangles, rings, etc)
+        for cycle in cconnections.get_cycles(self.max_cycle_elements):
+            cWithIDs = [cr if cr[0] in self.regions else (newIDs[cr[0]], cr[1]) for cr in cycle][:-1]
+            ids = [x[0] for x in cWithIDs]
+            
+            if len(set(ids)) > 1:
+                firstID = ids[0]
+                # add all regions as inner regions (of a bigger one)
+                for id in ids:
+                    self.add_inner_region(self.regions[id])
+
+                for cr1, cr2 in zip(cWithIDs, cWithIDs[1:]):
+                    id2 = newIDs[cr2[0]] if cr2[0] in newIDs else cr2[0]
+
+                    # add paths that connects cr1 and cr2
+                    self.regions[firstID].add_paths([x[0] for x in cr2[1]])
+                    if id2 != firstID:
+                        self.regions[firstID].merge([self.regions[id2]])
+                        del self.regions[id2]
+                        newIDs[id2] = firstID
+                        for nid in newIDs:
+                            if newIDs[nid] == id2:
+                                newIDs[nid] = firstID
+
+        # merge bi-connected crossings
+        for pairs in cconnections.get_pairs():
+            id1 = pairs[0] if pairs[0] in self.regions else newIDs[pairs[0]]
+            id2 = pairs[1] if pairs[1] in self.regions else newIDs[pairs[1]]
+            if id1 != id2:
+                # add the two regions to the inner regions (of a bigger one)
+                self.add_inner_region(self.regions[id1])
+                self.add_inner_region(self.regions[id2])
+                # add paths that are connecting these two regions
+                self.regions[id1].add_paths([x[0] for x in pairs[2]])
+                # merge the two regions
+                self.regions[id1].merge([self.regions[id2]])
+                # remove the old one
+                del self.regions[id2]
+                # update IDs
+                newIDs[id2] = id1
+                for nid in newIDs:
+                    if newIDs[nid] == id2:
+                        newIDs[nid] = id1
+
+
+    def add_inner_region(self, region):
+        # clone the given region and add it to the inner_regions structure
+        newRegion = rf.RegionFactory.clone(region)
+        self.inner_regions[newRegion.id] = newRegion
+
+    def add_missing_paths(self, boundaries = True, scale = 2):
+        for rid in self.regions:
+            region = self.regions[rid]
+            if region.is_crossroad():
+                region.add_missing_paths(boundaries = boundaries, scale = scale)
+
+
+    def in_crossroad_region(self, e):
+        tag = self.G[e[0]][e[1]][0][rg.Region.label_region]
+        if tag == -1:
+            False
+        else: 
+            return self.regions[tag].is_crossroad()
+
+
+    def get_adjacent_crossroad_regions(self, n):
+        result = []
+        for nb in self.G.neighbors(n):
+            e = (n, nb)
+            tag = self.G[e[0]][e[1]][0][rg.Region.label_region]
+            if tag != -1 and self.regions[tag].is_crossroad():
+                result.append(self.regions[tag].id)
+            else:
+                result.append(-1)
+        return result
+
+    def is_crossroad_node(self, n):
+        tag = self.G.nodes[n][rg.Region.label_region]
+        if tag == -1:
+            False
+        else: 
+            return self.regions[tag].is_crossroad()
+
+    ######################### Functions used to prepare the graph ########################
+
+    def remove_footways_and_parkings(G, keep_all_components):
+            
+
+        # remove footways and parkings
+        to_remove = []
+        for u, v, a in G.edges(data = True):
+            if "footway" in a or ("highway" in a and a["highway"] in ["footway"]):
+                to_remove.append((u, v))
+                # add missing crossings
+                if not "highway" in G.nodes[u]:
+                    G.nodes[u]["highway"] = "crossing"
+                if not "highway" in G.nodes[v]:
+                    G.nodes[v]["highway"] = "crossing"
+            if ("highway" in a and a["highway"] in ["cycleway", "path", "pedestrian", "steps"]):
+                to_remove.append((u, v))
+            #elif "service" in a and a["service"] in ["parking_aisle"]:
+            #    to_remove.append((u, v))                
+        G.remove_edges_from(to_remove)
+        G = ox.utils_graph.remove_isolated_nodes(G)
+        if not keep_all_components and len(G.nodes) != 0:
+            G = ox.utils_graph.get_largest_component(G)
+        return G
+
+        
+    ######################### Functions related to graph rendering (colors) ########################
+
+    # return edge colors according to the region label
+    def get_regions_colors(self):
+        result = {}
+        color = {}
+        for e in self.G.edges:
+            tag = self.G[e[0]][e[1]][e[2]][rg.Region.label_region]
+            if tag == -1:
+                result[e] = (0.5, 0.5, 0.5, 0.5)
+            else:
+                if not tag in color:
+                    color[tag] = Segmentation.random_color()
+                result[e] = color[tag]
+        return pd.Series(result)
+
+    # return edge colors according to the region class label
+    def get_regions_class_colors(self):
+        result = {}
+        for e in self.G.edges:
+            tag = self.G[e[0]][e[1]][e[2]][rg.Region.label_region]
+            if tag == -1:
+                result[e] = (0.5, 0.5, 0.5, 0.5)
+            elif self.regions[tag].is_crossroad():
+                result[e] = (0.8, 0, 0, 1)
+            elif self.regions[tag].is_branch():
+                result[e] = (0.6, 0.6, 0, 1)
+            else:
+                result[e] = (0.3, 0.3, 0.3, 1)
+
+        return pd.Series(result)
+
+
+    def random_color(only_bg = False):
+        r1 = math.pi * random.random()
+        r2 = math.pi * random.random()
+        start = 0.2
+        coef = 0.8
+        return (0 if only_bg else (start + coef * abs(math.sin(r1)) * abs(math.sin(r2))), \
+                start + coef * abs(math.cos(r1)) * abs(math.sin(r2)), \
+                start + coef * abs(math.sin(r1)) * abs(math.cos(r2)), 
+                1)
+
+    # return edge colors using one random color per label
+    def get_edge_random_colors_by_attr(G, label, values = {}):
+        result = {}
+        for e in G.edges:
+            tag = G[e[0]][e[1]][e[2]][label]
+            if not tag in values:
+                if tag == -1:
+                    values[tag] = (0.5, 0.5, 0.5, 0.5)
+                else:
+                    values[tag] = Segmentation.random_color()
+            result[e] = values[tag]
+        return pd.Series(result)
+        
+    def get_nodes_reliability_on_regions_colors(self):
+
+        result = {}
+        for n in self.G.nodes:
+            r_class = rel.Reliability.get_best_reliability_node(self.G, n)
+            r_value = self.G.nodes[n][r_class]
+            coef = (r_value - rel.Reliability.strongly_no) / (rel.Reliability.strongly_yes - rel.Reliability.strongly_no)
+            coef = math.pow(coef, 2)
+            if r_class == rel.Reliability.crossroad_reliability:
+                result[n] = (0.8, 0, 0, coef)
+            elif r_class == rel.Reliability.boundary_reliability:
+                adj = self.get_adjacent_crossroad_regions(n)
+                # in the middle of a branch
+                if len([n for n in adj if n != -1]) == 0:
+                    result[n] = (0, 0, 0.6, coef)
+                # inside a region
+                elif len(list(set([n for n in adj if n != -1]))) == 1 and len([n for n in adj if n != -1]) != 1:
+                    result[n] = (1, 0.6, 0.6, coef)
+                # in a boundary
+                else:
+                    result[n] = (0.6, 0.6, 0, coef)
+
+            else: # branch
+                result[n] = (0, 0, 0, 0)
+
+
+        return pd.Series(result)
+
+
+    def get_edges_reliability_colors(self):
+        result = {}
+        for e in self.G.edges:
+            r_value = self.G[e[0]][e[1]][e[2]][rel.Reliability.crossroad_reliability]
+            coef = (r_value - rel.Reliability.strongly_no) / (rel.Reliability.strongly_yes - rel.Reliability.strongly_no)
+            coef = math.pow(coef, 2)
+            result[e] = (1, 1, 1, coef)
+        return pd.Series(result)
+
+    def get_nodes_reliability_colors(self):
+
+        result = {}
+        for n in self.G.nodes:
+            r_class = rel.Reliability.get_best_reliability_node(self.G, n)
+            r_value = self.G.nodes[n][r_class]
+            coef = (r_value - rel.Reliability.strongly_no) / (rel.Reliability.strongly_yes - rel.Reliability.strongly_no)
+            coef = math.pow(coef, 2)
+            if r_class == rel.Reliability.boundary_reliability:
+                result[n] = (0.1, 0, 0.8, coef)
+            else:
+                result[n] = (0.8, 0, 0, coef)
+
+
+        return pd.Series(result)
+
+    def get_boundary_node_colors(self):
+
+        result = {}
+        for n in self.G.nodes:
+            nb_adj_crossings = len(list(set([r for r in self.get_adjacent_crossroad_regions(n) if r != -1])))
+            nbnb = len(list(self.G.neighbors(n)))
+            nbAdj = len([ nb for nb in self.G.neighbors(n) if rg.Region.unknown_region_edge_in_graph(self.G, (n, nb))])
+            if nbnb == nbAdj:
+                if nbnb == 1: # dead end
+                    result[n] = (0.5, 0.5, 0.5, 0.1)
+                elif rg.Region.unknown_region_node_in_graph(self.G, n):
+                    result[n] = (0, 0, 0.5, 1) # node not taggued, possibly a missing crossing
+                else:
+                    if nb_adj_crossings >= 2:
+                        result[n] = (0.6, 0.6, 0, 1) # splitter in a crossroad
+                    elif nb_adj_crossings == 0 and self.is_crossroad_node(n):
+                        result[n] = (1, 0, 0, 1) # single-node crossroad
+                    else:
+                        result[n] = (0, 0, 0, 0)
+            else:
+                if nb_adj_crossings >= 2:
+                    result[n] = (0.6, 0.6, 0, 1) # splitter in a crossroad
+                elif nb_adj_crossings == 0 and self.is_crossroad_node(n):
+                    result[n] = (1, 0, 0, 1) # single-node crossroad
+                else:
+                    result[n] = (0, 0, 0, 0)
+        return pd.Series(result)
+
+    def get_nodes_regions_colors(self):
+        result = {}
+        for n in self.G.nodes:
+            if len(list(self.G.neighbors(n))) <= 2:
+                result[n] = (0, 0, 0, 0)
+            else:
+                label = self.G.nodes[n][rg.Region.label_region]
+                if label < 0:
+                    result[n] = (0, 0, 0, 0)
+                else:
+                    nb_edge_in_region = len([nb for nb in self.G[n] if self.G[n][nb][0][rg.Region.label_region] == label])
+                    if nb_edge_in_region == 0:
+                        result[n] = Segmentation.random_color()
+                    else:
+                        result[n] = (0, 0, 0, 0)
+
+        return pd.Series(result)
+
+    # input: a list of crossroads (main crossroad and possibly contained crossroads)
+    def get_regions_colors_from_crossroad(self, cr):
+        crids = [ c.id for c in cr ]
+        mainCR = max(cr, key=lambda x: len(x.nodes))
+        maxID = max(crids)
+        result = {}
+        color = {}
+        for e in self.G.edges:
+            tag = self.G[e[0]][e[1]][e[2]][rg.Region.label_region]
+            if not tag in crids:
+                # check if it's a branch of the main crossroad
+                bid = mainCR.get_branch_id(e)
+                if bid == -1:
+                    result[e] = (0.5, 0.5, 0.5, 0.1)
+                else:
+                    tag = maxID + bid + 1
+                    if not tag in color:
+                        color[tag] = Segmentation.random_color()
+                    result[e] = color[tag]
+            else:
+                ncrs = len([ c for c in cr if c.has_edge(e) ])
+                result[e] = (math.sqrt(ncrs / len(crids)), 0, 0, 1)
+        return pd.Series(result)
+
+    # input: a list of crossroads (main crossroad and possibly contained crossroads)
+    def get_nodes_regions_colors_from_crossroad(self, cr):
+        crids = [ c.id for c in cr ]
+        mainCR = max(cr, key=lambda x: len(x.nodes))
+        result = {}
+        for n in self.G.nodes:
+            if len(list(self.G.neighbors(n))) <= 2:
+                result[n] = (0, 0, 0, 0)
+            else:
+                label = self.G.nodes[n][rg.Region.label_region]
+                if not label in crids:
+                    result[n] = (0, 0, 0, 0)
+                else:
+                    # get regions that contains this node with no adjacent edge
+                    ncrn = len([ r for r in cr if r.has_node(n) and len(r.edges_with_node(n)) == 0])
+                    if ncrn == 0:
+                        result[n] = (0, 0, 0, 0)
+                    else:
+                        result[n] = (math.sqrt(ncrn / len(crids)), 0, 0, 1)
+
+        return pd.Series(result)
+
+    ######################### text descriptions ########################
+
+    # return a list of crossroads (main crossroad and possibly contained crossroads)
+    def get_crossroad(self, longitude, latitude, multiscale = False):
+        distance = -1
+        middle = -1
+        for rid in self.regions:
+            if self.regions[rid].is_crossroad():
+                region = self.regions[rid]
+                x1 = self.G.nodes[region.get_center()]["x"]
+                y1 = self.G.nodes[region.get_center()]["y"]
+                d = ox.distance.great_circle_vec(lat1=y1, lng1=x1, lat2=latitude, lng2=longitude)
+                if distance < 0 or d < distance:
+                    distance = d
+                    middle = rid
+
+        if multiscale:
+            result = []
+            result.append(self.regions[middle])
+            for rid in self.inner_regions:
+                if self.regions[middle].contains(self.inner_regions[rid]):
+                    result.append(self.inner_regions[rid])
+            return result
+        else:
+            return [self.regions[middle]]
+
+    def to_text(self, longitude, latitude, multiscale = False):
+        cs = self.get_crossroad(longitude, latitude, multiscale)
+        result = ""
+        for i, c in enumerate(cs):
+            if i != 0:
+                result += "\n\n"
+            result += c.to_text()
+        return result
+
+    def to_text_all(self, multiscale = False):
+        result = ""
+        for rid in self.regions:
+            if self.regions[rid].is_crossroad():
+                result += self.regions[rid].to_text()
+                result += "\n"
+                result += "\n"
+
+        if multiscale:
+            result = "Inner crossroads:"
+            result += "\n"
+            for rid in self.inner_regions:
+                result += self.inner_regions[rid].to_text()
+                result += "\n"
+                result += "\n"
+        return result
+
+    ######################### json descriptions ########################
+
+    def to_json(self, filename, longitude, latitude, multiscale = False):
+        data = [x.to_json_data() for x in self.get_crossroad(longitude, latitude, multiscale)]
+
+        with open(filename, 'w') as outfile:
+            json.dump(data, outfile)
+
+
+    def to_json_all(self, filename, multiscale = False):
+        data = []
+        for rid in self.regions:
+            if self.regions[rid].is_crossroad():
+                entry = self.regions[rid].to_json_data()
+                data.append(entry)
+        
+        if multiscale:
+            for rid in self.inner_regions:
+                entry = self.inner_regions[rid].to_json_data()
+                data.append(entry)
+
+        with open(filename, 'w') as outfile:
+            json.dump(data, outfile)
+
+    
+    ######################### geopackage description ########################
+
+    def to_geopackage(self, filename):
+        # add new attributes
+        nx.set_node_attributes(self.G, values = "", name = "crossroad")
+        nx.set_node_attributes(self.G, values = "", name = "sub_crossroad")
+        nx.set_edge_attributes(self.G, values = "", name = "crossroad")
+        nx.set_edge_attributes(self.G, values = "", name = "sub_crossroad")
+        nx.set_edge_attributes(self.G, values = "", name = "branch")
+
+        # fill attributes according to the computed regions
+        for rid in self.regions:
+            if self.regions[rid].is_crossroad():
+                self.regions[rid].set_graph_attributes("crossroad", "branch")
+
+        for rid in self.inner_regions:
+            self.inner_regions[rid].set_graph_attributes("sub_crossroad")
+
+        # save graph
+        ox.io.save_graph_geopackage(self.G, filename)
+
+
diff --git a/crseg/utils.py b/crseg/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..8d739d123cf6ca5af8464cdcc75e9c4017ee1a98
--- /dev/null
+++ b/crseg/utils.py
@@ -0,0 +1,180 @@
+import osmnx as ox
+
+from . import region as r
+
+class Util:
+
+
+    def centroid(G, points):
+        x = 0.0
+        y = 0.0
+        for p in points:
+            x += G.nodes[p]["x"]
+            y += G.nodes[p]["y"]
+        return (x / len(points), y / len(points))
+
+    def coords_distance(point1, point2):
+        x1 = point1[0]
+        y1 = point1[1]
+        x2 = point2[0]
+        y2 = point2[1]
+        return ox.distance.great_circle_vec(lat1=y1, lng1=x1, lat2=y2, lng2=x2)
+
+    def distance_to(G, node, point):
+        x1 = G.nodes[node]["x"]
+        y1 = G.nodes[node]["y"]
+        x2 = point[0]
+        y2 = point[1]
+        return ox.distance.great_circle_vec(lat1=y1, lng1=x1, lat2=y2, lng2=x2)
+
+    # links are shorter than real paths
+    def distance_with_shortcut(G, node1, node2):
+        gEdge = G[node1][node2][0]
+        coef = 1
+        if "highway" in gEdge and gEdge["highway"] in ["primary_link", "secondary_link", "tertiary_link", "trunk_link", "motorway_link"]:
+            coef = 0.5
+        if "junction" in gEdge:
+            coef = 0.5
+        return Util.distance(G, node1, node2) * coef
+
+    def distance(G, node1, node2):
+        x1 = G.nodes[node1]["x"]
+        y1 = G.nodes[node1]["y"]
+        x2 = G.nodes[node2]["x"]
+        y2 = G.nodes[node2]["y"]
+        return ox.distance.great_circle_vec(lat1=y1, lng1=x1, lat2=y2, lng2=x2)
+
+    def bearing(G, node1, node2):
+        x1 = G.nodes[node1]["x"]
+        y1 = G.nodes[node1]["y"]
+        x2 = G.nodes[node2]["x"]
+        y2 = G.nodes[node2]["y"]
+        return ox.bearing.get_bearing((y1, x1), (y2, x2))
+
+    def length(G, path):
+        return sum([Util.distance(G, p1, p2) for p1, p2 in zip(path, path[1:])])
+
+    def length_with_shortcut(G, path):
+        return sum([Util.distance_with_shortcut(G, p1, p2) for p1, p2 in zip(path, path[1:])])
+
+    def angular_distance(angle1, angle2):
+        a = angle1 - angle2
+        if a > 180:
+            a -= 360
+        if a < -180:
+            a += 360 
+        return abs(a)
+
+    def is_inside_parking(G, node):
+        for nb in G.neighbors(node):
+            if (not "service" in G[node][nb][0]) or (G[node][nb][0]["service"] != "parking_aisle"):
+                return False
+        return True
+
+    def get_adjacent_streetnames(G, node):
+        streetnames = set()
+        for nb in G.neighbors(node):
+            if "name" in G[node][nb][0]:
+                streetnames.add(G[node][nb][0]["name"])
+            elif "ref" in G[node][nb][0]:
+                streetnames.add(G[node][nb][0]["ref"])
+            else:
+                streetnames.add(None)
+        return list(streetnames)
+
+    def is_biffurcation(G, n):
+        return len(list(G.neighbors(n))) > 2
+
+    def is_middle_polyline(G, n):
+        return len(list(G.neighbors(n))) == 2
+
+    def get_opposite_node(G, n, other):
+        for nb in G.neighbors(n):
+            if nb != other:
+                return nb
+        # will not append
+        return None
+
+
+    def get_path_to_biffurcation(G, n1, n2, max = -1):
+        path = [n1, n2]
+        length = Util.distance(G, n1, n2)
+
+        while (max < 0 or length < max) and Util.is_middle_polyline(G, path[len(path) - 1]):
+            path.append(Util.get_opposite_node(G, path[len(path) - 1], path[len(path) - 2]))
+            length += Util.distance(G, path[len(path) - 2], path[len(path) - 1])
+        
+        return path
+
+    # return true if two the node is part of 3 edges, and
+    # if two of them are one-way
+    def is_street_separation(G, n):
+        if len(G[n]) != 3:
+            return False
+        
+        return len([nb for nb in G.neighbors(n) if "oneway" in G[n][nb][0] and G[n][nb][0]["oneway"]]) >= 2
+
+    def is_part_of_local_triangle(G, n, max_perimeter = 150):
+
+        paths = [ Util.get_path_to_biffurcation(G, n, nb) for nb in G.neighbors(n)]
+
+        for i1, p1 in enumerate(paths):
+
+            p1_end = p1[len(p1) - 1]
+            p1_end_paths = [ Util.get_path_to_biffurcation(G, p1_end, nb) for nb in G.neighbors(p1_end)]
+            p1_end_neighbors = [ p[len(p) - 1] for p in p1_end_paths]
+
+            for i2 in range(i1, len(paths)):
+                p2 = paths[i2]
+                p2_end = p2[len(p2) - 1]
+                if p2_end in p1_end_neighbors:
+                    p = [ path for path in p1_end_paths if path[len(path) - 1] == p2_end][0]
+                    l = Util.length(G, p1) + Util.length(G, p2) + Util.length(G, p)
+                    if l < max_perimeter:
+                        return True
+
+        return False
+
+    def has_non_labeled_adjacent_edge(G, n):
+        for nb in G.neighbors(n):
+            if G[n][nb][0][r.Region.label_region] == -1:
+                return True
+        return False
+
+    def estimate_edge_width(G, edge):
+        gEdge = G[edge[0]][edge[1]][0]
+        import re
+        if "width" in gEdge and not re.match(r'^-?\d+(?:\.\d+)$', gEdge["width"]) is None:
+            return float(gEdge["width"])
+        elif "lanes" in gEdge:
+            nb = int(gEdge["lanes"])
+        else:
+            if "oneway" in gEdge and gEdge["oneway"]:
+                nb = 1
+            else:
+                nb = 2
+        
+        if "highway" in gEdge:
+            if gEdge["highway"] in ["motorway", "trunk"]:
+                width = 3.5
+            elif gEdge["highway"] in ["primary"]:
+                width = 3
+            elif gEdge["highway"] in ["secondary"]:
+                width = 2.75
+            elif gEdge["highway"] in ["service"]:
+                width = 2.25
+            else:
+                width = 2.75
+        else:
+            width = 3
+        
+        result = 0
+        # TODO: improve integration of cycleways in this computation
+        if ("cycleway" in gEdge and gEdge["cycleway"] == "track") or \
+            ("cycleway:left" in gEdge and gEdge["cycleway:left"] == "track") or \
+            ("cycleway:right" in gEdge and gEdge["cycleway:right"] == "track"):
+            result = (nb + 1) * width # ~ COVID tracks
+        else:
+            result = nb * width
+
+        return result
diff --git a/dist/crossroads-segmentation-0.1.1.tar.gz b/dist/crossroads-segmentation-0.1.1.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..8557a2321c8d1afb981a41a71a879620052ea813
Binary files /dev/null and b/dist/crossroads-segmentation-0.1.1.tar.gz differ
diff --git a/dist/crossroads_segmentation-0.1.1-py3-none-any.whl b/dist/crossroads_segmentation-0.1.1-py3-none-any.whl
new file mode 100644
index 0000000000000000000000000000000000000000..ba7191ef164be52a30da75cca3b2b395ef3eb126
Binary files /dev/null and b/dist/crossroads_segmentation-0.1.1-py3-none-any.whl differ
diff --git a/src/.gitignore b/examples/.gitignore
similarity index 100%
rename from src/.gitignore
rename to examples/.gitignore
diff --git a/src/crossroads-by-name.json b/examples/crossroads-by-name.json
similarity index 100%
rename from src/crossroads-by-name.json
rename to examples/crossroads-by-name.json
diff --git a/src/get-crossroad-description.py b/examples/get-crossroad-description.py
similarity index 100%
rename from src/get-crossroad-description.py
rename to examples/get-crossroad-description.py
diff --git a/src/get-paris-streets.py b/examples/get-paris-streets.py
similarity index 100%
rename from src/get-paris-streets.py
rename to examples/get-paris-streets.py
diff --git a/src/stats-intersection-info.py b/examples/stats-intersection-info.py
similarity index 100%
rename from src/stats-intersection-info.py
rename to examples/stats-intersection-info.py
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..457e40cc41cb3674781e7215a24f45b3e39e6dee
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1 @@
+osmnx
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..5e42f21a2bec68fbf096dfb5c047f5696230ae9b
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,28 @@
+import pathlib
+from setuptools import setup
+
+# The directory containing this file
+HERE = pathlib.Path(__file__).parent
+
+# The text of the README file
+README = (HERE / "README.md").read_text()
+
+# This call to setup() does all the work
+setup(
+    name="crossroads-segmentation",
+    version="0.1.1",
+    description="Crossroads segmentation is a python tool that produces automatic segmentations of data from OpenStreetMap.",
+    long_description=README,
+    long_description_content_type="text/markdown",
+    url="https://gitlab.limos.fr/jmafavre/crossroads-segmentation/",
+    author="Jean-Marie Favreau",
+    author_email="j-marie.favreau@uca.fr",
+    license="AGPL-3.0",
+    classifiers=[
+        "License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)",
+        "Programming Language :: Python :: 3",
+        "Programming Language :: Python :: 3.7",
+    ],
+    packages=["crseg"],
+    include_package_data=True,
+)
diff --git a/src/crseg/__init__.py b/src/crseg/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000