# Algorithm to find out points of inflection for a polyline

I am trying to find out the points of inflection , i.e. points where the curves in a line start and end. If you look at the image, the green line may be a road or a stream, and the black points are the points where the curves start and end.

What would be the high level steps to automate the generation of these points? I have ArcGIS desktop and am quite handy with ArcObjects.

• The source data is a polyline made of line segments and you want to aproximate it with curves, or does it already has arc segments? Oct 18, 2012 at 12:47
• Currently, it is made of line segments. Oct 18, 2012 at 12:53
• The illustration in this question looks remarkably like one published at esri.com/news/arcuser/0110/turning.html. Sep 23, 2013 at 18:05
• @whuber: Very astute observation. That was exactly the data source I had used to create the image. Sep 24, 2013 at 3:33

When the curve is comprised of line segments, then all interior points of those segments are inflection points, which is not interesting. Instead, the curve should be thought of as being approximated by the vertices of those segments. By splining a piecewise twice-differentiable curve through those segments, we can then compute the curvature. An inflection point, strictly speaking, is then a place where the curvature is zero.

In the example there are longish stretches where the curvature is nearly zero. This suggests that the indicated points ought to approximate the ends of such stretches of low-curvature regions.

An effective algorithm will therefore spline the vertices, compute the curvature along a dense set of intermediate points, identify ranges of near-zero curvature (using some reasonable estimate of what it means to be "near"), and mark the endpoints of those ranges.

Here is working `R` code to illustrate these ideas. Let's begin with a line string expressed as a sequence of coordinates:

``````xy <- matrix(c(5,20, 3,18, 2,19, 1.5,16, 5.5,9, 4.5,8, 3.5,12, 2.5,11, 3.5,3,
2,3, 2,6, 0,6, 2.5,-4, 4,-5, 6.5,-2, 7.5,-2.5, 7.7,-3.5, 6.5,-8), ncol=2, byrow=TRUE)
``````

Spline the x and y coordinates separately to achieve a parametrization of the curve. (The parameter will be called `time`.)

``````n <- dim(xy)[1]
fx <- splinefun(1:n, xy[,1], method="natural")
fy <- splinefun(1:n, xy[,2], method="natural")
``````

Interpolate the splines for plotting and computation:

``````time <- seq(1,n,length.out=511)
uv <- sapply(time, function(t) c(fx(t), fy(t)))
``````

We need a function to compute the curvature of a parameterized curve. It needs to estimate the first and second derivatives of the spline. With many splines (such as cubic splines) this is an easy algebraic calculation. `R` provides the first three derivatives automatically. (In other environments, one might want to compute the derivatives numerically.)

``````curvature <- function(t, fx, fy) {
# t is an argument to spline functions fx and fy.
xp <- fx(t,1); yp <- fy(t,1)            # First derivatives
xpp <- fx(t,2); ypp <- fy(t,2)          # Second derivatives
v <- sqrt(xp^2 + yp^2)                  # Speed
(xp*ypp - yp*xpp) / v^3                 # (Signed) curvature
# (Left turns have positive curvature; right turns, negative.)
}

kappa <- abs(curvature(time, fx, fy))     # Absolute curvature of the data
``````

I propose to estimate a threshold for zero curvature in terms of the extent of the curve. This at least is a good starting point; it ought to be adjusted according to the tortuosity of the curve (that is, increased for longer curves). This will later be used for coloring the plots according to curvature.

``````curvature.zero <- 2*pi / max(range(xy[,1]), range(xy[,2])) # A small threshold
i.col <- 1 + floor(127 * curvature.zero/(curvature.zero + kappa))
palette(terrain.colors(max(i.col)))                        # Colors
``````

Now that the vertices have been splined and the curvature computed, it remains only to find the inflection points. To show them we may plot the vertices, plot the spline, and mark the inflection points on it.

``````plot(xy, asp=1, xlab="x",ylab="y", type="n")
tmp <- sapply(2:length(kappa), function(i) lines(rbind(uv[,i-1],uv[,i]), lwd=2, col=i.col[i]))
points(t(sapply(time[diff(kappa < curvature.zero/2) != 0],
function(t) c(fx(t), fy(t)))), pch=19, col="Black")
points(xy)
``````

The open points are the original vertices in `xy` and the black points are the inflection points automatically identified with this algorithm. Because curvature cannot reliably be computed at the endpoints of the curve, those points are not specially marked.

• Maybe the terminology I used was wrong. What you assumed is exactly what I wanted. Your answer looks promising, and I will have to work with R to process my Shapefile. Oct 18, 2012 at 16:45

You can use the Densify tool. For this case, you choose to densify by angle, Next, choose the the maximum angle accepted in a straight line. Then apply to result line to the tool Split line at vertices. Finally, delete the lines having shape_length smaller the minimum road length.

In this picture, we see three steps:

1- Densify the line using angle. I have used 10 degrees as the parameter, and we used splitline. In the picture, the curved line is in its initial phase.

``````arcpy.Densify_edit("line" , "ANGLE" , "","",10)
arcpy.SplitLine_management("line" , "line_split")
``````

2- Select the segments where having shape_length not redundant. As we see in the table, I have not selected those redundant lengths. Then, I select them into a new feature class.

``````arcpy.Select_analysis("line_split" , "line_split_selected")
``````

3- We have extracted the vertices located in the edges of the lines, which are inflection points.

``````arcpy.FeatureVerticesToPoints_management("line_split_selected" , "line_split_pnt" , "DANGLE")
``````
• I have the same comments and questions concerning your other answer: it's a nice idea, but at the same time it's not clear that it will produce the desired result, nor how one should choose the threshold angle. Could you provide an illustration of the output so readers can evaluate what this proposal really does? Providing worked examples is especially important when recommending ESRI software as part of a solution, because their algorithms usually are not documented, making it impossible to know exactly what they are doing. Oct 18, 2012 at 17:15
• to be very sure that is a working solution , i need to test it, but i can't test it, i'm missing the data,so i suppose that the proposed tools by ESRI gonna work like expected, but this answers need to be tested further. Oct 18, 2012 at 17:25
• we could name them ideas and not answers Oct 18, 2012 at 17:25
• Would you like me to move them into comments, then? BTW, if you want data to test, you could--for a start--use the coordinates I posted in my answer, because they're close to the illustration in the question. But why not just use any geographic data you have handy? Oct 18, 2012 at 18:33
• yes really this solution is working better at extracting just straight lines. Oct 18, 2012 at 20:36

You can use the Generalize tool which has the maximum offset from the original line as a parameter, so you can choose the offset which fits your case.

If we name the original line "line_cur", and the generalized one "line_gen", then we could clip "line_cur" by "line_gen". The result will be the straight segment of "line_cur". Then we could clean some very short segment by deleting them with an sql query which selects the Shape_length greater that minimum road length.

• This is a nice idea. It's not clear how well it would work in practice, though. Could you perhaps show an example showing the inflection points that are found? Oct 18, 2012 at 16:46
• i have made an edit to include a picture, the picture explains how this tool could make a line sticking with straight segments, so we have to make a clip to the old line, to extract just the old straight lines segments Oct 18, 2012 at 16:50
• is anything not clear i'm available to respond your questions ? Oct 18, 2012 at 16:51
• I don't see any inflection points identified in the illustration. Where would they be, exactly? And how should one choose the tolerance for generalization? Oct 18, 2012 at 16:53
• i need to some data to perform the test, but i think we should choose the tolerance by experimentation Oct 18, 2012 at 16:55