Jump to content
A.M. Hoornweg

Recommended Posts

A.M. Hoornweg    126

A.M. Hoornweg
Posted

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.

 

Share this post

Link to post

ULIK    6

ULIK
Posted

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

Share this post

Link to post

A.M. Hoornweg    126

A.M. Hoornweg
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 ?

 

 

 

 

 

some_random_contour_map_.jpg

Edited by A.M. Hoornweg

Share this post

Link to post

A.M. Hoornweg    126

A.M. Hoornweg
Posted
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:

https://github.com/skia4delphi/skia4delphi/blob/eed4afbf8a34137a9bfa308bcb5ef87cee84abcb/Samples/Demo/FMX/Source/Sample.Form.PathsAndEffects.pas#L193

 

Result:

407181DE-5F0C-4E29-B067-A0BBFA6D6EBB.thumb.png.3b51235a5895e96571f33e9ec5603d9b.png

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

Share this post

Link to post

Sherlock    472

Sherlock
Posted

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.

Share this post

Link to post

vfbb    216

vfbb
Posted
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[0]);
    LPathBuilder.LineTo(APoints[1]);
    Exit(LPathBuilder.Detach);
  end;
  LSegments := Length(APoints) - 1;
  SetLength(LScratches, LSegments);
  LScratches[0].a := PointF(0, 0);
  LScratches[0].b := PointF(2, 2);
  LScratches[0].c := PointF(1, 1);
  LScratches[0].r := PointF(APoints[0].X + 2 * APoints[1].X, APoints[0].Y + 2 * APoints[1].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[0]);
  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:

image.thumb.png.8e591b43296bb2c8cd20f3fe88c773b2.png

 

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

  • Like 6

Share this post

Link to post

Anders Melander    1129

Anders Melander
Posted
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.

https://towardsdatascience.com/numerical-interpolation-natural-cubic-spline-52c1157b98ac

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 😕

Share this post

Link to post

Anders Melander    1129

Anders Melander
Posted

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).

  • Like 2

Share this post

Link to post

angusj    90

angusj
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[0];
  for i := 0 to len -1 do
  begin
    if i = len -1 then
    begin
      vec := GetAvgUnitVector(unitVecs[i], unitVecs[0]);
      d2 := pl[0]*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[1] := Result[0];
    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;

 

test20.thumb.png.d97635270def988873c388657a5bd5ff.png

  

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;

test80.thumb.png.f2e5c34379525f78150f0aef978729bb.png

 

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
  • Like 3

Share this post

Link to post

A.M. Hoornweg    126

A.M. Hoornweg
Posted
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!

  • Like 1

Share this post

Link to post

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Go To Topic Listing Algorithms, Data Structures and Class Design
×