Jump to content

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.

 

Share this post


Link to post

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

Share this post


Link to post

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

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
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 10

Share this post


Link to post
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

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

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 4

Share this post


Link to post
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

Hi Angus, I just discovered your GetSmoothPath() routine. I'm very interested to use it - but I can't make it work.

 

Contrary to what you write the source for the function MakePath is not included.

 

From what I see on the Image32 pages (http://www.angusj.com/delphi/image32/Docs/Units/Img32.Vector/Routines/MakePath.htm), this again calls MakePathI and CSpline which I can't find anywhere. I also wonder why do you convert an array of Doubles to an array of integers and then assign the result again to an array of doubles?

 

How should I use the the tGPGraphicsPath.Flatten of GDI+? It doesn't accept a tPathD as an argument.

 

Could you please elaborate?

Share this post


Link to post
12 hours ago, xstrider said:

Hi Angus, I just discovered your GetSmoothPath() routine. I'm very interested to use it - but I can't make it work.

Hi xstrider.

I've actually modified this function again since posting the link above, but I didn't want to abuse this thread since it's not specifically about my routine.

Anyhow, here's the link to the up to date documentation on this function that's been renamed SmoothToCubicBezier (which I think better describes what the function does).

This can be found in the Img32.Extra unit that part of my Image32 Library on GitHub.

Both MakePath and FlattenCBezier are found in Img32.Vector.

 

12 hours ago, xstrider said:

I also wonder why do you convert an array of Doubles to an array of integers and then assign the result again to an array of doubles?

I'm not sure what you're referring to here. The MakePath function parameter is array of double, and the entire Image32 library uses double coordinate values.

 

12 hours ago, xstrider said:

How should I use the the tGPGraphicsPath.Flatten of GDI+? It doesn't accept a tPathD as an argument.

It appears you're using GDIPlus and I'm not familiar with it. However, while GraphicsPath does accept cubic bezier input (and will flatten these), these beziers appear to be cubic bezier splines, which are very different to cubic bezier arrays. Nevertheless you could copy FlattenCBezier from my library and simply add the flattened path using GraphicsPath's AddPolygon method.

 

I hope that helps.

 

ps: Here's a link to a short video that demonstrates this smoothing (together with vectorizing monochrome raster images, and path simplification).

Edited by angusj
  • Like 1

Share this post


Link to post
8 hours ago, angusj said:

Hi xstrider.

I've actually modified this function again since posting the link above, but I didn't want to abuse this thread since it's not specifically about my routine.

Anyhow, here's the link to the up to date documentation on this function that's been renamed SmoothToCubicBezier (which I think better describes what the function does).

This can be found in the Img32.Extra unit that part of my Image32 Library on GitHub.

Both MakePath and FlattenCBezier are found in Img32.Vector.

 

I'm not sure what you're referring to here. The MakePath function parameter is array of double, and the entire Image32 library uses double coordinate values.

 

It appears you're using GDIPlus and I'm not familiar with it. However, while GraphicsPath does accept cubic bezier input (and will flatten these), these beziers appear to be cubic bezier splines, which are very different to cubic bezier arrays. Nevertheless you could copy FlattenCBezier from my library and simply add the flattened path using GraphicsPath's AddPolygon method.

 

I hope that helps.

 

ps: Here's a link to a short video that demonstrates this smoothing (together with vectorizing monochrome raster images, and path simplification).

 

Hi Angus,

 

thanks for your time and I have great respect for expertise which is way above my head.

 

But your earlier post promised a stand alone solution ("It requires no specific graphics library to use, just a few extra functions (also included below)") but this doesn't seem to be true. Well, it was worth a try.

 

 

Share this post


Link to post
6 minutes ago, xstrider said:

 

But your earlier post promised a stand alone solution ("It requires no specific graphics library to use, just a few extra functions (also included below)") but this doesn't seem to be true.

I wasn't suggesting that you download Image32, just a couple of extra functions. And MakePath was only a shortcut to demonstrate path construction, which of course you'd generate in your own way. And even with FlattenCBezier, I was (I think quite reasonably) expecting that your graphics library could handle this since cubic bezier paths are a very common graphics data structure. Anyhow, I'm sorry I dissapointed you.

Share this post


Link to post

I particularly like your Iucundae_21_B. Wonderful colours.

Edited by angusj

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

×