score:15
create a graph in which the nodes correspond to delaunay triangles and in which there is a graph edge between two triangles if and only if they share two vertices.
compute the connected components of the graph.
for each connected component, find all of the nodes that have less than three adjacent nodes (that is, those that have degree 0, 1, or 2). these correspond to the boundary triangles. we define the edges of a boundary triangle that are not shared with another triangle to be the boundary edges of that boundary triangle.
as an example, i have highlighted these boundary triangles in your example "question mark" delaunay triangulation:
by definition, every boundary triangle is adjacent to at most two other boundary triangles. the boundary edges of boundary triangles form cycles. you can simply traverse those cycles to determine polygon shapes for the boundaries. this will also work for polygons with holes if you keep them in mind in your implementation.
score:-1
alpha shapes is defined as a delaunay triangulation without edges exceeding alpha. first of remove all interior triangles and then all edges exceeding alpha.
score:0
it turns out topojson has a merge algorithm which performs just this task: https://github.com/mbostock/topojson/wiki/api-reference#merge
there's even an example showing it in action: http://bl.ocks.org/mbostock/9927735
in my case, it was easy for me to generate topojson data, and this library function accomplished the task perfectly for me.
score:0
building up on @timothy's answer, i used the following algorithm to calculate the boundary rings of a delaunay triangulation.
from matplotlib.tri import triangulation
import numpy as np
def get_boundary_rings(x, y, elements):
mpl_tri = triangulation(x, y, elements)
idxs = np.vstack(list(np.where(mpl_tri.neighbors == -1))).t
unique_edges = list()
for i, j in idxs:
unique_edges.append((mpl_tri.triangles[i, j],
mpl_tri.triangles[i, (j+1) % 3]))
unique_edges = np.asarray(unique_edges)
ring_collection = list()
initial_idx = 0
for i in range(1, len(unique_edges)-1):
if unique_edges[i-1, 1] != unique_edges[i, 0]:
try:
idx = np.where(
unique_edges[i-1, 1] == unique_edges[i:, 0])[0][0]
unique_edges[[i, idx+i]] = unique_edges[[idx+i, i]]
except indexerror:
ring_collection.append(unique_edges[initial_idx:i, :])
initial_idx = i
continue
# if there is just one ring, the exception is never reached,
# so populate ring_collection before returning.
if len(ring_collection) == 0:
ring_collection.append(np.asarray(unique_edges))
return ring_collection
score:1
slightly revise hanniel's answer for 3d point case (tetrahedron).
def alpha_shape(points, alpha, only_outer=true):
"""
compute the alpha shape (concave hull) of a set of points.
:param points: np.array of shape (n, 3) points.
:param alpha: alpha value.
:param only_outer: boolean value to specify if we keep only the outer border
or also inner edges.
:return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
the indices in the points array.
"""
assert points.shape[0] > 3, "need at least four points"
def add_edge(edges, i, j):
"""
add a line between the i-th and j-th points,
if not in the set already
"""
if (i, j) in edges or (j, i) in edges:
# already added
if only_outer:
# if both neighboring triangles are in shape, it's not a boundary edge
if (j, i) in edges:
edges.remove((j, i))
return
edges.add((i, j))
tri = delaunay(points)
edges = set()
# loop over triangles:
# ia, ib, ic, id = indices of corner points of the tetrahedron
print(tri.vertices.shape)
for ia, ib, ic, id in tri.vertices:
pa = points[ia]
pb = points[ib]
pc = points[ic]
pd = points[id]
# computing radius of tetrahedron circumsphere
# http://mathworld.wolfram.com/circumsphere.html
pa2 = np.dot(pa, pa)
pb2 = np.dot(pb, pb)
pc2 = np.dot(pc, pc)
pd2 = np.dot(pd, pd)
a = np.linalg.det(np.array([np.append(pa, 1), np.append(pb, 1), np.append(pc, 1), np.append(pd, 1)]))
dx = np.linalg.det(np.array([np.array([pa2, pa[1], pa[2], 1]),
np.array([pb2, pb[1], pb[2], 1]),
np.array([pc2, pc[1], pc[2], 1]),
np.array([pd2, pd[1], pd[2], 1])]))
dy = - np.linalg.det(np.array([np.array([pa2, pa[0], pa[2], 1]),
np.array([pb2, pb[0], pb[2], 1]),
np.array([pc2, pc[0], pc[2], 1]),
np.array([pd2, pd[0], pd[2], 1])]))
dz = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], 1]),
np.array([pb2, pb[0], pb[1], 1]),
np.array([pc2, pc[0], pc[1], 1]),
np.array([pd2, pd[0], pd[1], 1])]))
c = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], pa[2]]),
np.array([pb2, pb[0], pb[1], pb[2]]),
np.array([pc2, pc[0], pc[1], pc[2]]),
np.array([pd2, pd[0], pd[1], pd[2]])]))
circum_r = math.sqrt(math.pow(dx, 2) + math.pow(dy, 2) + math.pow(dz, 2) - 4 * a * c) / (2 * abs(a))
if circum_r < alpha:
add_edge(edges, ia, ib)
add_edge(edges, ib, ic)
add_edge(edges, ic, id)
add_edge(edges, id, ia)
add_edge(edges, ia, ic)
add_edge(edges, ib, id)
return edges
score:3
there now exists a python package alphashape which is extremely easy to use, and can be installed by pip
or conda
.
the main function has similar inputs to the answer given by @iddo hanniel, adjusting the second positional argument would give you the desired outline. alternatively, you could leave the seconda positional argument blank and the function would optimize that parameter for you to give you the best concave hull. beware, the computational time is increased greatly if you let the function optimize the value.
score:3
i know it is a delayed answer, but the methods posted here did not work for me for various reasons.
the package alphashape that is mentioned is generally good but its downside is that it uses shapely's cascade_union
and triangulation methods to give you the final concave hull which is super slow for large datasets and low alpha values (high precision). in this case the code posted by iddo hanniel (and the revision by harold) will keep a great number of edges on the interior and the only way to dissolve them is to use the aforementioned cascade_union
and triangulation from shapely.
i generally work with complex forms and the code below works fine and it is fast (2 seconds for 100k 2d points):
import numpy as np
from shapely.geometry import multilinestring
from shapely.ops import unary_union, polygonize
from scipy.spatial import delaunay
from collections import counter
import itertools
def concave_hull(coords, alpha): # coords is a 2d numpy array
# i removed the qbb option from the scipy defaults.
# it is much faster and equally precise without it.
# unless your coords are integers.
# see http://www.qhull.org/html/qh-optq.htm
tri = delaunay(coords, qhull_options="qc qz q12").vertices
ia, ib, ic = (
tri[:, 0],
tri[:, 1],
tri[:, 2],
) # indices of each of the triangles' points
pa, pb, pc = (
coords[ia],
coords[ib],
coords[ic],
) # coordinates of each of the triangles' points
a = np.sqrt((pa[:, 0] - pb[:, 0]) ** 2 + (pa[:, 1] - pb[:, 1]) ** 2)
b = np.sqrt((pb[:, 0] - pc[:, 0]) ** 2 + (pb[:, 1] - pc[:, 1]) ** 2)
c = np.sqrt((pc[:, 0] - pa[:, 0]) ** 2 + (pc[:, 1] - pa[:, 1]) ** 2)
s = (a + b + c) * 0.5 # semi-perimeter of triangle
area = np.sqrt(
s * (s - a) * (s - b) * (s - c)
) # area of triangle by heron's formula
filter = (
a * b * c / (4.0 * area) < 1.0 / alpha
) # radius filter based on alpha value
# filter the edges
edges = tri[filter]
# now a main difference with the aforementioned approaches is that we dont
# use a set() because this eliminates duplicate edges. in the list below
# both (i, j) and (j, i) pairs are counted. the reasoning is that boundary
# edges appear only once while interior edges twice
edges = [
tuple(sorted(combo)) for e in edges for combo in itertools.combinations(e, 2)
]
count = counter(edges) # count occurrences of each edge
# keep only edges that appear one time (concave hull edges)
edges = [e for e, c in count.items() if c == 1]
# these are the coordinates of the edges that comprise the concave hull
edges = [(coords[e[0]], coords[e[1]]) for e in edges]
# use this only if you need to return your hull points in "order" (i think
# its ccw)
ml = multilinestring(edges)
poly = polygonize(ml)
hull = unary_union(list(poly))
hull_vertices = hull.exterior.coords.xy
return edges, hull_vertices
score:17
here is some python code that does what you need. i modified the alpha-shape (concave hull) computation from here so that it doesn't insert inner edges (the only_outer
parameter). i also made it self-contained so it doesn't depend on an outside library.
from scipy.spatial import delaunay
import numpy as np
def alpha_shape(points, alpha, only_outer=true):
"""
compute the alpha shape (concave hull) of a set of points.
:param points: np.array of shape (n,2) points.
:param alpha: alpha value.
:param only_outer: boolean value to specify if we keep only the outer border
or also inner edges.
:return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
the indices in the points array.
"""
assert points.shape[0] > 3, "need at least four points"
def add_edge(edges, i, j):
"""
add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
assert (j, i) in edges, "can't go twice over same directed edge right?"
if only_outer:
# if both neighboring triangles are in shape, it is not a boundary edge
edges.remove((j, i))
return
edges.add((i, j))
tri = delaunay(points)
edges = set()
# loop over triangles:
# ia, ib, ic = indices of corner points of the triangle
for ia, ib, ic in tri.simplices:
pa = points[ia]
pb = points[ib]
pc = points[ic]
# computing radius of triangle circumcircle
# www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
s = (a + b + c) / 2.0
area = np.sqrt(s * (s - a) * (s - b) * (s - c))
circum_r = a * b * c / (4.0 * area)
if circum_r < alpha:
add_edge(edges, ia, ib)
add_edge(edges, ib, ic)
add_edge(edges, ic, ia)
return edges
if you run it with the following test code you will get this figure:
from matplotlib.pyplot import *
# constructing the input point data
np.random.seed(0)
x = 3.0 * np.random.rand(2000)
y = 2.0 * np.random.rand(2000) - 1.0
inside = ((x ** 2 + y ** 2 > 1.0) & ((x - 3) ** 2 + y ** 2 > 1.0) & ((x - 1.5) ** 2 + y ** 2 > 0.09))
points = np.vstack([x[inside], y[inside]]).t
# computing the alpha shape
edges = alpha_shape(points, alpha=0.25, only_outer=true)
# plotting the output
figure()
axis('equal')
plot(points[:, 0], points[:, 1], '.')
for i, j in edges:
plot(points[[i, j], 0], points[[i, j], 1])
show()
Source: stackoverflow.com
Related Query
- Calculate bounding polygon of alpha shape from the Delaunay triangulation
- calculate the date diff in javascript from JSON array
- nvD3 - multiBarChart - How to start from 0 and how to change the shape of controls
- Getting the centroid of each polygon from a GeoJson multipolygon in d3.js
- How do I calculate a bounding box for leaflet from a set of points in d3
- Calculate how many std deviations the values of certain keys are from the mean
- d3.js Lasso Drawing Polygon Shape Search Tool - on a google map (getting the coordinates/users in a given area)
- function to calculate the parent size from children's size
- Accessing d3.js element attributes from the datum?
- Calculate width of text before drawing the text
- clicking a node in d3 from a button outside the svg
- In d3, how to get the interpolated line data from a SVG line?
- Can d3.js draw two scatterplots on the same graph using data from the same source?
- Problems converting from shape to topojson
- Calculate the exact size of a font glyph
- Why can't I get the Bounding Box of this d3.js text?
- D3 bubble chart / pack layout - How to make bubbles radiate out from the largest bubbles to the smallest?
- How to calculate the scale's result?
- D3: Finding the area of a geo polygon in d3
- d3.js geo worldmap - merge russia (shift small part from the left next to america to the right)
- How to limit a d3 line-chart from showing the line outside of the range of the axis?
- Updating the data of a pack layout from JSON call and redrawing
- D3: get the bounding box of a selected element
- enclosing the nodes of a d3 force directed graph in a circle or a polygon or a cloud
- How to update elements of an HTML that the elements are created using data from a CSV file?
- How can I remove a line from the 110m TopoJson world map?
- How to assign the x-axis position of a node in a Sankey Diagram (D3) from the json file
- Display only values from the data set into X axis ticks
- d3 setting multiple attributes from the same function?
- How to get the data from the c3.js
More Query from same tag
- D3.js - how selections work - Need clarification on Mike's article
- Beginner in d3, making US county map
- d3.js scatter chart - clipping the sparklines
- How to get rid of hairline from overlapping svg shapes
- D3 Chart Legend, Bars not appearing
- Uncaught SyntaxError: Unexpected token <
- D3 zoomable sunburst order children
- Return value from JSON, not the children count
- d3js v4 x-axis-label is there but not visible
- Text inside Bar of the Grouped Bar Chart D3.js
- For loop does not work in D3.js
- Cannot plot polyline on D3 spinning globe
- Append different shapes
- Plot polygon using D3
- What HTML processing is used on observablehq?
- C3.js No Data option Not displaying
- How to animate elements for every swap in a array using d3.js
- Transition Between Curve Types Using D3.js
- How to have d3 tickFormat of (day, time) but only show time if day is the same?
- how to format time on xAxis use d3.js
- d3pie.Js Pie chart error to get the script started
- d3.js changing style of diagonal line
- D3 selection - reflecting data changes
- DC.js Stacked Line Chart not working
- Axis issue in a D3 punchcard - need data to extend past axis limit
- D3JS line chart inverted issue
- How do I resolve "[ERR_REQUIRE_ESM]: Must use import to load ES Module" when using D3.js 7.0.0 with Next.js 11.0.1?
- Using a scale in d3.js to fill up a vertical bar with 3 colors
- D3.js: rotate coordinates, not an entire element?
- c3.js how to set json datas in x/y Axis