Ofsetting 3D polygon

[lvjonok]

I am interested whether Clipper2 provides or plans to provide tool to offset polygon in 3D. For example I have a polygon in 3D which belongs to cone surface. There are ways to project it into 2D and offset using this library. But transforming result back to 3D will not have the same offset because of distortion.

I have found such paper describing ways to solve this: here and wondered if something is already done in Clipper2?

[mmassing]

If I understand you correctly, you're looking for a way to offset using a position- and direction-dependent magnitude? I'm not sure if that would be easy to accomodate in the current offsetter. It might be easier to subdivide the polygon and project using a local equidistant projection, until a distortion criterion is met; then assemble the pieces in a intermediate, global projection (e.g. by using a small "epsilon" offset to "fudge" over map discontinuities, then apply a union).

[lvjonok]

Wow, thank you for idea! I will definitely think about it. Seems like I would need to subdivide polygon into overlapping portions, so after offset with negative value and union it would merge without error)

Anyway, I understand it is not trivial task, thank you for sharing)

[AngusJohnson]

Lev, I'm glad Manuel's post was helpful because 3D stuff is above my pay grade 😜.

[mmassing]

Hi @lvjonok , I just saw the variable-width offset discussion, maybe this could be helpful for your case as well: #499

[lvjonok]

Hi @lvjonok , I just saw the variable-width offset discussion, maybe this could be helpful for your case as well: #499

Thank you! I will check

[fellafoo]

The approach I took was to rotate the 3D poly points onto the xy plane, offset them, then rotate them back.

function IsAlmostEqual(A, B, Epsilon: Double): Boolean;
begin
  Result := Abs(A - B) < Epsilon;
end;

function RotatePoint(const Pt, Axis: TPoint3D; Angle: TFloat): TPoint3D;
var
  SinAngle, CosAngle, OneMinusCosAngle: Double;
begin
  if IsAlmostEqual(Angle, 0.0, Epsilon) then begin
    Result := Pt;
    Exit;
  end;

  SinAngle := Sin(Angle);
  CosAngle := Cos(Angle);
  OneMinusCosAngle := 1 - CosAngle;

  Result.x := (CosAngle + Axis.x * Axis.x * OneMinusCosAngle) * Pt.x + (Axis.x * Axis.y * OneMinusCosAngle - Axis.z * SinAngle) * Pt.y +
     (Axis.x * Axis.z * OneMinusCosAngle + Axis.y * SinAngle) * Pt.z;
  Result.y := (Axis.y * Axis.x * OneMinusCosAngle + Axis.z * SinAngle) * Pt.x + (CosAngle + Axis.y * Axis.y * OneMinusCosAngle) * Pt.y +
     (Axis.y * Axis.z * OneMinusCosAngle - Axis.x * SinAngle) * Pt.z;
  Result.z := (Axis.z * Axis.x * OneMinusCosAngle - Axis.y * SinAngle) * Pt.x + (Axis.z * Axis.y * OneMinusCosAngle + Axis.x * SinAngle) * Pt.y +
     (CosAngle + Axis.z * Axis.z * OneMinusCosAngle) * Pt.z;
end; { RotatePoint }

function RotatePointsToXYPlane(const Points: TMaxPntArr; nPnts: Integer; out RotAxis, Normal: TPoint3D; out Angle: TFloat): TMaxPntArr;
var
  i: Integer;
begin
  // Calculate the normal vector of the plane formed by the points
  Normal := CrossProduct(SubtractPoints(Points[2], Points[1]), SubtractPoints(Points[3], Points[1]));
  NormalizeVector(Normal);

  // Calculate the angle between the normal vector and the Z-axis
  if DotProduct(Normal, NewPoint(0.0, 0.0, 1.0)) < 0 then
       Angle := -ArcCos(DotProduct(Normal, NewPoint(0.0, 0.0, 1.0)))
  else
       Angle := ArcCos(DotProduct(Normal, NewPoint(0.0, 0.0, 1.0)));

  // Calculate the rotation axis - cross product of the normal vector and the Z-axis
  RotAxis := CrossProduct(Normal, NewPoint(0.0, 0.0, 1.0));
  NormalizeVector(RotAxis);

  // Rotate each point around the rotation axis by the calculated angle
  for i := 1 to nPnts do begin
    Result[i] := RotatePoint(Points[i], RotAxis, Angle);
  end;
end; { RotatePointsToXYPlane }

function RotatePointsBackToOriginalPlane(const Points: TMaxPntArr; nPnts: Integer; {const} RotAxis, Normal: TPoint3D; Angle: TFloat): TMaxPntArr;
var
  i: Integer;
begin
  // Rotate each point around the original rotation axis by the negated angle
  for i := 1 to nPnts do begin
    Result[i] := RotatePoint(Points[i], RotAxis, -Angle);
  end;
end; { RotatePointsBackToOriginalPlane }

I've not fully tested it but it has been working with the 3D polygons I've test thus far. ;-)

[lvjonok]

Good idea! If your polygons are on some plane in 3D it is good solution and it is straightforward.

My issue was mainly that for different height of original polygons that lie on non-linear surface (cone surface) I want different offset on XY plane. Here comes a lot of problems with projection of such surfaces into 2D.

But anyway, thank you for sharing!

[fellafoo]

Here's a better implementation (using a 4x4 rotation matrix) for flattening, offsetting, then restoring polygon points.

type
  TPoint3D = packed record
    case aSInt of
      0: (x, y, z: TFloat);
      1: (v: array [0 .. 2] of TFloat);
  end;

  TMaxPntArr = array [1 .. MaxPoly] of TPoint3D;
  Vec4 = packed array [1 .. 4] of TFloat;
  Mat4x4 = packed array [1 .. 4] of Vec4;

function TranslationMatrix(const Tx, Ty, Tz: TFloat): Mat4x4;
begin
  SetIdent(Result);
  Result[1, 4] := Tx;
  Result[2, 4] := Ty;
  Result[3, 4] := Tz;
end; { TranslationMatrix }

function MultiplyMatPoint(const Mat: Mat4x4; const Point: TPoint3D): TPoint3D;
var
  i: Integer;
  temp: array [1 .. 4] of TFloat;
begin
  for i := 1 to 4 do begin
    temp[i] :=
       Mat[i, 1] * Point.x +
       Mat[i, 2] * Point.y +
       Mat[i, 3] * Point.z +
       Mat[i, 4]; { Add the translation component for points }
  end;

  { Convert back to 3D point by dividing by the homogeneous coordinate (temp[4]) }
  Result.x := temp[1] / temp[4];
  Result.y := temp[2] / temp[4];
  Result.z := temp[3] / temp[4];
end; { Mat4x4MultPoint }

function SubtractPoints( { const } A, B: TPoint3D): TPoint3D;
begin
  Result.x := A.x - B.x;
  Result.y := A.y - B.y;
  Result.z := A.z - B.z;
end; { SubtractPoints }

function CrossProduct( { const } A, B: TPoint3D): TPoint3D; 
begin
  Result.x := A.y * B.z - A.z * B.y;
  Result.y := A.z * B.x - A.x * B.z;
  Result.z := A.x * B.y - A.y * B.x;
end; { CrossProduct - Multiply }

function SurfaceNormal(const p1, p2, p3: TPoint3D): TPoint3D;
var
  u, v: TPoint3D;
begin
  u := SubtractPoints(p2, p1);
  v := SubtractPoints(p3, p1);
  Result := CrossProduct(u, v);
end; { SurfaceNormal }

function NormalizeVectorFunc(const v: TPoint3D): TPoint3D;
var
  Mag: TFloat;
begin
  Mag := Magnitude(v);

  if (Mag <> 0) then begin
    Result.x := v.x / Mag;
    Result.y := v.y / Mag;
    Result.z := v.z / Mag;
  end
  else begin
    Result := v;
  end;
end; { NormalizeVectorFunc }

function RotationMatrixToXYPlane(const Normal: TPoint3D): Mat4x4;
var
  xAxis, yAxis, zAxis: TPoint3D;
begin
  zAxis := NormalizeVectorFunc(Normal);

  if (Abs(zAxis.x) <= Abs(zAxis.y)) and (Abs(zAxis.x) <= Abs(zAxis.z)) then
       xAxis := CrossProduct(NewPoint(0.0, 1.0, 0.0), zAxis)
  else
       xAxis := CrossProduct(NewPoint(0.0, 0.0, 1.0), zAxis);

  xAxis := NormalizeVectorFunc(xAxis);
  yAxis := CrossProduct(zAxis, xAxis);

  SetIdent(Result);
  Result[1, 1] := xAxis.x;
  Result[1, 2] := xAxis.y;
  Result[1, 3] := xAxis.z;

  Result[2, 1] := yAxis.x;
  Result[2, 2] := yAxis.y;
  Result[2, 3] := yAxis.z;

  Result[3, 1] := zAxis.x;
  Result[3, 2] := zAxis.y;
  Result[3, 3] := zAxis.z;
end; { RotationMatrixToXYPlane }

function RotatePoint(const Pt: TPoint3D; RotMat: Mat4x4): TPoint3D;
begin
  Result := MultiplyMatPoint(RotMat, Pt);
end; { RotatePoint }

function RotatePoints(Points: TMaxPntArr; nPnts: Integer; RotMat: Mat4x4): TMaxPntArr;
var
  i: Integer;
begin
  for i := 1 to nPnts do
       Result[i] := RotatePoint(Points[i], RotMat);
end; { RotatePoints }

function FlattenInflateRestore(Poly1: TMaxPntArr; var PntCnt1: Integer): TMaxPntArr;
var
  i: Integer;
  TransMat1, RotMat1: Mat4x4;
  TransPnts1, RotPoints1: TMaxPntArr;
  Normal1: TPoint3D;

begin
  { Move Poly1 points to origin }
  TransMat1 := TranslationMatrix(-Poly1[1].x, -Poly1[1].y, -Poly1[1].z);
  for i := 1 to PntCnt1 do
       TransPnts1[i] := MultiplyMatPoint(TransMat1, Poly1[i]);
  { Calculate normal for translated points }
  Normal1 := SurfaceNormal(TransPnts1);
  { Calculate rotation matrix from normal }
  RotMat1 := RotationMatrixToXYPlane(Normal1);
  { Rotate translated points to xy plane }
  RotPoints1 := RotatePoints(TransPnts1, PntCnt1, RotMat1);
  { Draw 'flattened' polygon }
  AddPolygon(RotPoints1, PntCnt1, 3);
  { Offset polygon points }
  Inflate(RotPoints1, PntCnt1);
  { Draw new offset points (point count may have increased) }
  for i := 1 to PntCnt1 do
       AddMarker(RotPoints1[i], 14);
  { Rotate points back from xy plane }
  RotPoints1 := RotatePoints(RotPoints1, PntCnt1, TransposeMatrix(RotMat1));
  { Move points back from origin }
  for i := 1 to PntCnt1 do
       RotPoints1[i] := MultiplyMatPoint(TranslationMatrix(Poly1[1].x, Poly1[1].y, Poly1[1].z), RotPoints1[i]);
  { Draw new offset polygon in same location as original points }
  AddPolygon(RotPoints1, PntCnt1, 3);
end; { FlattenInflateRestore }