82
The undocumented Graphics`PolygonUtils`PointWindingNumber (if you're on versions < 10, use Graphics`Mesh`PointWindingNumber) does exactly this — it gives you the winding number of a point. A point lies inside the polygon if and only if its winding number is non-zero.
Using this, you can create a Boolean function to test if a point is inside the polygon
...
52
Using the function winding from Heike's answer to a related question
winding[poly_, pt_] :=
Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly),
2 Pi, -Pi]/2/Pi)]
to modify the test function in this Wolfram Demonstration by R. Nowak to
testpoint[poly_, pt_] :=
Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt -...
49
Edit: added Gradient -> grad[vars] option. Without this small option the code was several orders of magnitude slower.
Yes, it can! Unfortunately, not automatically.
There are different algorithms to do it (see special literature, e.g. Dziuk, Gerhard, and John E. Hutchinson. A finite element method for the computation of parametric minimal surfaces. ...
44
======= Update =========
Great question! It inspired this Wolfram Blog article and includes most of the code below plus some apps and fractal layouts like this:
I think it make sense to keep the older code blow for archival and historic purposes.
======= Older implementation =========
Excellent motivating creativity question. This is a bit big for a ...
44
Obtain the image:
i = Import["http://i.stack.imgur.com/iab6u.png"];
Compute the distance transform:
k = DistanceTransform[ColorNegate[i]] // ImageAdjust;
ReliefPlot[Reverse@ImageData[k]] (* To illustrate *)
Identify the "peaks," which must bound the Voronoi cells:
l = ColorNegate[Binarize[ColorNegate[LaplacianGaussianFilter[k, 2] // ...
41
We need quite a bit of preparation. In the first place we need methods to compute cell adjacency matrices from here. I copied the code for completeness.
CellAdjacencyMatrix[R_MeshRegion, d_, 0] := If[MeshCellCount[R, d] > 0,
Unitize[R["ConnectivityMatrix"[d, 0]]],
{}
];
CellAdjacencyMatrix[R_MeshRegion, 0, d_] := If[MeshCellCount[R, d] > 0,
...
40
Version 11 has both symbolic and numeric eigensolvers, see here for an overview
Here is a slightly different way to do it. We write a function that converts any PDE (1D/2D/3D) into discretized system matices:
Needs["NDSolve`FEM`"]
PDEtoMatrix[{pde_, Γ___}, u_, r__,
o : OptionsPattern[NDSolve`ProcessEquations]] :=
Module[{ndstate, feData, sd, bcData, ...
39
I guess the first step would always be to find an ordered list of points along the middle of the curve. That I can help with:
First binarize and thin the image of the curve, so you get a 1-pixel wide white line:
img = Import["http://i.stack.imgur.com/fEf1i.jpg"];
bin = Thinning@ColorNegate[Binarize[img]]
Finding the white pixels in this image is easy:
...
32
Well, you can use the undocumented RegionDistance which does exactly this as follows: (This answer, as written, only works for V9 as noted by Oska, for V10 see update below)
here is a triangle in 3D
region = Polygon[{{0, 0, 0}, {1, 0, 0}, {0, 1, 1}}];
Graphics3D[region]
Now suppose you want to find the shortest distance from the point {1, 1, 1} in 3D to ...
32
Note added 1/29/2020: the routines here have a bug where the mean curvature is sometimes computed with the opposite sign. I still need to work on how to fix this.
I guess I should not have been surprised that there are actually many ways to estimate the Gaussian and mean curvature of a triangular mesh. I shall present here a slightly compacted ...
answered Jun 13 '16 at 5:57
32
Geodesics in Heat Algorithm
At the suggestion of @user21 I am splitting up my answers to help make the code(s) for calculating geodesics distances easier to find for other people interested in these sorts of algorithms.
The Geodesics in Heat algorithm is a fast approximate algorithm for estimating geodesic distances on discrete meshes (but also a variety ...
31
There is now a built-in version of an algorithm in v10.1: WordCloud. I wonder whether any of your nice algorithms introduced here had any influence on the built-in function...
Individual words can be styled, annotated, rotated, etc., so I must assume that there is a polygon-intersection checking algorithm running under the hood. Would be useful to know more ...
29
You'll be interested in the (undocumented!) functions Graphics`Mesh`IntersectQ[] (for checking the intersections) and Graphics`Mesh`FindIntersections[] (for actually finding them). As a sample:
BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
lins = Table[{Line[RandomReal[1, {2, 2}]]}, {42}];]
Graphics`Mesh`...
answered Oct 27 '12 at 0:47
28
Fixed (see below)
Here's an approach:
r1 = Exp[-x^3 - y] - 1 == z;
r2 = y == z;
We create ImplicitRegions:
reg1 = ImplicitRegion[r1, {x, y, z}];
reg2 = ImplicitRegion[r2, {x, y, z}];
The intersection of these regions is the line you seek:
reg = RegionIntersection[reg1, reg2];
And here is the length (note the inclusion of the range of values in ...
27
Short answer
Yes, it is possible to speed up the Delaunay-triangulation and make it as fast as it is in Matlab.
If you are not afraid of some setup-work, then one possibility is to use a package which calls a c-implementation of the Delaunay-triangulation. One package I know is qh-math which is available in the Wolfram-library:
This package includes ...
27
Here is a method that utilizes $H^1$-gradient flows. This is far quicker than the $L^2$-gradient flow (a.k.a. mean curvature flow) or using FindMinimum and friends, in particular when dealing with finely discretized surfaces.
Background
For those who are interested: A major reason for numerical slowness of $L^2$-gradient flow is the Courant–Friedrichs ...
26
The second "Neat Example" in the documentation for SmoothKernelDistribution contains this compiled function:
(* A region function for a bounding polygon using winding numbers: *)
inPolyQ =
Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}},
Block[{polySides = Length[polygon], X = polygon[[All, 1]],
Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = ...
26
Sometimes speed is an issue if there are many polygons and or many points to check. There is an excellent reference on this issue under http://erich.realtimerendering.com/ptinpoly/ with the main conclusion that the angle summation algorithm should be avoided if speed is the objective.
Below is my Mathematica implementation of the point in polygon problem ...
24
Here's a possible approach.
First use TetGen to tetrahedralize the data:
Needs["TetGenLink`"]
{pts, tetrahedra} = TetGenDelaunay[data3D];
Next define a function to compute the radius of the circumsphere of a tetrahedron (formula from Wikipedia)
csr[{aa_, bb_, cc_, dd_}] :=
With[{a = aa - dd, b = bb - dd, c = cc - dd},
Norm[a.a Cross[b, c] + b.b Cross[...
24
Another approach to this problem is computing the winding number by integrating $1/z$ centered on the point of interest along the contour of the polygon in the complex plane. Sure this isn't exactly efficient, but still I think it's nice to see this working in action. And since complex integration is feasible in Mathematica, I just tried :)
PointToComplex[{...
22
You could use this package to triangulate your polygon, and then use this barycentric formula on each of the triangles.
inside[{{x1_, y1_}, {x2_, y2_}, r3 : {x3_, y3_}}, r : {_, _}] :=
# >= 0 && #2 >= 0 && # + #2 < 1 & @@
LinearSolve[{{x1 - x3, x2 - x3}, {y1 - y3, y2 - y3}}, r - r3]
Example for a single triangle:
tri = {{...
22
Because (a) RegionPlot3D does not render edges well and (b) detailed information about the vertices and faces could be worthwhile, I will offer a solution that finds this information and displays it clearly. (The first two lines of code produce the region plot, if you just want to stop there; the rest develop the improved solution.)
I am stuck at one thing:...
22
Here is a solution I can think of. Idea is to take the FullPolygon of a given country and then triangulate the region. Once that is done take the underlying Graph and do a FindShortestPath. Result will not be too bad.
fullPoly = CountryData["Vietnam", "FullPolygon"];
pts = Flatten[fullPoly[[1, 1]], 1];
line = Polygon[Range[##]] & @@@ Partition[{1}~Join~
...
22
There are many ways to do this, modifying, improving my method or doing a completely different thing. My goal here is to show a very basic idea that should give you a start. LocatorPane and Manipulate give means of interactive addition/deletion and dragging of points in 2D plane. The problem is how to add an edge -- there has to be interaction between 2 ...
21
This is just a quick sketching out of an answer (rescales galore!)
textOnCurve[text_, f_, n_, p_: 0.01] :=
Text[Rotate[text, ArcTan @@ (f[Rescale[n + p, {0, 1}, {p, 1 - p}]] -
f[Rescale[n - p, {0, 1}, {p, 1 - p}]])], f[n]]
textCurve[string_, f_, stylef_: (# &), range_: {0, 1}] :=
With[{chars = Characters@...
21
Note that ListContourPlot3D takes the coordinates to be the position indices by default. If you want to keep the coordinates used in generating the data, then you have to include it.
data = Flatten[
Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]},
{x, -2, 2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];
plot = ListContourPlot3D[data, Contours -&...
21
[Edit: I found this method a rather pleasing application of analytic geometry, so I rewrote the explanation hopefully to do it justice. In fact, it can be applied in any dimension, and I've updated the code to be general]
Here's my way:
faces[simplex_] := Partition[simplex, Length@simplex - 1, 1, 1];
(* outward-oriented unit normals to each of faces[a,b,...
Only top voted, non community-wiki answers of a minimum length are eligible
