## Recommended Posts

Hello all,

does anyone know a nice algorithm that will connect the points in a polygon in a rounded fashion? Its intended use is a contour map.

The polygon must go through all vertices but I'd like to avoid sharp corners for a more natural look;

it's OK if the lines between the vertices are slightly curved. The few solutions I found on the internet tend to bend without touching the vertices.

Have you already tried it with a combination of splines? They would touch the vertices.

From the Bézier curve you'll eventually have to move on to (Bézier) splines to get from a polygon to a bezigon. And not sure about VCL, but FMX has some Bézier procedures.

Edith says: Naturally VCL has it, because MS has it: https://docwiki.embarcadero.com/Libraries/Sydney/en/Vcl.Graphics.TCanvas.PolyBezier

Posted (edited)

It's not trivial at all.

I've downloaded a random contour map from the internet and attach it here. This is a map of a terrain, consisting of polygons that describe points of equal altitude above sea level.

So there's a polygon describing an altitude of 350 metres, another one that corresponds to 360 metres, another one for 370 metres etcetera.

The problem:  The polygons may never intersect (any point of the topography has only one altitude).  If an automated routine is to be used to create smooth connections between points of a polygon,  it is necessary to avoid "wild deflections" that would make these lines cross adjacent polygons.  A human being would have no problem drawing such connections by hand, intuitively. But an algorithm ?

Edited by A.M. Hoornweg

3 minutes ago, vfbb said:

If you don't mind third-party solutions, Skia has an effect that applies rounded edges automatically and works for Console, Vcl and FMX:

Result:

Such rounded edges "miss" the points of the polygon, they deviate before hitting the point.

Ok, so the points need to be hit, but the lines in between may be rounded. That might not be covered by splines to well. 😞

Well strike that, of course the Points will be hit. It is "just" a matter of choosing the controlpoints well to get a decent shape. But therein lies the rub I guess.

2 hours ago, A.M. Hoornweg said:

Such rounded edges "miss" the points of the polygon, they deviate before hitting the point.

Here's an example:

```uses
Skia;

function MakeCubicSplineInterpolation(const APoints: TArray<TPointF>): ISkPath;
var
LPathBuilder: ISkPathBuilder;
LSegments: Integer;
I: Integer;
mx: Single;
my: Single;
LScratches: array of
record
a, b, c, r, p: TPointF;
end;
begin
LPathBuilder := TSkPathBuilder.Create;
if Length(APoints) < 2 then
Exit(LPathBuilder.Detach);
if Length(APoints) = 2 then
begin
LPathBuilder.MoveTo(APoints);
LPathBuilder.LineTo(APoints);
Exit(LPathBuilder.Detach);
end;
LSegments := Length(APoints) - 1;
SetLength(LScratches, LSegments);
LScratches.a := PointF(0, 0);
LScratches.b := PointF(2, 2);
LScratches.c := PointF(1, 1);
LScratches.r := PointF(APoints.X + 2 * APoints.X, APoints.Y + 2 * APoints.Y);
for I := 1 to LSegments - 2 do
begin
LScratches[I].a := PointF(1, 1);
LScratches[I].b := PointF(4, 4);
LScratches[I].c := PointF(1, 1);
LScratches[I].r := PointF(4 * APoints[i].X + 2 * APoints[I + 1].X, 4 * APoints[I].Y + 2 * APoints[I + 1].Y);
end;
LScratches[LSegments - 1].a := PointF(2, 2);
LScratches[LSegments - 1].b := PointF(7, 7);
LScratches[LSegments - 1].c := PointF(0, 0);
LScratches[LSegments - 1].r := PointF(8 * APoints[LSegments - 1].X + APoints[LSegments].X, 8 * APoints[LSegments - 1].Y + APoints[LSegments].Y);
for I := 1 to LSegments - 1 do
begin
mx := LScratches[I].a.X / LScratches[I - 1].b.X;
my := LScratches[I].a.Y / LScratches[I - 1].b.Y;
LScratches[I].b := LScratches[I].b - PointF(mx * LScratches[I - 1].c.X, my * LScratches[I - 1].c.Y);
LScratches[I].r := LScratches[I].r - PointF(mx * LScratches[I - 1].r.X, my * LScratches[I - 1].r.Y);
end;
LScratches[LSegments - 1].p := PointF(LScratches[LSegments - 1].r.X / LScratches[LSegments - 1].b.X,
LScratches[LSegments - 1].r.Y / LScratches[LSegments - 1].b.Y);
for I := Length(APoints) - 3 downto 0 do
begin
LScratches[I].p := PointF((LScratches[I].r.X - LScratches[I].c.X * LScratches[I + 1].p.X) / LScratches[I].b.X,
(LScratches[I].r.Y - LScratches[I].c.Y * LScratches[I + 1].p.Y) / LScratches[I].b.Y);
end;
LPathBuilder.MoveTo(APoints);
for I := 0 to LSegments - 2 do
begin
LPathBuilder.CubicTo(LScratches[I].p,
PointF(2 * APoints[I + 1].X - LScratches[I + 1].p.X, 2 * APoints[I + 1].Y - LScratches[I + 1].p.Y), APoints[I + 1]);
end;
LPathBuilder.CubicTo(LScratches[LSegments - 1].p,
PointF(0.5 * (APoints[LSegments].X + LScratches[LSegments - 1].p.X),
0.5 * (APoints[LSegments].Y + LScratches[LSegments - 1].p.Y)), APoints[LSegments]);
Result := LPathBuilder.Detach;
end;

procedure TForm1.SkPaintBox1Draw(ASender: TObject; const ACanvas: ISkCanvas;
const ADest: TRectF; const AOpacity: Single);
var
LPaint: ISkPaint;
LMyPoints: TArray<TPointF>;
begin
LMyPoints := [PointF(62, 511), PointF(162, 605), PointF(262, 610),
PointF(362, 402), PointF(462, 959), PointF(562, 58), PointF(662, 272),
PointF(762, 99), PointF(862, 759), PointF(962, 945)];

LPaint := TSkPaint.Create(TSkPaintStyle.Stroke);
LPaint.Color := TAlphaColors.Red;
LPaint.AntiAlias := True;
LPaint.StrokeWidth := 3;
LPaint.StrokeCap := TSkStrokeCap.Round;
ACanvas.DrawPath(MakeCubicSplineInterpolation(LMyPoints), LPaint);
LPaint.StrokeWidth := 10;
LPaint.Color := TAlphaColors.Black;
ACanvas.DrawPoints(TSkDrawPointsMode.Points, LMyPoints, LPaint);
end;```

Result:

Note: You don't need to use Skia, it was just a facilitator for the example.

• 6

@vfbb: thanks for the example, I'll look into it!

11 hours ago, Rollo62 said:

So you need to fit all datapoints, but want to avoid large overshootings, then this paper is maybe also interesting for you.

Quote

We will use the top-down approach and make sure you visualize while you’re reading to understand it better.

I counted one photo of a Chinese paper lamp, two irrelevant meme pics, two general conceptual illustrations and 10 pages of math with no illustrations 😕

I just remembered that Graphics32 has two examples which demonstrates interpolation:

Both of these just uses Graphics32 for output. The curve generation is independent.

Built into Graphics32 there's also the TCanvas32.CurveTo method which does cubic Bézier interpolation (4 control points) and the TCanvas32.ConicTo method which does quadratic Bézier interpolation (3 control points).

• 2

Posted (edited)

Here's my GetSmoothPath() routine. It requires no specific graphics library to use, just a few extra functions (also included below).

This function generates an array of control points that's very easily converted into a flattened cubic bezier path using just about any 2D graphics library.

(nb: The code below has been written with simplicity as the focus rather than performance.)

```uses
SysUtils, Math;

type

TPointD = record
X, Y: double;
end;

TPathD = array of TPointD;
TArrayOfDouble = array of double;

function DistanceSqrd(const pt1, pt2: TPointD): double;
begin
result := Sqr(pt1.X - pt2.X) + Sqr(pt1.Y - pt2.Y);
end;

function Distance(const pt1, pt2: TPointD): double;
begin
Result := Sqrt(DistanceSqrd(pt1, pt2));
end;

function OffsetPoint(const pt: TPointD; dx, dy: double): TPointD;
begin
result.x := pt.x + dx;
result.y := pt.y + dy;
end;

function GetAvgUnitVector(const vec1, vec2: TPointD): TPointD;
var
inverseHypot: Double;
begin
Result.X := (vec1.X + vec2.X) * 0.5;
Result.y := (vec1.Y + vec2.Y) * 0.5;
inverseHypot := 1 / Hypot(Result.X, Result.Y);
Result.X := Result.X * inverseHypot;
Result.Y := Result.Y * inverseHypot;
end;

procedure MakeSymmetric(var val1, val2: double);
begin
val1 := (val1 + val2) * 0.5;
val2 := val1;
end;

function GetUnitVector(const pt1, pt2: TPointD): TPointD;
var
dx, dy, inverseHypot: Double;
begin
if (pt1.x = pt2.x) and (pt1.y = pt2.y) then
begin
Result.X := 0;
Result.Y := 0;
Exit;
end;
dx := (pt2.X - pt1.X);
dy := (pt2.Y - pt1.Y);
inverseHypot := 1 / Hypot(dx, dy);
dx := dx * inverseHypot;
dy := dy * inverseHypot;
Result.X := dx;
Result.Y := dy;
end;

// GetSmoothPath - returns cubic bezier control points
//   parameters: 1. path for smoothing
//               2. whether or not the smoothed path will closed
//               3. percent smoothness (0..100)
//               4. maximum dist control pts from path pts (0 = no limit)
//               5. symmetric vs asymmmetric control pts
function GetSmoothPath(const path: TPathD; pathIsClosed: Boolean;
percentOffset, maxCtrlOffset: double;
symmetric: Boolean): TPathD;
var
i, len, prev: integer;
vec: TPointD;
pl: TArrayOfDouble;
unitVecs: TPathD;
d, d1,d2: double;
begin
Result := nil;
len := Length(path);
if len < 3 then Exit;
d :=  Max(0, Min(100, percentOffset))/200;
if maxCtrlOffset <= 0 then maxCtrlOffset := MaxDouble;

SetLength(Result, len *3 +1);
prev := len-1;
SetLength(pl, len);
SetLength(unitVecs, len);
for i := 0 to len -1 do
begin
pl[i] := Distance(path[prev], path[i]);
unitVecs[i] := GetUnitVector(path[prev], path[i]);
prev := i;
end;

Result[len*3] := path;
for i := 0 to len -1 do
begin
if i = len -1 then
begin
vec := GetAvgUnitVector(unitVecs[i], unitVecs);
d2 := pl*d;
end else
begin
vec := GetAvgUnitVector(unitVecs[i], unitVecs[i+1]);
d2 := pl[i+1]*d;
end;
d1 := pl[i]*d;
if symmetric then MakeSymmetric(d1, d2);
if i = 0 then
Result[len*3-1] := OffsetPoint(path[i],
-vec.X * Min(maxCtrlOffset, d1), -vec.Y * Min(maxCtrlOffset, d1))
else
Result[i*3-1] := OffsetPoint(path[i],
-vec.X * Min(maxCtrlOffset, d1), -vec.Y * Min(maxCtrlOffset, d1));
Result[i*3] := path[i];
Result[i*3+1] := OffsetPoint(path[i],
vec.X * Min(maxCtrlOffset, d2), vec.Y * Min(maxCtrlOffset, d2));
end;
if not pathIsClosed then
begin
Result := Result;
dec(len);
Result[len*3-1] := Result[len*3];
SetLength(Result, Len*3 +1);
end;
end;```

And here's what it produces ...

the path to smooth (black),

the cubic bezier control path produced by GetSmoothPath() (blue)

and the flattened cubic bezier path (2D graphics library of you choice required) (red).

```var
TPathD path;
begin
path := MakePath([190,120, 260,270, 560,120, 190,490]);
path := GetSmoothPath(path, true, 20, 0, false);
path := ThirdParty2DGraphicsLibrary.FlattenCBezier(path);
end;```

```var
TPathD path;
begin
path := MakePath([190,120, 260,270, 560,120, 190,490]);
path := GetSmoothPath(path, true, 80, 0, false);
path := ThirdParty2DGraphicsLibrary.FlattenCBezier(path);
end;```

On 5/20/2022 at 8:45 PM, A.M. Hoornweg said:

The problem:  The polygons may never intersect (﻿﻿﻿﻿﻿

Edit: The best way to avoid intersections is to make sure you have enough data points before generating your curves.

Edited by angusj
Note to explain the code focuses on simplicity not performance
• 2

10 hours ago, Anders Melander said:

I just remembered that Graphics32 has two examples which demonstrates interpolation:

Built into Graphics32 there's also the TCanvas32.CurveTo method which does cubic Bézier interpolation (4 control points) and the TCanvas32.ConicTo method which does quadratic Bézier interpolation (3 control points).

Very interesting, thanks for the link!

• 1