unit Geomet;
interface
uses
	//	VCL
	Windows, SysUtils, Graphics, Classes, Forms,

	//	CORE
	ArrayUtils, XLists, AClasses;

const
	MaxZoomScale: Double = 1e9;
	MinZoomScale: Double = 1/1e15;
	
	kMinZone : Double = 1e-9;
	kMinZone2 = 1e-12;

	minConvexity = 1e-4;
	minConvexity2 = 1e-15;

	Epsilon = 1E-15;

type
  {$A-}
	TDPointBase = object
		x, y: Double;
	end;
  {$A+}

	PDPoint = ^TDPoint;
  {A-}
	TDPoint = object(TDPointBase)
		procedure Assign(aX, aY: Double);
		procedure Offset(aDx, aDy: Double);
		function IsEqual(const aP: TDPointBase): Boolean;
		function DistanceFrom(const aP: TDPointBase): Double;
		function CloselyTo(const aP: TDPointBase; aDelta: Double): Boolean;
		function DistanceFromSegment(const aP1, aP2: TDPointBase): Double;
		function DistanceFromEdge(const aP1, aP2: TDPoint; aConvexity: Double = 0): Double;

		function DistanceFromLine(const aP1, aP2: TDPointBase): Double;
	end;
  {$A+}

	function DPointMult(const pt: TDPoint;m: Double): TDPoint;
	function DPointAdd(const pt1,pt2: TDPoint): TDPoint;
	function DPointSub(const pt1,pt2: TDPoint): TDPoint;
	function DPointScalar(const pt1,pt2: TDPoint): Double;
	function DPointModul(const pt: TDPoint): Double;
	function DPointNorm(const pt: TDPoint): TDPoint;

type
	TDLine = object
		P1, P2: TDPoint;
		function IsSingular: Boolean;
	end;

	PDArc = ^TDArc;
	TDArcBase = object(TDLine)
		C: Double;
	end;

	TDArc = object(TDArcBase)
		function GetClosestPoint(const aPoint: TDPoint): TDPoint;
		function IsLine: Boolean;
		function IsArc: Boolean;
		function IsEqual(const anEdge: TDArcBase): Boolean;
		function Contains(const aPoint: TDPoint; anEpsilon: Double = 0): Boolean;
	end;

	PDEdge = ^TDEdge;
	TDEdge = TDArc;


	PDRect = ^TDRect;
  {$A-}
	TDRect = object
		x1, y1, x2, y2: Double;

		procedure Assign(aX1, aY1, aX2, aY2: Double);
		procedure AssignPoints(const aP1, aP2: TDPoint);
		procedure AssignFrom(const aPoints: array of TDPoint);

		function DX: Double;
		function DY: Double;

		function IsNormalized: Boolean;
		procedure Offset(aDx, aDy: double);
		procedure Inflate(aDx, aDy: double);
		procedure Expand(const aPoint: TDPoint);
		procedure Normalized;
		procedure Clear;

		function Center: TDPoint;
		function Contains(const aPoint: TDPoint): Boolean;
	end;
  {$A+}

	TTransformSize = function(aValue: Double): Double of object;
	TTransformPoint = function(aValue: TDPoint): TDPoint of object;

	{ TDPointList }

	TDPointList = class
	private
		List: TRecList;

		function GetPoint(I: Integer): TDPoint;
		procedure SetPoint(I: Integer; const aValue: TDPoint);
		function GetCount: Integer;
	public
		Closed: Boolean;			{ Для полигона Closed = True }

		property Points[I: Integer]: TDPoint read GetPoint write SetPoint; default;
		property Count: Integer read GetCount;

		function GetPrevIndex(anIndex: Integer): Integer;
		function GetNextIndex(anIndex: Integer): Integer;

		procedure Assign(aSrc: TDPointList);
		procedure AssignPoints(aPoints: PDPoint; aCount: Integer);

		procedure AddPoint(const aPoint: TDPoint);
		procedure AddPoints(aPointList: TDPointList);
		procedure InsertPoint(anIndex: Integer; const aPoint: TDPoint);
		procedure DeletePoint(Index: Integer);
		procedure Clear;

		function GetSquare: Double;
		function EdgeAtPoint(const aPoint: TDPoint; aDistance: Double): Integer;
		function LineTouched(aP1, aP2: TDPoint): Boolean;
		function GetDistanceFromPoint(aPoint: TDPoint): Double;
		function GetBounds(var aBounds: TDRect): Boolean;

		function GetIntersectionCount(aPoint: TDPoint): Integer;
		function ContainsPoint(aPoint: TDPoint): Boolean;

		procedure LoadFromStream(aStream: TStream);
		procedure SaveToStream(aStream: TStream);

		procedure DeleteSingularSegments;

		constructor Create;
		destructor Destroy; override;
	end;

	TDPointsCompare = function(const A, B: TDPoint; aData: Pointer): Integer;

	TDPointArray = object
	private
		FItems: array of TDPoint;
		FCount: Integer;

	public
		procedure Clear(aCapacity: Integer = 0);
		function  Add(const aPoint: TDPoint): Integer; overload;
		procedure Add(const aPoints: TDPointArray); overload;

		procedure Sort(aCompareFunction: TDPointsCompare; aData: Pointer);

		function GetItem(Index: Integer): TDPoint;
		function Count: Integer;
	end;


{ CoordSys }

const
	M_PI = Pi;
	M_PI_2 = M_PI / 2;
	M_PI_4 = M_PI / 4;

type
	TCoordSys = object
	private
		FInitByPoints: Boolean;
		FDX, FDY: Extended;
		FHypot: Extended;

	public
		Center, EndPoint: TDPointBase;
		AngleValue: Extended;
		AngleCalculated: Boolean;
		AngleSin: Extended;
		AngleCos: Extended;

		procedure InitPoints(const aCenter, anEnd: TDPointBase);
		procedure InitAngle(const aCenter: TDPoint; anAngle: Extended);

		procedure MakeLocal(var p: TDPoint);
		procedure MakeGlobal(var p: TDPoint);

		function ToLocal(const p: TDPointBase): TDPoint;
		function ToGlobal(const p: TDPoint): TDPoint;

		function GetXProjection(const aGlobalP: TDPoint): Extended;
		function GetYProjection(const aGlobalP: TDPoint): Extended;

		function MakeGlobalFromXY(x, y: double): TDPoint;

		function IsValid: Boolean;
		procedure SetInvalid;
		function Angle: Extended;
	end;

	function DPointFromPoint(aPoint: TPoint): TDPoint;
	function DRectFromRect(const aRect: TRect): TDRect;
	function RectFromDRect(const aR: TDRect): TRect;
	function PointFromDPoint(const aDPoint: TDPoint): TPoint;

	function DPoint(aX, aY: Double): TDPoint;
	function DLine(aP1, aP2: TDPoint): TDLine;
	function DArc(aP1, aP2: TDPoint; aC: Double): TDArc;
	function DRect(aX1, aY1, aX2, aY2: Double): TDRect;
	function DEdge(const aP1, aP2: TDPoint; aC: Double = 0): TDEdge;

	function DEdgeIsEqual(const anEdge1, anEdge2: TDEdge): Boolean;
	function DEdgeReverse(const anEdge: TDEdge): TDEdge;

	function DPointOffsetBy(const aP: TDPoint; aDX, aDY: Double): TDPoint;
	function DRectIsEqual(const aR1, aR2: TDRect): Boolean;
	procedure DRectIntersectWith(var aR1: TDRect; const aR2: TDRect);
	procedure DRectUnionWith(var aR1: TDRect; const aR2: TDRect);
	procedure DRectExpandByRect(var aR1: TDRect; const aR2: TDRect);
	procedure DRectExpandByExtentions(var aR1: TDRect; const aR2: TDRect);

	function DRectIntersectedWith(const aR1, aR2: TDRect): Boolean; overload;
	function DRectIntersectedWith(aRect: TDRect; const anEdge: TDEdge): Boolean; overload;
	function DRectContainsRect(const aR1, aR2: TDRect): Boolean;


	function LinePointPosition(aLine: TDLine; aPoint: TDPoint): Integer;

	function IsPointNearEdge(const p, a, b: TDPoint; aDistance: Double): Boolean;
	function GetXDIntersection(P1, P2: TDPoint; P: TDPoint; var X: Double): Boolean;

	procedure MakeOutlinePolygon(aSrc, aDst: TDPointList; aDelta: Double);
	function ZoomScaleInRange(aZoomScale, aMinScale, aMaxScale: Double): Boolean;

	function ArcContainsPoint(const ptC,pt1: TDPoint; delta: Double;
																				 const ptTest: TDPoint): Boolean;
	{	Возвращает True, если есть пересечение между двумя отрезками.
		По ссылке возвращает точку пересечения	}
	function Intersect(aBP1, aEP1, aBP2, aEP2: TDPoint; var aPoint: TDPoint): Boolean;

	function IntersectLineLine(const aLine1, aLine2: TDLine; var aP1, aP2: TDPoint): Integer; overload;
	function IntersectLineLine(const aLine1, aLine2: TDLine; var aP1, aP2: TDPoint;
		out aTouched: Boolean): Integer; overload;
	function IntersectLineArc(const aLine: TDLine; const anArc: TDArc; var aP1, aP2: TDPoint): Integer;
	function IntersectArcArc(const anArc1, anArc2: TDArc; var aP1, aP2: TDPoint): Integer;
	function IntersectEdgeEdge(const anEdge1, anEdge2: TDEdge; var aP1, aP2: TDPoint; out aTouched: Boolean): Integer;
	function InnerIntersectEdgeEdge(const anEdge1, anEdge2: TDEdge; var aP1, aP2: TDPoint): Integer;
	function IntersectLineLineCircle(const aLine: TDLine;
		const aCenter: TDPoint; aRadius: Double; var aP1, aP2: TDPoint): Integer;
	function IntersectCircleCircle(const aCenter1: TDPoint; aRadius1: Double;
		const aCenter2: TDPoint; aRadius2: Double; var aP1, aP2: TDPoint): Integer;
	function IntersectLineLineSegment(const aLine, aSegment: TDLine; var aP: TDPoint): Integer;



	procedure ClipLineByRect(const aRect: TDRect; const aLine: TDLine; out aClipLine: TDLine);

	function NormalizeDirection(anAngle: Extended): Extended;
	function GetAngle(const a, b: TDPoint): Double;
function GetCenter(const a, b, c: TDPoint; var delta: Double; out aResult: TDPoint): Boolean;
function GetCenter2(const a, b: TDPoint; aConvexity: Double; var delta: Double;
	out aResult: TDPoint): Boolean;
	function GetArcPoint(aP1: TDPoint; anAngle: Double; aP2: TDPoint): TDPoint;
	function GetArcPoint2(aP1, aCenter, aP2: TDPoint; aClockwise: Boolean): TDPoint;
function GetConvexity(aP1, aP, aP2: TDPoint; out aResult: Double): Boolean;
	function GetConvexity2(aP1, aCenter, aP2: TDPoint; aClockwise: Boolean): Double;
	function ConvexityFromAngle(anAngle: Double): Double;
	function AngleInRange(anAngle, aBeginAngle, aDelta: Double): Boolean;
	procedure CheckResult3(const aP1, aP2, aP3, aP4: TDPoint; numPoints: Integer;
		var pt1,pt2: TDPoint);

function ZoomScaleToStr(aZoomScale: Extended): String;
function ZoomScaleRangeToStr(aMinScale, aMaxScale: Double): String;
function AngleToStr(anAngleInRad: Double; anAccurateFormat: Boolean = False): String;
function StrToAngle(aText: String): Double;
function StrIsAngle(s: String): Boolean;
procedure DecodeAngle(v: Extended; out aDeg, aMin, aSec: Integer);
function EncodeAngle(aDeg, aMin, aSec: Extended): Extended;

function SetDPoint(const anOldValue, aNewValue: TDPoint; var aHasChanges: Boolean): TDPoint;
function SetDRect(const anOldValue, aNewValue: TDRect; var aHasChanges: Boolean): TDRect;


type
	TSPRelations = (sprLeft,  sprRight,  sprBeyond,  sprBehind, sprBetween, sprOrigin, sprDestination);
	TLPRelations = (lprLeft,  lprRight,  lprContain);

function DSegmentContaintPoint(const aSegment: TDLine; const aPoint: TDPoint; anEpsilon: Double = 0): Boolean;
function DSegmentPointRelation(const aSegment: TDLine; const aPoint: TDPoint): TSPRelations;
function DLinePointRelation(const aLine: TDLine; const aPoint: TDPoint): TLPRelations;


var
	dMinZone: Double;
	
implementation

uses
	Math, Versions,
	InSysUtils;


// корректность не проверяется, на входе обязательно должен быть не вырожденный отрезок

function DSegmentContaintPoint(const aSegment: TDLine; const aPoint: TDPoint; anEpsilon: Double = 0): Boolean;
var
	aDX, aDY: Double;
	aP: Double;
begin
	aDX := aSegment.P1.X - aSegment.P2.X;
	aDY := aSegment.P1.Y - aSegment.P2.Y;
	if Abs(aDX) > Abs(aDY) then begin
		aP := (aPoint.x - aSegment.P2.x)/aDX;
		if (aP >= 0) and (aP <= 1) and	(Abs(aP*aDY - (aPoint.y - aSegment.P2.y)) <= anEpsilon) then
			Result := True
		else
			Result := False
	end else begin
		aP := (aPoint.y - aSegment.P2.y)/aDY;
		if (aP >= 0) and (aP <= 1) and (Abs(aP*aDX - (aPoint.x - aSegment.P2.x)) <= anEpsilon) then
			Result := True
		else
			Result := False;
	end;
end;

function DArcContaintPoint(const anArc: TDArc; const aPoint: TDPoint; anEpsilon: Double = 0): Boolean;
var
	aCenter: TDPoint;
	aDelta: Double;
begin
	if GetCenter2(anArc.P1, anArc.P2, anArc.C, aDelta, aCenter) then begin
		if ArcContainsPoint(aCenter, anArc.P1, aDelta, aPoint) then
			Result := Abs(aCenter.DistanceFrom(aPoint) - aCenter.DistanceFrom(anArc.P1)) <= anEpsilon
		else
			Result := False;
	end else begin
		Assert(False,
			format('Invalid params, can''t calculate arccenter. Convexity: %f', [anArc.C]));

		Result := False;
	end;
end;

{ TDLine }

function TDLine.IsSingular: Boolean;
begin
	Result := P1.IsEqual(P2);
end;

{ TDArc }

function TDArc.Contains(const aPoint: TDPoint; anEpsilon: Double = 0): Boolean;
begin
	if IsSingular then
		Result := aPoint.IsEqual(aPoint)
	else if IsLine then
		Result := DSegmentContaintPoint(Self, aPoint, anEpsilon)
	else if IsArc then
		Result := DArcContaintPoint(Self, aPoint, anEpsilon)
	else
		Result := False
end;

function TDArc.GetClosestPoint(const aPoint: TDPoint): TDPoint;
var
	aCoordSys: TCoordSys;
	anAlfa: Double;
	aSin, aCos: Extended;
	aCenter: TDPoint;
	aDelta: Double;
	aRadius: Double;
begin
	if IsArc then begin
		if GetCenter2(P1, P2, C, aDelta, aCenter) then
			aRadius := aCenter.DistanceFrom(P2)
		else
			aRadius := 0;
		anAlfa := GetAngle(aCenter, aPoint);
		if AngleInRange(anAlfa, GetAngle(aCenter, P1), aDelta) then begin
			SinCos(anAlfa, aSin, aCos);
			Result := DPoint(aCenter.x + aRadius*aCos, aCenter.y + aRadius*aSin)
		end else begin
			Result := P1;
			if aPoint.DistanceFrom(P1) > aPoint.DistanceFrom(P2) then
				Result := P2
		end
	end else begin
		aCoordSys.InitPoints(P1, P2);
		Result := aPoint;
		aCoordSys.MakeLocal(Result);
		if Result.X < 0 then
			Result := P1
		else if Result.X > aCoordSys.ToLocal(P2).X then
			Result := P2
		else
			Result := aCoordSys.ToGlobal(DPoint(Result.x, 0));
	end;
end;

function TDArc.IsArc: Boolean;
begin
	Result := Abs(C) > minConvexity2;
end;

function TDArc.IsLine: Boolean;
begin
	Result := Abs(C) <= minConvexity2;
end;

function TDArc.IsEqual(const anEdge: TDArcBase): Boolean;
var
	aC1, aC2: Double;
begin {todo: ???}
	aC1 := C;
	aC2 := anEdge.C;
	Result := (aC1 = aC2) and P2.IsEqual(anEdge.P2) and P1.IsEqual(anEdge.P1) or
						(aC1 = -aC2) and P1.IsEqual(anEdge.P2) and P2.IsEqual(anEdge.P1);
end;


{ выдает положение точки относительно прямой}

function DSegmentPointRelation(const aSegment: TDLine; const aPoint: TDPoint): TSPRelations;
var
	a, b: TDPoint;
	s: double;
begin
	a.Assign(aSegment.P2.x - aSegment.P1.x, aSegment.P2.y - aSegment.P1.y);
	b.Assign(aPoint.x - aSegment.P1.x, aPoint.y - aSegment.P1.y);
	s := a.x * b.y - b.x * a.y;
	if s > 0 then
		Result := sprLeft
	else if (s < 0) then
		Result := sprRight
	else if ((a.x * b.x < 0) or (a.y * b.y < 0)) then
		Result := sprBehind
	else if (sqr(a.x) + sqr(a.y)) < (sqr(b.x) + sqr(b.y)) then
		Result := sprBeyond
	else if aPoint.IsEqual(aSegment.P1) then
		Result := sprOrigin
	else if aPoint.IsEqual(aSegment.P2) then
		Result := sprDestination
	else
		Result := sprBetween;
end;

function DLinePointRelation(const aLine: TDLine; const aPoint: TDPoint): TLPRelations;
var
	a, b: TDPoint;
	s: double;
begin
	a.Assign(aLine.P2.x - aLine.P1.x, aLine.P2.y - aLine.P1.y);
	b.Assign(aPoint.x - aLine.P1.x, aPoint.y - aLine.P1.y);
	s := a.x * b.y - b.x * a.y;
	if (s > 0) then
		Result := lprLeft
	else if (s < 0) then
		Result := lprRight
	else
		Result := lprContain;
end;

function ZoomScaleToStr(aZoomScale: Extended): String;
begin
	if aZoomScale = 0 then
		Result := '0'
	else if aZoomScale < 1 then
		Result := FloatToStrF(1/aZoomScale, ffFixed, 10, 0)
	else
		Result := FloatToStrF(1/aZoomScale, ffGeneral, 1, 4)
end;

function ZoomScaleRangeToStr(aMinScale, aMaxScale: Double): String;
begin
	if aMinScale < 0 then
		Result := '('
	else
		Result := '[';

	Result := Result + ZoomScaleToStr(Abs(aMinScale)) + '-' + ZoomScaleToStr(Abs(aMaxScale));

	if aMaxScale < 0 then
		Result := Result + ')'
	else
		Result := Result + ']';
end;

procedure DecodeAngle(v: Extended; out aDeg, aMin, aSec: Integer);
begin
	if v < 0 then
		v := v + 2*PI*(Int(abs(v/(2*PI)))+1);

	v := RadToDeg(v);
	aDeg := Trunc(v);
	aMin := Trunc(60*Frac(v));
	aSec := Round(60*Frac(60*Frac(v)));
	if aSec >= 60 then begin
		aMin := aMin + aSec div 60;
		aSec := aSec mod 60;
	end;
	if aMin >= 60 then begin
		aDeg := aDeg + aMin div 60;
		aMin := aMin mod 60;
	end;
	aDeg := aDeg mod 360;
end;

function EncodeAngle(aDeg, aMin, aSec: Extended): Extended;
begin
	Result := DegToRad(aDeg + aMin/60 + aSec/3600);
end;

function AngleToStr(anAngleInRad: Double; anAccurateFormat: Boolean = False): String;
var
	g, m, s: Integer;
begin
	DecodeAngle(anAngleInRad, g, m, s);
	if anAccurateFormat then
		Result := Format('%3.3d°', [g])
	else
		Result := Format('%d°', [g]);

	if (m <> 0) and (s = 0) then
		Result := Result + Format('%2.2d''', [m])
	else if (s <> 0) then
		Result := Result + Format('%2.2d''%2.2d"', [m, s]);

	if anAccurateFormat and (g >= 0) then
		Result := '+' + Result;
end;

function StrToAngle(aText: String): Double;
var
	aPos: Integer;
	d, m, s: Extended;
begin
	aText := Trim(aText);
	try
		m := 0;
		s := 0;
		aPos := Pos('°', aText);
		if aPos = 0 then aPos := Pos(' ', aText);
		if aPos = 0 then
			d := StringToFloat(aText)
		else begin
			d := StringToFloat(Trim(Copy(aText, 1, aPos - 1)));
			aText := Trim(Copy(aText, aPos + 1, Length(aText) - aPos));
			aPos := Pos('''', aText);
			if aPos = 0 then aPos := Pos(' ', aText);
			if aPos = 0 then
				if aText <> '' then
					m := StringToFloat(aText)
				else
			else begin
				m := StringToFloat(Trim(Copy(aText, 1, aPos - 1)));
				aText := Trim(Copy(aText, aPos + 1, Length(aText) - aPos));
				aPos := Pos('"', aText);
				if aPos = 0 then aPos := Pos('''''', aText);
				if aPos <> 0 then
					aText := Copy(aText, 1, aPos - 1);
				if aText <> '' then
					s := StringToFloat(aText);
			end
		end;
		Result := EncodeAngle(d, m, s)
	except
		raise Exception.Create('Некорректное значение угла')
	end
end;

function StrIsAngle(s: String): Boolean;
var
	aPos: Integer;
begin
	s := Trim(s);
	try
		aPos := Pos('°', s);
		if aPos = 0 then aPos := Pos(' ', s);
		if aPos = 0 then
			Result := StrIsFloat(s)
		else begin
			Result := StrIsFloat(Trim(Copy(s, 1, aPos - 1)));
			if Result then begin
				s := Trim(Copy(s, aPos + 1, Length(s) - aPos));
				aPos := Pos('''', s);
				if aPos = 0 then aPos := Pos(' ', s);
				if aPos = 0 then
					Result := (s = '') or StrIsFloat(s)
				else begin
					Result := StrIsFloat(Trim(Copy(s, 1, aPos - 1)));
					if Result then begin
						s := Trim(Copy(s, aPos + 1, Length(s) - aPos));
						aPos := Pos('"', s);
						if aPos = 0 then aPos := Pos('''''', s);
						if aPos <> 0 then
							s := Trim(Copy(s, 1, aPos - 1));
						Result := (s = '') or StrIsFloat(s);
					end
				end
			end
		end;
	except
		Result := False
	end
end;
{
procedure PolyPolyline(DC: HDC; var Points; var PolyCounts; Count: Integer);
var
	aPolyCounts: ^Integer;
	aPPoint: ^TPoint;
	i: Integer;
begin
	aPPoint := @Points;
	aPolyCounts := @PolyCounts;
	for i := 0 to Count - 1 do begin
		Polyline(DC, aPPoint^, aPolyCounts^);
		Inc(aPPoint, aPolyCounts^);
		Inc(aPolyCounts);
	end
end;
 }

{ Возвращает угол в радианах от 0 до 2*Pi }
function NormalizeDirection(anAngle: Extended): Extended;
var
	Pi2: Extended;
begin
	Pi2 := Pi * 2;
	Result := anAngle;
	while Result > Pi2 do
		Result := Result - Pi2;
	while Result < 0 do
		Result := Result + Pi2;
end;

{	Возвращает угол вектора от точки a к точке b в радианах	от 0 до 2*Pi	}
function GetAngle(const a, b: TDPoint): Double;
begin
	Result := ArcTan2((b.y - a.y), (b.x - a.x));
	if Result < 0 then Result := Result + 2*Pi;
end;

{	Возвращает центр окружности, заданной тремя точками. delta - это угол,	}
{	на который надо нарастить угол вектора из центра к точке окружности,		}
{	чтобы прийти из точки a в точку c через точку b	}
function GetCenter(const a, b, c: TDPoint; var delta: Double; out aResult: TDPoint): Boolean;
var
	aCoordSys: TCoordSys;
	aB, aC: TDPoint;
	alfa, beta, gama: Double;
begin
	aCoordSys.InitPoints(DPoint((a.x + b.x)/2, (a.y + b.y)/2), b);
	aB := aCoordSys.ToLocal(b);
	aC := aCoordSys.ToLocal(c);
	Result := aC.y <> 0;
	if Result then begin
		aResult := DPoint(0, (Sqr(aC.x) + Sqr(aC.y) - Sqr(aB.x))/2/aC.y);
		alfa := GetAngle(aResult, aCoordSys.ToLocal(a));
		gama := GetAngle(aResult, aB);
		beta := GetAngle(aResult, aC);
		if alfa > beta then beta := beta + 2*Pi;
		if (gama > alfa)and(gama < beta)or(gama+2*Pi > alfa)and(gama+2*Pi < beta) then
			delta := beta - alfa
		else
			delta := beta - 2*Pi - alfa;
		aCoordSys.MakeGlobal(aResult);
	end;
end;

{	Возвращает центр дуги, заданной двумя точками и выпуклостью. delta - это угол,	}
{	на который надо нарастить угол вектора из центра к точке окружности,						}
{	чтобы прийти из точки a в точку c через точку b																	}
function GetCenter2(const a, b: TDPoint; aConvexity: Double; var delta: Double;
	out aResult: TDPoint): Boolean;
var
	aCoordSys: TCoordSys;
	anAngle, aDistance: Double;
begin
	Result := (Abs(aConvexity) > minConvexity2) and not a.IsEqual(b);
	if Result then begin
		aCoordSys.InitPoints(a, b);
		aDistance := aCoordSys.ToLocal(b).x / 2;
		if Abs(aConvexity) > 1e-2 then
			aResult := aCoordSys.ToGlobal(DPoint(aDistance, aDistance * (1 - Sqr(aConvexity))/ (2*aConvexity)))
		else begin
			anAngle := 2 * ArcTan(-aConvexity);
			if anAngle > 0 then anAngle := anAngle - Pi/2
				else anAngle := anAngle + Pi/2;
			aResult := aCoordSys.ToGlobal(DPoint(aDistance, aDistance * Tan(anAngle)))
		end;
		aCoordSys.InitPoints(aResult, a);
		delta := GetAngle(DPoint(0, 0), aCoordSys.ToLocal(b));
		if aConvexity > 0 then delta := delta - 2*Pi;
		if delta > 0 then delta := delta - 2*Pi
			else delta := delta + 2*Pi;
	end
end;

{	Возвращает точку, лежащую на дуге из точки aP1 в точку aP2 (по середине), }
{	причем из точки aP1 дуга выходит под углом anAngle	}
function GetArcPoint(aP1: TDPoint; anAngle: Double; aP2: TDPoint): TDPoint;
var
	anY, alfa: Double;
	aP: TDPoint;
	aCoordSys: TCoordSys;
	aSin, aCos: Extended;
begin
	aCoordSys.InitAngle(aP1, anAngle);
	aP := aCoordSys.ToLocal(aP2);
	if Abs(aP.y) < kMinZone2 then begin
		Result := DPoint((aP1.x + aP2.x)/2, (aP1.y + aP2.y)/2);
		Exit
	end;
	anY := (Sqr(aP.x) + Sqr(aP.y))/2/aP.y;
	if anY > 0 then
		alfa := Pi/2 + GetAngle(DPoint(0, anY), aP)
	else
		alfa := Pi/2 - GetAngle(DPoint(0, anY), aP);
	if alfa < 0 then alfa := alfa + 2*Pi;
	if alfa >= 2*Pi then alfa := alfa - 2*Pi;
	SinCos(alfa / 2, aSin, aCos);
	Result := aCoordSys.ToGlobal(DPoint(Abs(anY)*aSin, anY*(1 - aCos)));
end;

{	Возвращает точку, лежащую на дуге из точки aP1 в точку aP2 (по середине), }
{	с центром в точке aCenter. Причем если aClockwise = True, то из точки aP1 }
{	дуга идет в точку aP2 по часовой стрелке, иначе - против часовой стрелки.	}
{	Считается, что точка aP2 лежит на том же расстоянии от центра, что и aP1.	}
{	В противном случае радиус усредняется.	}
function GetArcPoint2(aP1, aCenter, aP2: TDPoint; aClockwise: Boolean): TDPoint;
var
	alfa, aRadius: Double;
	aP: TDPoint;
	aCoordSys: TCoordSys;
begin
	aCoordSys.InitPoints(aCenter, aP1);
	aP := aCoordSys.ToLocal(aP2);
	alfa := GetAngle(DPoint(0, 0), aP) / 2;
	if aClockwise then
		alfa := alfa + Pi;
	aRadius := (Sqrt(Sqr(aP.x) + Sqr(aP.y)) + aCoordSys.ToLocal(aP1).x)/2;
	Result := aCoordSys.ToGlobal(DPoint(aRadius*cos(alfa), aRadius*sin(alfa)));
end;

{	Возвращает выпуклость дуги, заданной тремя точками	}
function GetConvexity(aP1, aP, aP2: TDPoint; out aResult: Double): Boolean;
var
	aDelta: Double;
	aCenter: TDPoint;
	aCoordSys: TCoordSys;
begin
	Result := True;
	if aP1.IsEqual(aP) or aP.IsEqual(aP2) or aP2.IsEqual(aP1) then
		aResult := 0
	else begin
		aCoordSys.InitPoints(aP1, aP2);
		Result := GetCenter(aP1, aP, aP2, aDelta, aCenter);
		if Result then
			aResult := GetConvexity2(aP1, aCenter, aP2, aDelta < 0)
	end
end;

{	Возвращает выпуклость дуги, заданной двумя точками и центром }
function GetConvexity2(aP1, aCenter, aP2: TDPoint; aClockwise: Boolean): Double;
var
	anAngle: Double;
	aCoordSys: TCoordSys;
	aSin, aCos: Extended;
begin
	aCoordSys.InitPoints(aP1, aP2);
	aCenter := aCoordSys.ToLocal(aCenter);
	if aCenter.y <> 0 then begin
		anAngle := ArcTan(aCenter.x / aCenter.y);
		if (aCenter.y > 0) = aClockwise then
			if anAngle > 0 then anAngle := anAngle - Pi
				else anAngle := anAngle + Pi;
	end else
		if aClockwise then
			anAngle := -Pi/2
		else
			anAngle := Pi/2;
	anAngle := anAngle / 2;
	SinCos(anAngle, aSin, aCos);
	if aCos <> 0 then
		Result := aSin / aCos
	else
		Result := 0
end;

function ConvexityFromAngle(anAngle: Double): Double;
var
	aSin, aCos: Extended;
begin
	SinCos(anAngle / 2, aSin, aCos);
	if aCos <> 0 then
		Result := aSin / aCos
	else
		Result := 0;
end;

//	Попадает ли угол anAngle в диапазон дуги { BeginPoint, EndPoint, ExtraData	}
//	aCenter = GetCenter2(BeginPoint, EndPoint, PDouble(ExtraData)^, aDelta)
//	aBeginAngle = GetAngle(aCenter, BeginPoint)
function AngleInRange(anAngle, aBeginAngle, aDelta: Double): Boolean;
var
	anEndAngle: Double;
begin
	anEndAngle := aBeginAngle + aDelta;
	if aDelta > 0 then
		Result := (anAngle >= aBeginAngle) and (anAngle <= anEndAngle) or
							(anAngle + 2*Pi >= aBeginAngle) and (anAngle + 2*Pi <= anEndAngle) or
							(anAngle - 2*Pi >= aBeginAngle) and (anAngle - 2*Pi <= anEndAngle)
	else
		Result := (anAngle <= aBeginAngle) and (anAngle >= anEndAngle) or
							(anAngle + 2*Pi <= aBeginAngle) and (anAngle + 2*Pi >= anEndAngle) or
							(anAngle - 2*Pi <= aBeginAngle) and (anAngle - 2*Pi >= anEndAngle);
end;

function Intersect(aBP1, aEP1, aBP2, aEP2: TDPoint; var aPoint: TDPoint): Boolean;
var
	aCoordSys: TCoordSys;
	aP1, aP2: TDPoint;
begin
	Result := False;
	aCoordSys.InitPoints(aBP1, aEP1);
	aP1 := aCoordSys.ToLocal(aBP2);
	aP2 := aCoordSys.ToLocal(aEP2);
	if aP1.y * aP2.y > 0 then Exit;
	aPoint.y := 0;
	if aP1.y = aP2.y then
		aPoint.x := aP1.x
	else
		aPoint.x := aP1.x + aP1.y*(aP2.x - aP1.x)/(aP1.y - aP2.y);
	Result := (aPoint.x >= 0) and (aPoint.x <= aCoordSys.ToLocal(aEP1).x);
	aCoordSys.MakeGlobal(aPoint);
end;

function IntersectLineLine(const aLine1, aLine2: TDLine; var aP1, aP2: TDPoint): Integer;
var
	aTouched: Boolean;
begin
	Result := IntersectLineLine(aLine1, aLine2, aP1, aP2, aTouched);
end;

function IntersectLineLine(const aLine1, aLine2: TDLine; var aP1, aP2: TDPoint;
	out aTouched: Boolean): Integer;
var
	aDx1, aDx2, aDy1, aDy2, aT1, aT2, aNumerator, aDenominator: Double;

	function GetPoint: Integer;
	begin
		if Abs(aDx1) < kMinZone then
			aT1 := (aLine2.P1.y - aLine1.P1.y + aDy2*aT2)/aDy1
		else
			aT1 := (aLine2.P1.x - aLine1.P1.x + aDx2*aT2)/aDx1;
		if (aT2 >= 0) and (aT2 <= 1) and (aT1 >= 0) and (aT1 <= 1) then begin
			Result := 1;
			if Abs(aT2) < 0.0000000001 then begin
				aTouched := True;
				aP1.x := aLine2.P1.x;
				aP1.y := aLine2.P1.y;
			end else if aT2 > 0.9999999999 then begin
				aTouched := True;
				aP1.x := aLine2.P2.x;
				aP1.y := aLine2.P2.y;
			end else
				if Abs(aT1) < 0.0000000001 then begin
					aTouched := True;
					aP1.x := aLine1.P1.x;
					aP1.y := aLine1.P1.y;
				end else if aT1 > 0.9999999999 then begin
					aTouched := True;
					aP1.x := aLine1.P2.x;
					aP1.y := aLine1.P2.y;
				end else begin
					aTouched := True;
					aP1.x := aLine1.P1.x + aT1*aDx1;
					aP1.y := aLine1.P1.y + aT1*aDy1;
				end
		end else
			Result := 0;
	end;

	function Match: Integer;

		function PointInLine(aPoint: TDPoint; aLine: TDLine): Boolean;
		var
			aT: Double;
		begin
			if Abs(aLine.P2.x - aLine.P1.x) < kMinZone then
				aT := (aPoint.y - aLine.P1.y)/(aLine.P2.y - aLine.P1.y)
			else
				aT := (aPoint.x - aLine.P1.x)/(aLine.P2.x - aLine.P1.x);
			Result := (aT >= 0) and (aT <= 1)
		end;

	begin
		aTouched := True;
		Result := 0;
		if PointInLine(aLine1.P1, aLine2) then begin
			Inc(Result);
			aP1 := aLine1.P1;
		end;
		if PointInLine(aLine1.P2, aLine2) then begin
			Inc(Result);
			if Result = 1 then aP1 := aLine1.P2 else aP2 := aLine1.P2
		end;
		if Result = 2 then Exit;
		if PointInLine(aLine2.P1, aLine1) then begin
			Inc(Result);
			if Result = 1 then aP1 := aLine2.P1 else aP2 := aLine2.P1
		end;
		if Result = 2 then Exit;
		if PointInLine(aLine2.P2, aLine1) then begin
			Inc(Result);
			if Result = 1 then aP1 := aLine2.P2 else aP2 := aLine2.P2
		end;
	end;

var
	aSingular1, aSingular2: Boolean;
begin
	Result := 0;
	aTouched := False;

	aSingular1 := aLine1.IsSingular;
	aSingular2 := aLine2.IsSingular;
	if aSingular1 and aSingular2 then begin
		if aLine1.P1.IsEqual(aLine2.P1) then begin
			aTouched := True;
			aP1 := aLine1.P1;
			Result := 1;
		end;
	end else if aSingular1 then begin
		if DSegmentContaintPoint(aLine2, aLine1.P1, kMinZone) then begin
			aTouched := True;
			aP1 := aLine1.P1;
			Result := 1;
		end;
	end else if aSingular2 then begin
		if DSegmentContaintPoint(aLine1, aLine2.P1, kMinZone) then begin
			aTouched := True;
			aP1 := aLine2.P1;
			Result := 1;
		end;
	end
	else begin
		aDx1 := aLine1.P2.x - aLine1.P1.x;
		aDx2 := aLine2.P2.x - aLine2.P1.x;
		aDy1 := aLine1.P2.y - aLine1.P1.y;
		aDy2 := aLine2.P2.y - aLine2.P1.y;
		if Abs(aDx1) < kMinZone then
			if Abs(aDx2) < kMinZone then
				if Abs(aLine1.P1.x - aLine2.P1.x) < kMinZone then
					Result := Match
				else
					Result := 0
			else begin
				aT2 := (aLine1.P1.x - aLine2.P1.x)/aDx2;
				Result := GetPoint
			end
		else begin
			aNumerator := aLine1.P1.y - aLine2.P1.y + aDy1*(aLine2.P1.x - aLine1.P1.x)/aDx1;
			aDenominator := aDy2 - aDy1*aDx2/aDx1;
			if Abs(aDenominator) < kMinZone then
				if Abs(aNumerator) < kMinZone then
					Result := Match
				else
					Result := 0
			else begin
				aT2 := aNumerator/aDenominator;
				Result := GetPoint
			end
		end
	end;
end;

procedure CheckResult(aP11, aP12, aP21, aP22: TDPoint; var aP1, aP2: TDPoint; var aResult: Integer);

	function GetZone: Double;
	begin
		Result := Abs(aP12.y - aP11.y) + Abs(aP12.x - aP11.x) + Abs(aP22.y - aP21.y) + Abs(aP22.x - aP21.x);
		Result := Result / 100000000.0
	end;

var
	aZone: Double;
begin
	if aP11.IsEqual(aP21) or aP11.IsEqual(aP22) then begin
		aZone := GetZone;
		if aP11.CloselyTo(aP1, aZone) then
			aP1 := aP11
		else if aP11.CloselyTo(aP2, aZone) then
			aP2 := aP11
		else begin
			case aResult of
				0: 	aP1 := aP11;
				1: 	aP2 := aP11;
				2: 	if aP1.DistanceFrom(aP11) < aP2.DistanceFrom(aP11) then
							aP1 := aP11
						else
							aP2 := aP11
			end;
			Inc(aResult)
		end
	end else
	if aP12.IsEqual(aP21) or aP12.IsEqual(aP22) then begin
		aZone := GetZone;
		if aP12.CloselyTo(aP1, aZone) then
			aP1 := aP12
		else if aP12.CloselyTo(aP2, aZone) then
			aP2 := aP12
		else begin
			case aResult of
				0: 	aP1 := aP12;
				1: 	aP2 := aP12;
				2: 	if aP1.DistanceFrom(aP12) < aP2.DistanceFrom(aP12) then
							aP1 := aP12
						else
							aP2 := aP12
			end;
			Inc(aResult)
		end
	end;
end;

function CheckResult2(const aPt:array of TDPoint; minZone: Double; var pt: TDPoint): Boolean;
var
	dis, min: Double;
	index, i: Integer;
begin
	Result :=False;
	min := aPt[0].DistanceFrom(pt);
	index :=0;
	for i := 1 to 3 do begin
		dis := aPt[i].DistanceFrom(pt);
		if dis < min then begin
			min := dis;
			index := i;
    end;
  end;
  if aPt[index].CloselyTo(pt, minZone) then begin
    pt := aPt[index];
    Result := True;
  end;
end;

procedure CheckResult3(const aP1, aP2, aP3, aP4: TDPoint; numPoints: Integer;
																													var pt1,pt2: TDPoint);
var
	minZone: Double;
begin
	if numPoints < 1 then Exit;
	minZone := dMinZone; //kMinZone;
  CheckResult2([aP1, aP2, aP3, aP4],minZone,pt1);
  if numPoints < 2 then Exit;
	CheckResult2([aP1, aP2, aP3, aP4],minZone,pt2);
	if pt1.IsEqual(pt2) then
		raise Exception.Create('CheckResult3: Точки склеелись');
end;

function ArcContainsPoint(const ptC,pt1: TDPoint; delta: Double;
																				 const ptTest: TDPoint): Boolean;
var
	zone, angle1, angle2, angleTest, tmp: Double;
begin
	zone := dMinZone/(ptC.DistanceFrom(pt1));
	angle1 := GetAngle(ptC, pt1);
	angleTest := GetAngle(ptC, ptTest);
	angle2 := angle1 + delta;
	if delta = 0 then raise Exception.Create('arc: delta = 0');
	if delta < 0 then begin
		tmp := angle1;
		angle1 := angle2;
		angle2 := tmp;
	end;

	if angle1 < 0 then begin
		angle1 := angle1 + 2 * Pi;
		angle2 := angle2 + 2 * Pi;
	end;
	if angle2 < 2 * Pi then begin
		Result := (angleTest >= angle1 - zone ) and (angleTest <= angle2 + zone );
	end
	else
		if angleTest >= angle1 - zone then Result := True
		else begin
			angleTest := angleTest + 2 * Pi;
			Result := (angleTest <= angle2 + zone );
		end;
end;


function IntersectLineArc(const aLine: TDLine; const anArc: TDArc; var aP1, aP2: TDPoint): Integer;
var
	aCenter, aC: TDPoint;
	aRadius, aH, aDelta: Double;
	aCoordSys: TCoordSys;
	pt2Loc: TDPoint;

	function Check(aPoint: TDPoint): Boolean;
	begin
		Result := (aPoint.x >= -kMinZone) and
			(aPoint.x <= aCoordSys.ToLocal(aLine.P2).x + kMinZone);
		if Result then
			Result := ArcContainsPoint(aC, aCoordSys.ToLocal(anArc.P1), aDelta, aPoint);
	end;

var
	aSingular1, aSingular2: Boolean;
begin
	Result := 0;

	aSingular1 := aLine.IsSingular;
	aSingular2 := anArc.IsSingular;
	if aSingular1 and aSingular2 then begin
		if aLine.P1.IsEqual(anArc.P1) then begin
			aP1 := aLine.P1;
			Result := 1;
		end;
	end else if aSingular1 then begin
		if DArcContaintPoint(anArc, aLine.P1, kMinZone) then begin
			aP1 := aLine.P1;
			Result := 1;
		end;
	end else if aSingular2 then begin
		if DSegmentContaintPoint(aLine, anArc.P1, kMinZone) then begin
			aP1 := anArc.P1;
			Result := 1;
		end;
	end else begin
		if GetCenter2(anArc.P1, anArc.P2, anArc.C, aDelta, aCenter) then begin
			aRadius := aCenter.DistanceFrom(anArc.P1);
			aCoordSys.InitPoints(aLine.P1, aLine.P2);
			pt2Loc := aCoordSys.ToLocal(aLine.P2);

			aC := aCoordSys.ToLocal(aCenter);
			aH := Abs(aC.y);
			if abs(aRadius - aH)< kMinZone then begin
				aP1.x := aC.x;
				aP1.y := 0;
				if Check(aP1) then begin
					aCoordSys.MakeGlobal(aP1);
					Inc(Result);
				end;
			end
			else if aRadius >= aH then begin
				aP1.x := aC.x - Sqrt(Sqr(aRadius) - Sqr(aH));
				aP1.y := 0;
				aP2.x := aC.x + Sqrt(Sqr(aRadius) - Sqr(aH));
				aP2.y := 0;
				if Check(aP1) then begin
					aCoordSys.MakeGlobal(aP1);
					Inc(Result);
				end;
				if Check(aP2) and (aRadius > aH) then begin
					aCoordSys.MakeGlobal(aP2);
					if Result = 0 then
						aP1 := aP2;
					Inc(Result)
				end;
			end;
			CheckResult(aLine.P1, aLine.P2, anArc.P1, anArc.P2, aP1, aP2, Result);
			CheckResult3(aLine.P1, aLine.P2, anArc.P1, anArc.P2, Result, aP1, aP2);
		end
	end;
end;

function DPointAdd(const pt1,pt2: TDPoint): TDPoint;
begin
	Result.X := pt1.X + pt2.X;
	Result.Y := pt1.Y + pt2.Y;
end;

function DPointSub(const pt1,pt2: TDPoint): TDPoint;
begin
	Result.X := pt1.X - pt2.X;
	Result.Y := pt1.Y - pt2.Y;
end;

function DPointMult(const pt: TDPoint;m: Double): TDPoint;
begin
	Result.X := pt.X * m;
	Result.Y := pt.Y * m;
end;

function DPointScalar(const pt1,pt2: TDPoint): Double;
begin
	result:= pt1.x*pt2.x + pt1.y*pt2.y;
end;

function DPointModul(const pt: TDPoint): Double;
begin
	result:= sqrt(DPointScalar(pt, pt));
end;

function DPointNorm(const pt: TDPoint): TDPoint;
var
  modul: Double;
begin
  Modul := DPointModul(pt);
  Result.X := pt.X / Modul;
  Result.Y := pt.Y / Modul;
end;


function IntersectArcArc(const anArc1, anArc2: TDArc; var aP1, aP2: TDPoint): Integer;
var
	aCenter1, aCenter2: TDPoint;
	aRadius1, aRadius2, aH, aDelta1, aDelta2, aC1, aDistance: Double;
	aCoordSys: TCoordSys;
	aSingular1, aSingular2: Boolean;

	function Check(aPoint: TDPoint): Boolean;
	begin
//			Result := AngleInRange(GetAngle(aCenter1, aPoint), GetAngle(aCenter1, anArc1.P1), aDelta1) and
//								AngleInRange(GetAngle(aCenter2, aPoint), GetAngle(aCenter2, anArc2.P1), aDelta2);
		Result := ArcContainsPoint(aCenter1,  anArc1.P1, aDelta1 , aPoint) and
							ArcContainsPoint(aCenter2,  anArc2.P1, aDelta2 , aPoint);
	end;

begin
	Result := 0;

	aSingular1 := anArc1.IsSingular;
	aSingular2 := anArc2.IsSingular;

	if aSingular1 and aSingular2 then begin
		if anArc1.P1.IsEqual(anArc2.P1) then begin
			aP1 := anArc1.P1;
			Result := 1;
		end
	end
	else if aSingular1 then begin
		if DArcContaintPoint(anArc2, anArc1.P1) then begin
			aP1 := anArc1.P1;
			Result := 1;
		end;
	end
	else if aSingular2 then begin
		if DArcContaintPoint(anArc1, anArc2.P1) then begin
			aP1 := anArc2.P1;
			Result := 1;
		end;
	end
	else if GetCenter2(anArc1.P1, anArc1.P2, anArc1.C, aDelta1, aCenter1) then begin
		aRadius1 := aCenter1.DistanceFrom(anArc1.P1);
		if GetCenter2(anArc2.P1, anArc2.P2, anArc2.C, aDelta2, aCenter2) then begin
			aRadius2 := aCenter2.DistanceFrom(anArc2.P1);
			aDistance := aCenter1.DistanceFrom(aCenter2);

			if aCenter1.IsEqual(aCenter2) then begin
				if Abs(aRadius1 - aRadius2) < Epsilon then
					Result := -1
				else
					Result := 0
			end
			else if	(abs(aDistance - (aRadius1 + aRadius2)) < kMinZone) or
							(abs(aRadius1 - (aDistance + aRadius2)) < kMinZone) or
							(abs(aRadius2 - (aDistance + aRadius1)) < kMinZone)
			then begin
				if aRadius1 >= aRadius2 then
					aP1 := DPointAdd(aCenter1,
													DPointMult(DPointNorm(DPointSub(aCenter2,aCenter1)), aRadius1))
				else
					aP1 := DPointAdd(aCenter2,
													DPointMult(DPointNorm(DPointSub(aCenter1,aCenter2)), aRadius2));

				if Check(aP1) then
						Inc(Result);
			end
			else if (aDistance <= aRadius1 + aRadius2) and
							(aRadius1 <= aDistance + aRadius2) and
							(aRadius2 <= aDistance + aRadius1)
			then
				if aCenter1.IsEqual(aCenter2) then begin
					aRadius1 := GetAngle(aCenter2, anArc2.P1);
					aRadius2 := GetAngle(aCenter1, anArc1.P1);
					if AngleInRange(GetAngle(aCenter2, anArc1.P1), aRadius1, aDelta2) then begin
						Inc(Result);
						aP1 := anArc1.P1;
					end;
					if AngleInRange(GetAngle(aCenter2, anArc1.P2), aRadius1, aDelta2) then begin
						Inc(Result);
						if Result = 1 then aP1 := anArc1.P2 else aP2 := anArc1.P2
					end;
					if AngleInRange(GetAngle(aCenter1, anArc2.P1), aRadius2, aDelta1) then begin
						Inc(Result);
						if Result = 1 then aP1 := anArc2.P1 else aP2 := anArc2.P1
					end;
					if AngleInRange(GetAngle(aCenter1, anArc2.P2), aRadius2, aDelta1) then begin
						Inc(Result);
						if Result = 1 then aP1 := anArc2.P2 else aP2 := anArc2.P2
					end;
					if Result > 0 then Result := 3
				end else begin
					aC1 := (Sqr(aRadius1) - Sqr(aRadius2) + Sqr(aDistance))/2/aDistance;
					aH := Sqrt(Abs(Sqr(aRadius1) - Sqr(aC1)));
					aP1.x := aC1;
					aP2.x := aC1;
					aP1.y := aH;
					aP2.y := -aH;
					aCoordSys.InitPoints(aCenter1, aCenter2);
					aCoordSys.MakeGlobal(aP1);
					aCoordSys.MakeGlobal(aP2);
					if Check(aP1) then
						Inc(Result);
					if Check(aP2) then begin
						if Result = 0 then
							aP1 := aP2;
						Inc(Result)
					end
				end;
			CheckResult(anArc1.P1, anArc1.P2, anArc2.P1, anArc2.P2, aP1, aP2, Result);
			CheckResult3(anArc1.P1, anArc1.P2, anArc2.P1, anArc2.P2, Result, aP1, aP2);
		end
	end
end;

function IntersectEdgeEdge(const anEdge1, anEdge2: TDEdge; var aP1, aP2: TDPoint; out aTouched: Boolean): Integer;
begin
	if anEdge1.IsLine and anEdge2.IsLine then
		Result := IntersectLineLine(anEdge1, anEdge2, aP1, aP2, aTouched)
	else if anEdge1.IsArc and  anEdge2.IsArc then
		Result := IntersectArcArc(anEdge1, anEdge2, aP1, aP2)
	else if anEdge1.IsLine then
		Result := IntersectLineArc(anEdge1, anEdge2, aP1, aP2)
	else
		Result := IntersectLineArc(anEdge2, anEdge1, aP1, aP2);
end;

function IsSegmentVertex(const aSegment: TDLine; const aP: TDPoint): Boolean;
begin
	Result := aSegment.P1.IsEqual(aP) or aSegment.P2.IsEqual(aP);
end;

function InnerIntersectEdgeEdge(const anEdge1, anEdge2: TDEdge; var aP1, aP2: TDPoint): Integer;
var
	aTouched: Boolean;
begin
	Result := IntersectEdgeEdge(anEdge1, anEdge2, aP1, aP2, aTouched);
	if aTouched and (Result = 2) then begin
		Result := 0;
		Exit;
	end;

	if (Result > 1) and (IsSegmentVertex(anEdge1, aP2) or IsSegmentVertex(anEdge2, aP2)) then
		Dec(Result);

	if (Result > 0) and (IsSegmentVertex(anEdge1, aP1) or IsSegmentVertex(anEdge2, aP1)) then begin
		aP1 := aP2;
		Dec(Result);
	end;
end;

function IntersectLineLineCircle(const aLine: TDLine; const aCenter: TDPoint; aRadius: Double; var aP1, aP2: TDPoint): Integer;
var
	aCoordSys: TCoordSys;
	aLocalCenter: TDPoint;
	x: Double;
begin
	aCoordSys.InitPoints(aLine.P1, aLine.P2);
	aLocalCenter := aCoordSys.ToLocal(aCenter);
	x := Sqr(aRadius) - Sqr(aLocalCenter.y);
	if x > 0 then begin
		if x < Epsilon then begin
			aP1 := aCoordSys.MakeGlobalFromXY(aLocalCenter.x, 0);
			Result := 1;
		end
		else begin
			aP1 := aCoordSys.MakeGlobalFromXY(aLocalCenter.x - x, 0);
			aP2 := aCoordSys.MakeGlobalFromXY(aLocalCenter.x + x, 0);
			Result := 2;
		end;
	end
	else
		Result := 0;
end;

function IntersectCircleCircle(const aCenter1: TDPoint; aRadius1: Double;
	const aCenter2: TDPoint; aRadius2: Double; var aP1, aP2: TDPoint): Integer;
var
	aCoordSys: TCoordSys;

	a, h, d: Double;
begin
	d := aCenter1.DistanceFrom(aCenter2);

	Result := 0;
	if (d <= aRadius1 + aRadius2) and (d >= Abs(aRadius1 - aRadius2)) and (d > 0) then begin
		aCoordSys.InitPoints(aCenter1, aCenter2);

		a := (Sqr(d) + Sqr(aRadius1) - Sqr(aRadius2))/(2*d);
		h := Sqrt(Sqr(aRadius1) - Sqr(a));
		if h < Epsilon then begin
			aP1 := aCoordSys.MakeGlobalFromXY(a, 0);
			Result := 1;
		end else begin
			aP1 := aCoordSys.MakeGlobalFromXY(a, h);
			aP2 := aCoordSys.MakeGlobalFromXY(a, -h);
			Result := 2;
		end;
	end
	else if (d < Epsilon/2) and  (Abs(aRadius1 - aRadius2) < Epsilon/2) then
		Result := -1;
end;

function IntersectLineLineSegment(const aLine, aSegment: TDLine; var aP: TDPoint): Integer;
var
	aCoordSys: TCoordSys;
	aP1, aP2: TDPoint;
	DY: Double;
begin
	aCoordSys.InitPoints(aLine.P1, aLine.P2);

	aP1 := aCoordSys.ToLocal(aSegment.P1);
	aP2 := aCoordSys.ToLocal(aSegment.P2);
	DY := aP2.y - aP1.y;
	if ((aP1.y >= 0) and (aP2.y <= 0)) or ((aP1.y <= 0) and (aP2.y >= 0)) then begin
		if Abs(DY) < Epsilon then
			Result := -1
		else begin
			aP := aCoordSys.MakeGlobalFromXY(aP1.x - aP1.y*(aP2.x - aP1.x)/DY, 0);
			Result := 1;
		end;
	end
	else
		Result := 0;
end;

function DPoint(aX, aY: Double): TDPoint;
begin
	Result.x := aX;
	Result.y := aY;
end;

function DLine(aP1, aP2: TDPoint): TDLine;
begin
	Result.P1 := aP1;
	Result.P2 := aP2
end;

function DArc(aP1, aP2: TDPoint; aC: Double): TDArc;
begin
	Result.P1 := aP1;
	Result.P2 := aP2;
	Result.C := aC
end;

function DEdge(const aP1, aP2: TDPoint; aC: Double): TDEdge;
begin
	with Result do begin
		P1 := aP1;
		P2 := aP2;
		C := aC;
	end;
end;

function DRect(aX1, aY1, aX2, aY2: Double): TDRect;
begin
	Result.x1 := aX1;
	Result.y1 := aY1;
	Result.x2 := aX2;
	Result.y2 := aY2;
end;

{	TDPoint }

procedure TDPoint.Assign(aX, aY: Double);
begin
	x := aX;
	y := aY;
end;

procedure TDPoint.Offset(aDx, aDy: Double);
begin
	x := x + aDx;
	y := y + aDy;
end;

function TDPoint.IsEqual(const aP: TDPointBase): Boolean;
begin
	Result := (aP.x = x) and (aP.y = y);
end;

function TDPoint.DistanceFrom(const aP: TDPointBase): Double;
begin
	Result := Sqrt(Sqr(x - aP.x) + Sqr(y - aP.y))
{		Result := Hypot(x - aP.x, y - aP.y);
}	end;

function TDPoint.CloselyTo(const aP: TDPointBase; aDelta: Double): Boolean;
begin
	Result := (abs(x - aP.x) <= aDelta) and (abs(y - aP.y) <= aDelta)
end;

function TDPoint.DistanceFromSegment(const aP1, aP2: TDPointBase): Double;
var
	aCoordSys: TCoordSys;
	aPoint: TDPoint;
begin
	if TDPoint(aP1).IsEqual(aP2) then
		Result := DistanceFrom(aP1)
	else begin
		aCoordSys.InitPoints(aP1, aP2);
		aPoint := aCoordSys.ToLocal(Self);
		if aPoint.X < 0 then
			Result := aPoint.DistanceFrom(DPoint(0, 0))
		else if aPoint.X > aCoordSys.ToLocal(aP2).X then
			Result := aPoint.DistanceFrom(aCoordSys.ToLocal(aP2))
		else
			Result := abs(aPoint.Y);
	end
end;

function TDPoint.DistanceFromEdge(const aP1, aP2: TDPoint;
	aConvexity: Double = 0): Double;
var
	aDelta: Double;
	aCenter: TDPoint;
begin
	if Abs(aConvexity) < minConvexity2 then
		Result := DistanceFromSegment(aP1, aP2)
	else begin
		GetCenter2(aP1, aP2, aConvexity, aDelta, aCenter);
		Result := Min(Min(DistanceFrom(aP1), DistanceFrom(aP2)),
			Abs(aCenter.DistanceFrom(Self) - aCenter.DistanceFrom(aP1)));
	end;
end;

function TDPoint.DistanceFromLine(const aP1, aP2: TDPointBase): Double;
var
	Dx: Extended;
	Dy: Extended;
begin
	Dx := aP2.x - aP1.x;
	Dy := aP2.y - aP1.y;

	if (Dx = 0) and (Dy = 0) then
		raise Exception.CreateFmt('Singular line (%f, %f) (TDPoint.DistanceFromLine)',
			[aP1.X, aP1.Y]);

	Result := Abs((-(x - aP1.x)*DY  + (y - aP1.y)*DX) / Sqrt(Sqr(Dy) + Sqr(Dx)));
end;


function DPointOffsetBy(const aP: TDPoint; aDX, aDY: Double): TDPoint;
begin
	Result.x := aP.x + aDX;
	Result.y := aP.y + aDY;
end;

{	TDRect }

procedure TDRect.Assign(aX1, aY1, aX2, aY2: Double);
begin
	x1 := aX1;
	y1 := aY1;
	x2 := aX2;
	y2 := aY2;
end;

procedure TDRect.AssignPoints(const aP1, aP2: TDPoint);
begin
	x1 := aP1.x;
	y1 := aP1.y;
	x2 := aP2.x;
	y2 := aP2.y;
	Normalized;
end;

procedure TDRect.AssignFrom(const aPoints: array of TDPoint);
var
	i: Integer;
begin
	AssignPoints(aPoints[Low(aPoints)], aPoints[Low(aPoints)]);
	for i := Low(aPoints) + 1 to High(aPoints) do
		Expand(aPoints[i]);
end;

function TDRect.DX: Double;
begin
	Result := x2 - x1;
end;

function TDRect.DY: Double;
begin
	Result := y2 - y1;
end;

function DRectIsEqual(const aR1, aR2: TDRect): Boolean;
begin
	with aR1 do
		Result := (x1 = aR2.x1) and (y1 = aR2.y1) and
							(x2 = aR2.x2) and (y2 = aR2.y2);
end;

procedure DRectIntersectWith(var aR1: TDRect; const aR2: TDRect);
var
	d: Double;
begin
	with aR1 do begin
		d := Max(x1, aR2.x1);
		x2 := Min(x2, aR2.x2);
		x1 := d;

		d := Max(y1, aR2.y1);
		y2 := Min(y2, aR2.y2);
		y1 := d;
	end;
end;

procedure DRectUnionWith(var aR1: TDRect; const aR2: TDRect);
var
	d: Double;
begin
	with aR1 do begin
		d := Min(x1, aR2.x1);
		x2 := Max(x2, aR2.x2);
		x1 := d;

		d := Min(y1, aR2.y1);
		y2 := Max(y2, aR2.y2);
		y1 := d;
	end;
end;

procedure DRectExpandByRect(var aR1: TDRect; const aR2: TDRect);
begin
	with aR1 do begin
		if aR2.x1 < x1 then x1 := aR2.x1;
		if aR2.x2 > x2 then x2 := aR2.x2;
		if aR2.y1 < y1 then y1 := aR2.y1;
		if aR2.y2 > y2 then y2 := aR2.y2;
	end;
end;

procedure DRectExpandByExtentions(var aR1: TDRect; const aR2: TDRect);
begin
	aR1.x1 := aR1.x1 + aR2.x1;
	aR1.x2 := aR1.x2 + aR2.x2;
	aR1.y1 := aR1.y1 + aR2.y1;
	aR1.y2 := aR1.y2 + aR2.y2;
end;

function DRectIntersectedWith(const aR1, aR2: TDRect): Boolean;
begin
	Result := False;
	with aR1 do begin
		if x1 > aR2.x2 then Exit;
		if y1 > aR2.y2 then Exit;
		if x2 < aR2.x1 then Exit;
		if y2 < aR2.y1 then Exit;
	end;
	Result := True;
end;

function DRectIntersectedWith(aRect: TDRect; const anEdge: TDEdge): Boolean;

	function IsIntercectedWithLine: Boolean;
	var
		P, K: Double;
	begin
		Result := False;
		if Abs(anEdge.P2.y - anEdge.P1.y) > Epsilon then begin
			K := (anEdge.P2.x - anEdge.P1.x)/(anEdge.P2.y - anEdge.P1.y);

			P := k*(aRect.y1 - anEdge.P1.y) + anEdge.P1.x;
			if (aRect.x1 <= P) and (P <= aRect.x2) then begin
				Result := True;
				Exit;
			end;

			P := k*(aRect.y2 - anEdge.P1.y) + anEdge.P1.x;
			if (aRect.x1 <= P) and (P <= aRect.x2) then begin
				Result := True;
				Exit;
			end;
		end;

		if Abs(anEdge.P2.x - anEdge.P1.x) > Epsilon then begin
			K := (anEdge.P2.y - anEdge.P1.y)/(anEdge.P2.x - anEdge.P1.x);

			P := k*(aRect.x1 - anEdge.P1.x) + anEdge.P1.y;
			if (aRect.y1 <= P) and (P <= aRect.y2) then begin
				Result := True;
				Exit;
			end;

			P := k*(aRect.x2 - anEdge.P1.x) + anEdge.P1.y;
			if (aRect.y1 <= P) and (P <= aRect.y2) then begin
				Result := True;
				Exit;
			end;
		end
	end;

	function IsIntercectedWithArc: Boolean;
	var
		aCenter: TDPoint;
		aRadius: Double;
		aDelta: Double;

		P, r2, d: Double;
	begin
		Result := False;
		if GetCenter2(anEdge.P1, anEdge.P2, anEdge.C, aDelta, aCenter) then begin
			aRadius := aCenter.DistanceFrom(anEdge.P1);

			r2 := Sqr(aRadius);
			if Abs(aCenter.x  - aRect.x1) <= aRadius then begin
				d := Sqrt(r2 - Sqr(aRect.x1 - aCenter.x));

				P := d + aCenter.y;
				if (aRect.y1 <= P) and (P <= aRect.y2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(aRect.x1, P))
				then begin
					Result := True;
					Exit;
				end;

				P := -d + aCenter.y;
				if (aRect.y1 <= P) and (P <= aRect.y2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(aRect.x1, P))
				then begin
					Result := True;
					Exit;
				end;
			end;

			if Abs(aCenter.x  - aRect.x2) <= aRadius then begin
				d := Sqrt(r2 - Sqr(aRect.x2 - aCenter.x));

				P := d + aCenter.y;
				if (aRect.y1 <= P) and (P <= aRect.y2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(aRect.x2, P))
				then begin
					Result := True;
					Exit;
				end;

				P := -d + aCenter.y;
				if (aRect.y1 <= P) and (P <= aRect.y2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(aRect.x2, P))
				then begin
					Result := True;
					Exit;
				end;
			end;


			if Abs(aCenter.y  - aRect.y1) <= aRadius then begin
				d := Sqrt(r2 - Sqr(aRect.y1 - aCenter.y));

				P := d + aCenter.x;
				if (aRect.x1 <= P) and (P <= aRect.x2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(P, aRect.y1))
				then begin
					Result := True;
					Exit;
				end;

				P := -d + aCenter.y;
				if (aRect.x1 <= P) and (P <= aRect.x2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(P, aRect.y1))
				then begin
					Result := True;
					Exit;
				end;
			end;

			if Abs(aCenter.y  - aRect.y2) <= aRadius then begin
				d := Sqrt(r2 - Sqr(aRect.y2 - aCenter.y));

				P := d + aCenter.x;
				if (aRect.x1 <= P) and (P <= aRect.x2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(P, aRect.y2))
				then begin
					Result := True;
					Exit;
				end;

				P := -d + aCenter.y;
				if (aRect.x1 <= P) and (P <= aRect.x2) and
					ArcContainsPoint(aCenter, anEdge.P1, aDelta, DPoint(P, aRect.y2))
				then begin
					Result := True;
					Exit;
				end;
			end;
		end;
	end;

begin
	Result := False;
	if aRect.Contains(anEdge.P1) or aRect.Contains(anEdge.P1) then
		Result := True
	else if not anEdge.IsSingular then begin
		aRect.Normalized;
		if anEdge.IsLine then
			Result := IsIntercectedWithLine
		else
			Result := IsIntercectedWithArc;
	end
end;


function LinePointPosition(aLine: TDLine; aPoint: TDPoint): Integer;
var
	DX, DY: double;
	Delta: Double;
begin
	DX := aLine.P2.x - aLine.P1.x;
	DY := aLine.P2.y - aLine.P1.y;

	Assert(not((DX = 0) and (DY = 0)));

	if DX >= DY then
		Delta := (DY/DX)*(aPoint.x - aLine.P1.x) - (aPoint.y - aLine.P1.x)
	else
		Delta := (aPoint.x - aLine.P1.x) - (DX/DY)*(aPoint.y - aLine.P1.x);

	Result := Sign(Delta);
end;


function DRectContainsRect(const aR1, aR2: TDRect): Boolean;
begin
	with aR1 do
		Result :=
			(aR2.x1 >= x1) and (aR2.x1 <= x2) and
			(aR2.x2 >= x1) and (aR2.x2 <= x2) and
			(aR2.y1 >= y1) and (aR2.y1 <= y2) and
			(aR2.y2 >= y1) and (aR2.y2 <= y2);
end;

function TDRect.Center: TDPoint;
begin
	Result.Assign((x1 + x2) / 2, (y1 + y2) / 2);
end;

procedure TDRect.Offset(aDx, aDy: double);
begin
	x1 := x1 + aDx;
	x2 := x2 + aDx;
	y1 := y1 + aDy;
	y2 := y2 + aDy;
end;

procedure TDRect.Inflate(aDx, aDy: double);
begin
	x1 := x1 - aDx;
	x2 := x2 + aDx;
	y1 := y1 - aDy;
	y2 := y2 + aDy;
end;

function TDRect.IsNormalized: Boolean;
begin
	Result := (x1 < x2 ) and (y1 < y2);
end;

procedure TDRect.Expand(const aPoint: TDPoint);
begin
	if aPoint.x < x1 then x1 := aPoint.x;
	if aPoint.x > x2 then x2 := aPoint.x;
	if aPoint.y < y1 then y1 := aPoint.y;
	if aPoint.y > y2 then y2 := aPoint.y;
end;

procedure TDRect.Clear;
begin
	FillChar(Self, sizeof(Self), 0);
end;

function TDRect.Contains(const aPoint: TDPoint): Boolean;
begin
	Result := (aPoint.x >= x1) and (aPoint.x <= x2) and
						(aPoint.y >= y1) and (aPoint.y <= y2);
end;

procedure TDRect.Normalized;
var
	tmp: Double;
begin
	if x1 > x2 then begin
		tmp := x1;
		x1 := x2;
		x2 := tmp;
	end;
	if y1 > y2 then begin
		tmp := y1;
		y1 := y2;
		y2 := tmp;
	end;
end;

function DPointFromPoint(aPoint: TPoint): TDPoint;
begin
	Result.X := aPoint.X;
	Result.Y := aPoint.Y;
end;

function DRectFromRect(const aRect: TRect): TDRect;
begin
	Result.X1 := aRect.Left;
	Result.Y1 := aRect.Top;
	Result.X2 := aRect.Right;
	Result.Y2 := aRect.Bottom;
end;

function RectFromDRect(const aR: TDRect): TRect;
begin
	Result.Left := Round(aR.X1);
	Result.Top := Round(aR.Y1);
	Result.Right := Round(aR.X2);
	Result.Bottom := Round(aR.Y2);
end;

function PointFromDPoint(const aDPoint: TDPoint): TPoint;
begin
	Result.x := Round(aDPoint.x);
	Result.y := Round(aDPoint.y);
end;
	

{ TDPointList }

function TDPointList.GetPoint(I: Integer): TDPoint;
begin
	Result := PDPoint(List.Items[I])^;
end;

procedure TDPointList.SetPoint(I: Integer; const aValue: TDPoint);
begin
	PDPoint(List.Items[I])^ := aValue;
end;

function TDPointList.GetCount: Integer;
begin
	Result := List.Count;
end;

procedure TDPointList.AddPoint(const aPoint: TDPoint);
begin
	List.Add(aPoint);
end;

procedure TDPointList.AddPoints(aPointList: TDPointList);
var
	i: Integer;
begin
	for i := 0 to aPointList.Count - 1 do
		AddPoint(aPointList[i]);
end;

procedure TDPointList.InsertPoint(anIndex: Integer; const aPoint: TDPoint);
begin
	List.Insert(anIndex, aPoint);
end;

procedure TDPointList.DeletePoint(Index: Integer);
begin
	List.Delete(Index);
end;

procedure TDPointList.LoadFromStream(aStream: TStream);
begin
	List.LoadFromStream(aStream);
	if DataVersion >= $0101 then
		aStream.Read(Closed, sizeof(Closed))
	else
		Closed := False;
end;

procedure TDPointList.SaveToStream(aStream: TStream);
begin
	List.SaveToStream(aStream);
	aStream.Write(Closed, sizeof(Closed));
end;

constructor TDPointList.Create;
begin
	inherited Create;
	List := TRecList.Create(sizeof(TDPoint));
end;

destructor TDPointList.Destroy;
begin
	List.Free;
	inherited Destroy;
end;

procedure TDPointList.Assign(aSrc: TDPointList);
begin
	List.Assign(aSrc.List);
	Closed := aSrc.Closed;
end;

procedure TDPointList.Clear;
begin
	List.Clear;
end;

procedure TDPointList.AssignPoints(aPoints: PDPoint; aCount: Integer);
begin
	List.LoadFromMemory(aPoints, aCount);
end;

function TDPointList.ContainsPoint(aPoint: TDPoint): Boolean;
begin
	if (List.Count < 3) or not Closed then begin
		Result := False;
		Exit
	end;
	Result := Odd(GetIntersectionCount(aPoint));
end;

function TDPointList.GetIntersectionCount(aPoint: TDPoint): Integer;
var
	Intersected: Boolean;
	IntersectionX: Double;
	I: Integer;
begin
	Result := 0;
	for I := 1 to List.Count - 1 do
	begin
		Intersected := GetXDIntersection(Points[I - 1], Points[I], aPoint, IntersectionX);
		if Intersected and (IntersectionX <= aPoint.X) then
			Inc(Result);
	end;
	Intersected := GetXDIntersection(Points[List.Count - 1], Points[0], aPoint, IntersectionX);
	if Intersected and (IntersectionX <= aPoint.X) then
		Inc(Result);
end;

function TDPointList.EdgeAtPoint(const aPoint: TDPoint; aDistance: Double): Integer;
var
	i: Integer;
	aFirst,
	aPrev,
	aNext: PDPoint;
begin
	if List.Count > 1 then begin
		aFirst := List.Items[0];
		aPrev := aFirst;
		for i := 1 to List.Count - 1 do begin
			aNext := List.Items[i];
			if IsPointNearEdge(aPoint, aPrev^, aNext^, aDistance) then begin
				Result := i - 1;
				Exit;
			end;
			aPrev := aNext;
		end;
		if Closed then
			if IsPointNearEdge(aPoint, aPrev^, aFirst^, aDistance) then begin
				Result := List.Count - 1;
				Exit;
			end;
	end;
	Result := -1;
end;

function TDPointList.GetSquare: Double;
var
	P1, P2: PDPoint;
	I: Integer;
begin
	Result := 0;
	if List.Count > 2 then
	begin
		for I := 0 to List.Count - 2 do
		begin
			P2 := List.Items[I + 1];
			P1 := List.Items[I];
			Result := Result + (P2^.X - P1^.X) * (P2^.Y + P1^.Y) / 2;
		end;

		P2 := List.Items[0];
		P1 := List.Items[List.Count - 1];
		Result := Result + (P2^.X - P1^.X) * (P2^.Y + P1^.Y) / 2;
	end;
end;

function TDPointList.LineTouched(aP1, aP2: TDPoint): Boolean;
var
	i: Integer;
	aPoint: TDPoint;
begin
	Result := True;
	for i := 0 to Count - 2 do
		if Intersect(aP1, aP2, Points[i], Points[i + 1], aPoint) then Exit;
	if Closed then
		if Intersect(aP1, aP2, Points[0], Points[Count - 1], aPoint) then Exit;
	Result := False;
end;

function TDPointList.GetDistanceFromPoint(aPoint: TDPoint): Double;
var
	aDistance: Double;
	i: Integer;
begin
	Result := 0;
	if Count < 1 then Exit;
	Result := aPoint.DistanceFrom(Points[0]);
	for i := 0 to Count - 2 do begin
		aDistance := aPoint.DistanceFromSegment(Points[i], Points[i + 1]);
		if aDistance < Result then
			Result := aDistance
	end;
	if Closed then begin
		aDistance := aPoint.DistanceFromSegment(Points[0], Points[Count - 1]);
		if aDistance < Result then
			Result := aDistance
	end;
end;

function TDPointList.GetBounds(var aBounds: TDRect): Boolean;
var
	i: Integer;
	aP: TDPoint;
begin
	Result := False;
	for i := 0 to Count - 1 do begin
		aP := Points[i];
		if Result then
			aBounds.Expand(aP)
		else begin
			Result := True;
			aBounds.Assign(aP.x, aP.y, aP.x, aP.y);
		end;
	end;
end;

function TDPointList.GetPrevIndex(anIndex: Integer): Integer;
begin
	if anIndex = -1 then
		Result := -1
	else if anIndex > 0 then
		Result := anIndex - 1
	else if Closed then
		Result := Count - 1
	else
		Result := -1;
end;

function TDPointList.GetNextIndex(anIndex: Integer): Integer;
begin
	if anIndex = -1 then
		Result := -1
	else if anIndex < Count - 1 then
		Result := anIndex + 1
	else if Closed then
		Result := 0
	else
		Result := -1;
end;

procedure TDPointList.DeleteSingularSegments;
var
	aPrev: TDPoint;
	aNext: TDPoint;
	i, j: Integer;
begin
	if Count < 1 then
		Exit;

	if Closed then begin
		aPrev := Points[0];
		i := Count - 1;
		j := 1;
	end
	else begin
		aPrev := Points[Count - 1];
		i := Count - 2;
		j := 0;
	end;

	while i >= j do begin
		aNext := Points[i];
		if aNext.IsEqual(aPrev) then
			DeletePoint(i)
		else
			aPrev := aNext;

		Dec(i);
	end
end;

{ TDPointArray }

procedure TDPointArray.Clear(aCapacity: Integer = 0);
begin
	SetLength(FItems, Max(aCapacity, 0));
	FCount := 0;
end;

function TDPointArray.Add(const aPoint: TDPoint): Integer;
begin
	if FCount >= Length(FItems) then
		SetLength(FItems, Max(1, Length(FItems))*2);

	Result := FCount;
	FItems[FCount] := aPoint;
	Inc(FCount);
end;

procedure TDPointArray.Add(const aPoints: TDPointArray);
var
	i, j: Integer;
begin
	if aPoints.Count = 0 then
		Exit;

	j := Count;

	SetLength(FItems, Count + aPoints.Count);
	for i := 0 to aPoints.Count - 1 do
		FItems[i + j] := aPoints.FItems[i];

	Inc(FCount, aPoints.Count)
end;

procedure TDPointArray.Sort(aCompareFunction: TDPointsCompare; aData: Pointer);

	procedure QuickSort(L, R: Integer);
	var
		x, w: TDPoint;
		i, j: Integer;
	begin
		i := L;
		j := R;
		x := FItems[(i+j) div 2];

		repeat
			while aCompareFunction(FItems[i], x, aData) < 0 do inc(i);
			while aCompareFunction(FItems[j], x, aData) > 0 do dec(j);

			if i <= j then begin
				w := FItems[i];
				FItems[i] := FItems[j];
				FItems[j] := w;

				inc(i);
				dec(j);
			end;
		until i > j;

		if L < j then
			QuickSort(L, j);

		if i < R then
			QuickSort(i, R);
	end;

begin
	if (FCount > 0) then
		QuickSort(0, FCount - 1);
end;

function TDPointArray.GetItem(Index: Integer): TDPoint;
begin
	CheckBounds(Index, FCount);
	Result := FItems[Index];
end;

function TDPointArray.Count: Integer;
begin
	Result := FCount;
end;


{ TCoordSys }

procedure TCoordSys.InitPoints(const aCenter, anEnd: TDPointBase);
var
	aHypot2: Extended;
begin
	Center := aCenter;
	EndPoint := anEnd;

	FDx := anEnd.x - aCenter.x;
	FDy := anEnd.y - aCenter.y;
	aHypot2 := Sqr(FDy) + Sqr(FDx);
	Assert(aHypot2 > 0);
	if aHypot2 = 0 then
		raise Exception.Create('Invalid params (TCoordSys.InitPoints)');

{	if aHypot2 = 0 then begin
		FInitByPoints := False;
		AngleValue := ArcTan2(FDy, FDx);
		SinCos(AngleValue, AngleSin, AngleCos);
		AngleCalculated := True
	end
	else begin}
		FInitByPoints := True;
		FHypot := Sqrt(aHypot2);
		AngleSin := FDY / FHypot;
		AngleCos := FDX / FHypot;
		AngleCalculated := False
//	end;
end;

procedure TCoordSys.InitAngle(const aCenter: TDPoint; anAngle: Extended);
begin
	FInitByPoints := False;
	Center := aCenter;
	AngleValue := anAngle;
	AngleCalculated := True;
	SinCos(anAngle, AngleSin, AngleCos);
end;

procedure TCoordSys.MakeGlobal(var p: TDPoint);
begin
	p := ToGlobal(p);
end;

function TCoordSys.MakeGlobalFromXY(x, y: double): TDPoint;
begin
	if FInitByPoints then begin
		Result.x := (x * FDX  - y * FDY)/FHypot + Center.x;
		Result.y := (x * FDY  + y * FDX)/FHypot + Center.y;
	end
	else begin
		Result.x := x * AngleCos - y * AngleSin + Center.x;
		Result.y := x * AngleSin + y * AngleCos + Center.y;
	end;
end;

procedure TCoordSys.MakeLocal(var p: TDPoint);
begin
	p := ToLocal(p);
end;

function TCoordSys.ToLocal(const p: TDPointBase): TDPoint;
var
	x, y: double;
begin
	x := p.x - Center.x;
	y := p.y - Center.y;

	if FInitByPoints then begin
		Result.x := (x * FDX + y * FDY)/FHypot;
		Result.y := (-x * FDY + y * FDX)/FHypot;
	end
	else begin
		Result.x := x * AngleCos + y * AngleSin;
		Result.y := - x * AngleSin + y * AngleCos;
	end;
end;

function TCoordSys.GetXProjection(const aGlobalP: TDPoint): Extended;
begin
	if FInitByPoints then
		Result := ((aGlobalP.x - Center.x) * FDX + (aGlobalP.y - Center.y) * FDY)/FHypot
	else
		Result := (aGlobalP.x - Center.x) * AngleCos + (aGlobalP.y - Center.y) * AngleSin;
end;

function TCoordSys.GetYProjection(const aGlobalP: TDPoint): Extended;
begin
	if FInitByPoints then
		Result := (-(aGlobalP.x - Center.x)*FDY + (aGlobalP.y - Center.y)*FDX)/FHypot
	else
		Result := -(aGlobalP.x - Center.x)*AngleSin + (aGlobalP.y - Center.y)*AngleCos;
end;

function TCoordSys.ToGlobal(const p: TDPoint): TDPoint;
begin
	with p do begin
		if FInitByPoints then begin
			Result.x := (x * FDX - y * FDY)/FHypot + Center.x;
			Result.y := (x * FDY + y * FDX)/FHypot + Center.y;
		end
		else begin
			Result.x := x * AngleCos - y * AngleSin + Center.x;
			Result.y := x * AngleSin + y * AngleCos + Center.y;
		end;
	end
end;

function TCoordSys.IsValid: Boolean;
begin
	Result := (AngleSin >= -1) and (AngleSin <= 1) and (AngleCos >= -1)
		and (AngleCos <= 1);
end;

procedure TCoordSys.SetInvalid;
begin
	AngleSin := 2;
	AngleCos := 2
end;

function TCoordSys.Angle: Extended;
begin
	if AngleCalculated then
		Result := AngleValue
	else begin
		AngleValue := ArcTan2(FDY, FDX);
		AngleCalculated := True;
		Result := AngleValue
	end
end;

{ geometric functions }

procedure NormalizeLineX(var P1, P2: TDPoint);
var
	Tmp: TDPoint;
begin
	if P1.X > P2.X then begin
		Tmp := P1;
		P1 := P2;
		P2 := Tmp;
	end;
end;

procedure NormalizeLineY(var P1, P2: TDPoint);
var
	Tmp: TDPoint;
begin
	if P1.Y > P2.Y then begin
		Tmp := P1;
		P1 := P2;
		P2 := Tmp;
	end;
end;

function IsPointNearEdge(const p, a, b: TDPoint;aDistance: Double): Boolean;
var
	aCoordSys: TCoordSys;
	aLocalP,
	aLocalB: TDPoint;
begin
	if a.IsEqual(b) then
		Result := p.DistanceFrom(a) < aDistance
	else begin
		aCoordSys.InitPoints(a, b);
		aLocalP := aCoordSys.ToLocal(p);
		aLocalB := aCoordSys.ToLocal(b);
		Result := (Abs(aLocalP.y) < aDistance) and
							(aLocalP.x >= -aDistance) and (aLocalP.x <= aLocalB.x + aDistance);
	end;
end;

function GetXDIntersection(P1, P2: TDPoint; P: TDPoint; var X: Double): Boolean;
begin
	Result := False;
	if not P.IsEqual(P2) and (P1.Y <> P2.Y) then begin
		NormalizeLineY(P1, P2);
		if (P.Y >= P1.Y) and (P.Y <= P2.Y) then begin
			X := (P.Y - P1.Y) / (P2.Y - P1.Y) * (P2.X - P1.X) + P1.X;
			Result := True;
		end;
	end;
end;

procedure MakeOutlinePolygon(aSrc, aDst: TDPointList; aDelta: Double);
var
	aCoordSys, aCoordSys1: TCoordSys;
	p, p1, p2, p3: TDPoint;
	dx, dy: Double;
	i: Integer;
begin
	if aSrc.Count = 0 then Exit;

	if aSrc.Count = 1 then
	begin
		p := aSrc[0];
		aDst.AddPoint(DPoint(p.x - aDelta, p.y - aDelta));
		aDst.AddPoint(DPoint(p.x + aDelta, p.y - aDelta));
		aDst.AddPoint(DPoint(p.x + aDelta, p.y + aDelta));
		aDst.AddPoint(DPoint(p.x - aDelta, p.y + aDelta));
		Exit;
	end;

	if aSrc.Closed then
		aSrc.AddPoint(aSrc[0]);

	p1 := aSrc[0];
	p2 := aSrc[1];

	{	Поместить первую точку }
	aCoordSys.InitPoints(p1, p2);
	p.x := 0;
	p.y := aDelta;
	aCoordSys.MakeGlobal(p);
	aDst.AddPoint(p);

	{	Поместить промежуточные точки }
	for i := 2 to aSrc.Count - 1 do begin
		p1 := aSrc[i - 2];
		p2 := aSrc[i - 1];
		p3 := aSrc[i];
		aCoordSys.InitPoints(p1, p2);
		aCoordSys.MakeLocal(p2);
		aCoordSys.MakeLocal(p3);

		aCoordSys1.InitPoints(p2, p3);

		dx := p3.x-p2.x;
		dy := p3.y-p2.y;
		if dy <> 0 then begin
			p.x := p2.x + (aDelta*(1-aCoordSys1.AngleCos)-p2.y)*dx/dy - aDelta*aCoordSys1.AngleSin;
			p.y := aDelta;
			aCoordSys.MakeGlobal(p);
			aDst.AddPoint(p);
		end;
	end;

	{	Поместить последнюю точку }
	p := aSrc[aSrc.Count - 1];
	aCoordSys.InitPoints(aSrc[aSrc.Count - 2], p);
	aCoordSys.MakeLocal(p);
	p.y := aDelta;
	aCoordSys.MakeGlobal(p);
	aDst.AddPoint(p);

	{--------------------------------------------------}
	p1 := aSrc[aSrc.Count - 1];
	p2 := aSrc[aSrc.Count - 2];

	{	Поместить первую точку }
	aCoordSys.InitPoints(p1, p2);
	p.x := 0;
	p.y := aDelta;
	aCoordSys.MakeGlobal(p);
	aDst.AddPoint(p);

	{	Поместить промежуточные точки }
	for i := aSrc.Count - 3 downto 0 do begin
		p1 := aSrc[i + 2];
		p2 := aSrc[i + 1];
		p3 := aSrc[i];
		aCoordSys.InitPoints(p1, p2);
		aCoordSys.MakeLocal(p2);
		aCoordSys.MakeLocal(p3);

		aCoordSys1.InitPoints(p2, p3);
		dx := p3.x-p2.x;
		dy := p3.y-p2.y;
		if dy <> 0 then begin
			p.x := p2.x + (aDelta*(1-aCoordSys1.AngleCos)-p2.y)*dx/dy - aDelta*aCoordSys1.AngleSin;
			p.y := aDelta;
			aCoordSys.MakeGlobal(p);
			aDst.AddPoint(p);
		end;
	end;

	{	Поместить последнюю точку }
	p := aSrc[0];
	aCoordSys.InitPoints(aSrc[1], p);
	aCoordSys.MakeLocal(p);
	p.y := aDelta;
	aCoordSys.MakeGlobal(p);
	aDst.AddPoint(p);

	if aSrc.Closed then
		aSrc.DeletePoint(aSrc.Count - 1);
	aDst.Closed := True
end;

function ZoomScaleInRange(aZoomScale, aMinScale, aMaxScale: Double): Boolean;
begin
	Result :=
		(aZoomScale = 0) or
		(((aMinScale = 0) or
		 ((aMinScale > 0) and (aZoomScale <= aMinScale)) or
		 ((aMinScale < 0) and (aZoomScale < -aMinScale)))
			and
		((aMaxScale = 0) or
		 ((aMaxScale > 0) and (aZoomScale >= aMaxScale)) or
		 ((aMaxScale < 0) and (aZoomScale > -aMaxScale))))
end;

function SetDPoint(const anOldValue, aNewValue: TDPoint; var aHasChanges: Boolean): TDPoint;
begin
	Result := aNewValue;
	if (aNewValue.X <> anOldValue.X) or (aNewValue.Y <> anOldValue.Y) then
		aHasChanges := True;
end;

function SetDRect(const anOldValue, aNewValue: TDRect; var aHasChanges: Boolean): TDRect;
begin
	Result := aNewValue;
	if (aNewValue.X1 <> anOldValue.X1) or (aNewValue.Y1 <> anOldValue.Y1) or
		(aNewValue.X2 <> anOldValue.X2) or (aNewValue.Y2 <> anOldValue.Y2)
	then
		aHasChanges := True;
end;

function DEdgeIsEqual(const anEdge1, anEdge2: TDEdge): Boolean;
begin
	Result := Abs(Abs(anEdge1.C) - Abs(anEdge2.C)) < minConvexity2;

	if Result then begin
		Result := False;
		if anEdge1.P1.IsEqual(anEdge2.P1) then
			Result := (Sign(anEdge1.C) = Sign(anEdge2.C)) and
				anEdge1.P2.IsEqual(anEdge2.P2)
		else if anEdge1.P1.IsEqual(anEdge2.P2) then
			Result := (Sign(anEdge1.C) = - Sign(anEdge2.C)) and
				anEdge1.P2.IsEqual(anEdge2.P1)
	end;
end;

function DEdgeReverse(const anEdge: TDEdge): TDEdge;
begin
	Result := DEdge(anEdge.P2, anEdge.P1, -anEdge.C);
end;

procedure ClipLineByRect(const aRect: TDRect; const aLine: TDLine; out aClipLine: TDLine);
var
	aCoordSys: TCoordSys;
	aX1, aX2: Double;
	aTestX: Double;
begin
	aCoordSys.InitPoints(aLine.P1, aLine.P2);

	aX1	:= 0;
	aX2 := aCoordSys.GetXProjection(aLine.P2);

	aTestX := aCoordSys.GetXProjection(DPoint(aRect.x1, 0));
	aX1 := Max(aTestX, aX1);
	aX2 := Min(aTestX, aX2);

	aTestX := aCoordSys.GetXProjection(DPoint(aRect.x2, 0));
	aX1 := Max(aTestX, aX1);
	aX2 := Min(aTestX, aX2);

	aTestX := aCoordSys.GetXProjection(DPoint(0, aRect.y1));
	aX1 := Max(aTestX, aX1);
	aX2 := Min(aTestX, aX2);

	aTestX := aCoordSys.GetXProjection(DPoint(0, aRect.y2));
	aX1 := Max(aTestX, aX1);
	aX2 := Min(aTestX, aX2);

	aClipLine := DLine(
		aCoordSys.ToGlobal(DPoint(aX1, 0)),
		aCoordSys.ToGlobal(DPoint(aX2, 0)));
end;


end.

