expertcad/POWERCAD30/UNITS/rrEllipses.pas

{
  rr.inc:	Include file, defines conditional compilation symbols
		for Borland Delphi / C++Builder RAD tools

  Author:       Robert Roßmair, Robert.Rossmair@t-online.de

  Version:      1.5, 24-Mar-2001

  Predefined compiler version symbols are:
  (VERxyz means compiler version xy.z)

  VER80		Delphi 1.0
  VER90		Delphi 2.0
  VER93		C++Builder 1.0
  VER100	Delphi 3.0
  VER110	C++Builder 3.0
  VER120	Delphi 4.0
  VER125	C++Builder 4.0
  VER130	Delphi 5.0
  VER135	C++Builder 5.0
  VER140        Kylix 1.0
}

{$B-} { Complete Boolean Evaluation }
{$R-} { Range-Checking }
{$V-} { Var-String Checking }
{$T-} { Typed @ operator }
{$X+} { Extended syntax }
{$P+} { Open string params }
{$IFDEF WIN32}
{$J+} { Writeable structured consts }
{$H+} { Use long strings by default }
{$ENDIF}

{$DEFINE VER90_up}
{$DEFINE VER93_up}
{$DEFINE VER100_up}
{$DEFINE VER110_up}
{$DEFINE VER120_up}
{$DEFINE VER125_up}
{$DEFINE VER130_up}
{$DEFINE VER135_up}
{$DEFINE VER140_up}

{$IFNDEF VER140}
  {$UNDEF VER140_up}
  {$IFNDEF VER135}
    {$UNDEF VER135_up}
    {$IFNDEF VER130}
      {$UNDEF VER130_up}
      {$IFNDEF VER125}
        {$UNDEF VER125_up}
        {$IFNDEF VER120}
          {$UNDEF VER120_up}
          {$IFNDEF VER110}
            {$UNDEF VER110_up}
            {$IFNDEF VER100}
              {$UNDEF VER100_up}
              {$IFNDEF VER93}
                {$UNDEF VER93_up}
                {$IFNDEF VER90}
                  {$UNDEF VER90_up}
                {$ENDIF}
              {$ENDIF}
            {$ENDIF}
          {$ENDIF}
        {$ENDIF}
      {$ENDIF}
    {$ENDIF}
  {$ENDIF}
{$ENDIF}


{$IFDEF VER93}  { Borland C++Builder 1.0 }
  {$DEFINE CBUILDER}
{$ENDIF}

{$IFDEF VER110} { Borland C++Builder 3.0 }
  {$DEFINE CBUILDER}
{$ENDIF}

{$IFDEF VER125} { Borland C++Builder 4.0 }
  {$DEFINE CBUILDER}
{$ENDIF}

{$IFDEF VER135} { Borland C++Builder 5.0 }
  {$DEFINE CBUILDER}
{$ENDIF}

{$IFDEF VER140} { Borland Kylix 1.0 }
  {$DEFINE Kylix}
{$ENDIF}

{$IFDEF CBUILDER}
  {$IFNDEF VER93}
  {$ObjExportAll On}
  {$ENDIF}
  {$DEFINE NoMath}      { No Math unit }
{$ENDIF}

{$IFDEF VER100_up}		{ Delphi 3 and above }
  {$DEFINE RR_PackageSupport}
  {$DEFINE RR_Interfaces}
  {$DEFINE RR_Resourcestring}
  {$DEFINE RR_EnhancedTBitmap}
  {$DEFINE RR_JPEG}
  {$IFDEF VER110_up}		{ C++Builder 3 and above }
    {$DEFINE RR_DefaultParams}
    {$IFDEF VER120_up}		{ Delphi 4 and above }
      {$DEFINE RR_Overloading}
      {$DEFINE RR_Int64}
      {$IFDEF VER140_up}
        {$DEFINE RR_CLXAvailable}
      {$ENDIF}
    {$ENDIF}
  {$ENDIF}
{$ENDIF}

{$DEFINE RR_PicEditBroken}    { PicEdit unit broken (missing LibConst) }

{$IFNDEF VER100 (not Delphi 3) }
  {$IFNDEF VER110 (not C++Builder 3) }
     {$IFNDEF VER120 (not Delphi 4) }
     {$UNDEF RR_PicEditBroken}
     {$ENDIF}
  {$ENDIF}
{$ENDIF}

{$IFDEF VER100}
  { Graphics.CopyBitmap bug: invalid parameters are passed to
    CreateDIBSection when HandleType = bmDIB and the value of the
    boolean expression (Width*Height=0) changes.
    This will cause an EOutOfResources exception. }
  {$DEFINE RR_CopyBitmapBug}
{$ENDIF}

{$IFDEF Win32}
//  {$DEFINE MSWindows}
{$ENDIF}


(***************************************************************************

     The contents of this file are subject to the Mozilla Public License
     Version 1.0 (the "License"); you may not use this file except in
     compliance with the License. You may obtain a copy of the License at
     http://www.mozilla.org/MPL/

     Software distributed under the License is distributed on an "AS IS"
     basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
     License for the specific language governing rights and limitations
     under the License.

     The Original Code is rrEllipse.pas.

     The Initial Developer of the Original Code is Robert Rossmair.

	mailto:Robert.Rossmair@t-online.de
	http://home.t-online.de/home/Robert.Rossmair/

     Portions created by Robert Rossmair are Copyright © 1999, 2000
     Robert Rossmair. All Rights Reserved.

     Contributor(s): ______________________________________.

     Compatibility: Delphi 4(?), Delphi 5, Kylix 1

     Last Modified: March 24, 2001

 ***************************************************************************)

{.$DEFINE Debug}
{$IFDEF Win32}
  {$DEFINE MSWindows}
{$ENDIF}

unit rrEllipses;

interface

uses
{$IFDEF MSWindows}
  Windows, Graphics,
{$ELSE}
  Types, Qt, QGraphics,
{$ENDIF}
  Math,
  SysUtils;

type
  T2DVector = packed record
    X, Y: Double;
  end;
  T2DPoint = T2DVector;

  T2DPointArray = array of T2DPoint;

  TCubicBezierCurve = array[0..3] of T2DPoint;
  PCubicBezierCurve = ^TCubicBezierCurve;

{$IFDEF MSWindows}
procedure CircularArc(DC: hDC; X, Y, Radius, BegAngle, EndAngle: Double);
{$ENDIF}
procedure DrawCircularArc(Canvas: TCanvas; X, Y, Radius, BegAngle, EndAngle: Double);

{ Draw elliptical arc rotated by angle Angle;

  X, Y: Center coordinates
  A, B: Half axis lengths

  BegAngle, EndAngle are relative to Angle }

{$IFDEF MSWindows}
procedure EllipticalArc(DC: hDC; X, Y, A, B, Angle, BegAngle, EndAngle: Double);
{$ENDIF}
procedure DrawEllipticalArc(Canvas: TCanvas; X, Y, A, B, Angle, BegAngle, EndAngle: Double);

{ Ellipse approximation by 8 cubic Bezier curves.

  Relative approximation error < 4.25E-6
  X, Y: Center coordinates
  A, B: Half axis lengths
  Angle: Angle in radians }

{$IFDEF MSWindows}
procedure EllipseEx(DC: hDC; X, Y, A, B, Angle: Double);
  {$IFDEF RR_Overloading} overload; {$ENDIF}
{$ENDIF}
procedure DrawEllipse(Canvas: TCanvas; X, Y, A, B, Angle: Double);
  {$IFDEF RR_Overloading} overload; {$ENDIF}

{ Ellipse approximation by N cubic Bezier curves, N >= 2

  Relative approximation error <= (1-cos(Pi/N))³/(54*(1+cos(Pi/N)))

  X, Y: Center coordinates
  A, B: Half axis lengths
  Angle: Angle in radians }

{$IFDEF RR_Overloading}
  {$IFDEF MSWindows}
  procedure EllipseEx(DC: hDC; X, Y, A, B, Angle: Double; N: Integer); overload;
  {$ENDIF}
procedure DrawEllipse(Canvas: TCanvas; X, Y, A, B, Angle: Double; N: Integer); overload;
{$ENDIF RR_Overloading}

{ Helper functions }

function FPoint(X, Y: Double): T2DPoint;
function NormAngle(const Angle: Double): Double; // fold angle to [-Pi, Pi[
procedure Rotate(var P: array of T2DVector; Angle: Double);
procedure Scale(var P: array of T2DVector;
  ScaleFactorX, ScaleFactorY: Double);
procedure Translate(var P: array of T2DVector; dX, dY: Double);

procedure BezierArcPoints(
  var Points: T2DPointArray;
  X, Y,
  Radius,
  BegAngle,
  EndAngle: Double);

procedure BezierElpArcPoints(
  var Points: T2DPointArray;
  X, Y,
  aLen,bLen,
  angle,
  BegAngle,
  EndAngle: Double);

procedure GetCircleThruPoints(
  const P1, P2, P3: T2DPoint;
  var Center: T2DPoint;
  var Radius: Double);
{$IFDEF RR_Overloading} overload;

procedure GetCircleThruPoints(
  const X1, Y1, X2, Y2, X3, Y3: Double;
  var CenterX, CenterY: Double;
  var Radius: Double); overload;

procedure GetCircleThruPoints(
  const P1, P2, P3: TPoint;
  var Center: T2DPoint;
  var Radius: Double); overload;

{$ENDIF RR_Overloading}
var
  DxfMode: Boolean = False;
  DxfMode16: Boolean = False;
  DxfMode32: Boolean = False;

implementation
//uses DrawObjects;
{$IFDEF RR_Resourcestring}
resourcestring
{$ELSE}
const
{$ENDIF}
  SPointsOnStraightLine = 'GetCircleThruPoints: Points lie on straight line';

const
  Pi2 = 2*Pi;
  Inv2Pi = 0.5/Pi;
  Epsilon = 1E-12;

//var
//  DxfMode: Boolean = False;

{$IFDEF RR_Overloading}
{$IFDEF MSWindows}
procedure DrawEllipse(Canvas: TCanvas; X, Y, A, B, Angle: Double; N: Integer);
begin
  EllipseEx(Canvas.Handle, X, Y, A, B, Angle, N);
end;

procedure EllipseEx(DC: hDC; X, Y, A, B, Angle: Double; N: Integer);
{$ELSE (no MSWindows)}
procedure DrawEllipse(Canvas: TCanvas; X, Y, A, B, Angle: Double; N: Integer);
{$ENDIF MSWindows}
var
  P: array of T2DPoint;
  Q: array of TPoint;
  Alpha, Beta: Double;
  SinA, CosA: Double;
  Temp: Double;
  I, J, K: Integer;
begin
  if N < 2 then N := 2;
  Alpha := Pi/N;
  Beta := 4 * (1-Cos(Alpha)) / (3*Sin(Alpha));
  SetLength(P, N*3+1);
  for I := 0 to N-1 do
  begin
    Temp := I * (2*Alpha);
    CosA := cos(Temp);
    SinA := sin(Temp);
    K := I*3;
    if K = 0
      then J := N*3-1
      else J := K-1;
    P[K].X := A * CosA;
    P[K].Y := B * SinA;
    SinA := Beta * A * SinA;
    CosA := Beta * B * CosA;
    P[J].X := P[K].X + SinA;
    P[J].Y := P[K].Y - CosA;
    P[K+1].X := P[K].X - SinA;
    P[K+1].Y := P[K].Y + CosA;
  end;
  P[N*3] := P[0];
  CosA := cos(Angle);
  SinA := sin(Angle);
  SetLength(Q, N*3+1);
  for I := Low(P) to High(P) do
  // Rotate and translate
  begin
    Q[I].X := Round(P[I].X*CosA - P[I].Y*SinA +X);
    Q[I].Y := Round(P[I].X*SinA + P[I].Y*CosA +Y);
  end;
  {$IFDEF MSWindows}
  // Draw ellipse and fill it using the current pen and brush
  BeginPath(DC);
  PolyBezier(DC, Q, High(Q)+1);
  EndPath(DC);
  StrokeAndFillPath(DC);
  {$ELSE}
  for I := 0 to N-1 do
    Canvas.PolyBezier(Q, I * 3);
  {$ENDIF MSWindows}
end;
{$ENDIF RR_Overloading}

{$IFDEF MSWindows}
procedure DrawEllipse(Canvas: TCanvas; X, Y, A, B, Angle: Double);
begin
  EllipseEx(Canvas.Handle, X, Y, A, B, Angle);
end;

procedure EllipseEx(DC: hDC; X, Y, A, B, Angle: Double);
{$ELSE (no MSWindows)}
procedure DrawEllipse(Canvas: TCanvas; X, Y, A, B, Angle: Double);
{$ENDIF MSWindows}
const
  N = 8; // N Bezier curves

  P: array[0..N*3] of T2DPoint = (  // precomputed unit circle data
  (X:  1;	      Y:  0),
  (X:  1; 	      Y:  0.26521648984),
  (X:  0.89464315963; Y:  0.51957040274),
  (X:  0.70710678119; Y:  0.70710678119),

  (X:  0.51957040274; Y:  0.89464315963),
  (X:  0.26521648984; Y:  1),
  (X:  0; 	      Y:  1),

  (X: -0.26521648984; Y:  1),
  (X: -0.51957040274; Y:  0.89464315963),
  (X: -0.70710678119; Y:  0.70710678119),

  (X: -0.89464315963; Y:  0.51957040274),
  (X: -1;	      Y:  0.26521648984),
  (X: -1;	      Y:  0),

  (X: -1;	      Y: -0.26521648984),
  (X: -0.89464315963; Y: -0.51957040274),
  (X: -0.70710678119; Y: -0.70710678119),

  (X: -0.51957040274; Y: -0.89464315963),
  (X: -0.26521648984; Y: -1),
  (X:  0;	      Y: -1),

  (X:  0.26521648984; Y: -1),
  (X:  0.51957040274; Y: -0.89464315963),
  (X:  0.70710678119; Y: -0.70710678119),

  (X:  0.89464315963; Y: -0.51957040274),
  (X:  1;	      Y: -0.26521648984),
  (X:  1;	      Y: 0));
var
  Q: array[0..N*3] of TPoint;
  SinA, CosA: Double;
  I: Integer;
begin
  CosA := cos(Angle);
  SinA := sin(Angle);
  for i := Low(P) to High(P) do
  // Scale, rotate and translate
  begin
    Q[I].X := Round(P[I].X*A*CosA - P[I].Y*B*SinA +X);
    Q[I].Y := Round(P[I].X*A*SinA + P[I].Y*B*CosA +Y);
  end;
  {$IFDEF MSWindows}
  // Draw ellipse and fill it using the current pen and brush
  BeginPath(DC);
  PolyBezier(DC, Q, High(Q)+1);
  EndPath(DC);
  StrokeAndFillPath(DC);
  {$ELSE}
  for I := 0 to N-1 do
    Canvas.PolyBezier(Q, I * 3);
  {$ENDIF MSWindows}
end;

function FPoint(X, Y: Double): T2DPoint;
begin
  Result.X := X;
  Result.Y := Y;
end;

function NormAngle(const Angle: Double): Double;
begin
  Result := Frac(Angle * Inv2Pi);
  if Result < -0.5 then Result := Result + 1 else
  if Result >= 0.5 then Result := Result - 1;
  Result := Result*Pi2;
end;

procedure Rotate(var P: array of T2DVector; Angle: Double);
var
  I: Integer;
  X, Y: Double;
  SinA: Double;
  CosA: Double;
begin
  SinA := sin(Angle);
  CosA := cos(Angle);
  for I := Low(P) to High(P) do
  begin
    X := P[I].X;
    Y := P[I].Y;
    P[I].X := X * CosA - Y * SinA;
    P[I].Y := X * SinA + Y * CosA;;
  end;
end;

procedure Scale(var P: array of T2DVector; ScaleFactorX, ScaleFactorY: Double);
var
  I: Integer;
begin
  for I := Low(P) to High(P) do
  begin
    P[I].X := P[I].X * ScaleFactorX;
    P[I].Y := P[I].Y * ScaleFactorY;
  end;
end;

procedure Translate(var P: array of T2DVector; dX, dY: Double);
var
  I: Integer;
begin
  for I := Low(P) to High(P) do
  begin
    P[I].X := P[I].X + dX;
    P[I].Y := P[I].Y + dY;
  end;
end;

procedure FitBezierToUnitArc(var Points: T2DPointArray; BegAngle, EndAngle: Double);
var
  I, N: Integer;
var
  CosA, SinA: Double;
  Alpha, Beta: Double;
begin
  BegAngle := NormAngle(BegAngle);
  EndAngle := NormAngle(EndAngle);
  //if Abs(BegAngle - EndAngle) < 0.001 then //07.10.2011 if BegAngle = EndAngle then
  // IGOR 2017-05-03 òàê îíî ïðàâèëüíåå, à òî áûâàþò àðêè ñ îîî÷åíü áëèçêèìè óãëàìè 
  if Abs(BegAngle - EndAngle) < 0.00000001 then
  begin
    BegAngle := EndAngle;
    EndAngle := EndAngle + 2 * Pi;
  end;

  Alpha := 0.5 * (EndAngle - BegAngle);
  if Alpha < 0 then
    Alpha := Alpha + Pi;
  if DxfMode then
    N := Ceil(Alpha * (128 / Pi))
  else if DxfMode32 then
    N := Ceil(Alpha * (32 / Pi))
  else if DxfMode16 then
    N := Ceil(Alpha * (16 / Pi))
  else
    N := Ceil(Alpha * (4 / Pi));
  if N > 1 then
    Alpha := Alpha/N;

  if Alpha < Epsilon then Exit;
  SetLength(Points, N * 3 + 1);

  CosA := cos(Alpha);
  SinA := sin(Alpha);
  Beta := 4 * (1 - CosA) / (3 * SinA);

  for I := 0 to N - 1 do
  begin
    Points[I*3].X := CosA;
    Points[I*3].Y := -SinA;
    Points[I*3+3].X := CosA;
    Points[I*3+3].Y := SinA;
    Points[I*3+1].X := Points[I*3+0].X + Beta * SinA;
    Points[I*3+1].Y := Points[I*3+0].Y + Beta * CosA;
    Points[I*3+2].X := Points[I*3+3].X + Beta * SinA;
    Points[I*3+2].Y := Points[I*3+3].Y - Beta * CosA;
    Rotate(Slice(PCubicBezierCurve(@Points[I*3])^, 4),
           BegAngle+(I*2+1)*Alpha);
  end;

end;

{$IFDEF MSWindows}
procedure CircularArc(
  DC: hDC;
  X, Y,
  Radius,
  BegAngle,
  EndAngle: Double);
begin
  EllipticalArc(DC, X, Y, Radius, Radius, 0, BegAngle, EndAngle);
end;
{$ENDIF}

procedure DrawCircularArc(
  Canvas: TCanvas;
  X, Y,
  Radius,
  BegAngle,
  EndAngle: Double);
begin
  DrawEllipticalArc(Canvas, X, Y, Radius, Radius, 0, BegAngle, EndAngle);
end;

procedure DrawEllipticalArc(
  Canvas: TCanvas;
  X, Y,
  A, B,
  Angle,
  BegAngle,
  EndAngle: Double);
{$IFDEF MSWindows}
begin
  EllipticalArc(Canvas.Handle, X, Y, A, B, Angle, BegAngle, EndAngle);
end;

procedure EllipticalArc(
  DC: hDC;
  X, Y,
  A, B,
  Angle,
  BegAngle,
  EndAngle: Double);
{$ENDIF MSWindows}

{$IFDEF Debug}
  procedure CrossAt(DC: hDC; const P: TPoint);
  begin
    MoveToEx(DC, P.X-3, P.Y-3, nil);
    LineTo(DC, P.X+4, P.Y+4);
    MoveToEx(DC, P.X-3, P.Y+3, nil);
    LineTo(DC, P.X+4, P.Y-4);
  end;
{$ENDIF Debug}

  procedure ScaleAngle(var Angle: Double);
  var
    SinA, CosA: Extended;
  begin
    SinCos(Angle, SinA, CosA);
    Angle := ArcTan2(A * SinA, B * CosA);
  end;

var
  I: Integer;
  Points: T2DPointArray;
  P: array of TPoint;
begin
  if A <> B then
  begin
    ScaleAngle(BegAngle);
    ScaleAngle(EndAngle);
  end;
  FitBezierToUnitArc(Points,
        BegAngle,
        EndAngle);
  Scale(Points, A, B);
  Rotate(Points, Angle);
  Translate(Points, X, Y);
  SetLength(P, Length(Points));
  for I := 0 to Length(Points)-1 do
  begin
    P[I].X := Round(Points[I].X);
    P[I].Y := Round(Points[I].Y);
{$IFDEF Debug}
    CrossAt(DC, P[I]);
{$ENDIF Debug}
  end;
  {$IFDEF MSWindows}
  PolyBezier(DC, P[0], Length(P));
  {$ELSE}
  for I := 0 to Length(P) div 3-1 do
    Canvas.PolyBezier(P, I * 3);
  {$ENDIF}
end;

procedure BezierElpArcPoints(var Points: T2DPointArray; X, Y, aLen, bLen, Angle, BegAngle, EndAngle: Double);

  procedure ScaleAngle(var Angle: Double);
  var
    SinA, CosA: Extended;
  begin
    SinCos(Angle, SinA, CosA);
    Angle := ArcTan2(alen * SinA, blen * CosA);
  end;
  
begin
  if aLen <> bLen then
  begin
    ScaleAngle(BegAngle);
    ScaleAngle(EndAngle);
  end;
  FitBezierToUnitArc(Points, BegAngle, EndAngle);
  Scale(Points, aLen, bLen);
  Rotate(Points, Angle);
  Translate(Points, X, Y);
end;

procedure BezierArcPoints(var Points: T2DPointArray; X, Y, Radius, BegAngle, EndAngle: Double);
begin
  FitBezierToUnitArc(Points, BegAngle, EndAngle);
  Scale(Points, Radius, Radius);
  Translate(Points, X, Y);
end;


procedure GetCircleThruPoints(
  const P1, P2, P3: T2DPoint;
  var Center: T2DPoint;
  var Radius: Double);

  procedure GetMiddleOrthogonal(
    P1, P2: T2DPoint;
    var P: T2DPoint;
    var N: T2DVector);
  begin
    P.X := 0.5*(P1.X+P2.X);
    P.Y := 0.5*(P1.Y+P2.Y);
    N.X := P1.Y-P2.Y;
    N.Y := P2.X-P1.X;
  end;

var
  P12, P23: T2DPoint;
  N12, N23: T2DVector;
var
  a12, b: Double;
begin
  GetMiddleOrthogonal(P1, P2, P12, N12);
  GetMiddleOrthogonal(P2, P3, P23, N23);

  b := (N12.Y*N23.X-N12.X*N23.Y);
  if Abs(b) < Epsilon then raise EZeroDivide.Create(SPointsOnStraightLine);
  a12 := ((P23.Y-P12.Y) * N23.X + (P12.X-P23.X) * N23.Y)/b;

  Center.X := P12.X + a12 * N12.X;
  Center.Y := P12.Y + a12 * N12.Y;

  Radius := Sqrt(Sqr(Center.X-P1.X)+Sqr(Center.Y-P1.Y));
end;

{$IFDEF RR_Overloading}
procedure GetCircleThruPoints(
  const X1, Y1, X2, Y2, X3, Y3: Double;
  var CenterX, CenterY: Double;
  var Radius: Double);
var
  C: T2DPoint;
begin
  GetCircleThruPoints(
  	FPoint(X1, Y1),
        FPoint(X2, Y2),
        FPoint(X3, Y3),
        C, Radius);
  CenterX := C.X;
  CenterY := C.Y;
end;

procedure GetCircleThruPoints(
  const P1, P2, P3: TPoint;
  var Center: T2DPoint;
  var Radius: Double);
begin
  GetCircleThruPoints(
  	FPoint(P1.X, P1.Y),
        FPoint(P2.X, P2.Y),
        FPoint(P3.X, P3.Y),
        Center, Radius);
end;
{$ENDIF RR_Overloading}

end.