
'GEOPLANAR.INC


type point
  float x,y
end type

type line
  point p1,p2
end type

type boundingrect
  point p1,p2
end type

type pointex
  float x,y
  int   n     'point count
  word  pp[32]'list of connected pointex
end type
'see moveableobject

type linex
  int     p1  'first pointex index
  int     p2  'second point index
  float   le  'length
  int     dl  'deletion flag
end type



'2D MATH GRAPHICS SUPPORT
'========================

sub drawline(sys q, n=1)
========================
float p at (q)
'float*p=q
sys i
glBegin GL_LINES
for i=1 to n
  glVertex2fv p
  @p+=8
  glVertex2fv p
  @p+=8
next 
glEnd
end sub

sub drawpoint(sys q, n=1)
=========================
float p at (q)
sys i
glBegin GL_POINTS
for i=1 to n
  glVertex2fv p : @p+=8
next
glEnd
end sub


=================================

uses generics


'function samepos(point *p,*q) as sys
'====================================
'if (p.x = q.x)and(p.y = q.y)
'  return 1
'endif
'end function

macro samepos int(r,p,q,  pa,qa)
================================
sys pa,qa
pa=@p
qa=@q
mov rcx,pa
mov rdx,qa
mov esi,1
mov edi,0
mov eax,[rcx]
cmp eax,[rdx]
cmovnz esi,edi 
mov eax,[rcx+4]
cmp eax,[rdx+4]
cmovnz esi,edi 
mov r,esi
end macro


macro sqdist float(r,p,q, pa,qa)
================================
'for float,point,pointex,moveableObject
'assumes float variables x,y
sys pa,qa
pa=@p
qa=@q
mov rcx,pa
mov rdx,qa
fld dword[rcx]
fsub dword[rdx]
fmul st0
fld dword[rcx+4]
fsub dword[rdx+4]
fmul st0
faddp st1
fstp dword r
end macro

macro dist (a,b)
================
hypot( b.x-a.x, b.y-a.y)
end macro

macro inrange int(r,m1,m2,a,  k)
================================
'uses simple float terms
'check m1 & m2 for min-max order

fld dword a
fld dword m1
fld dword m2  'onstack m2 m1 a
fcomi         'm2-m1
jae fwd k     'swap if (jb) m2<m1
fxch st1      'm1 m2 a
k:
mov eax,-1    'true
mov ecx,0     'to falsify
              'max min a
fxch st2      'a min max
fxch st1      'min a max
fcomip        'min-a
cmova eax,ecx 'falsify if min>a
fxch st1      'max a
fcomip        'max-a
cmovb eax,ecx 'falsify if max<a
fstp st0      'dump a
mov r,eax
end macro

'function dist(point *a,*b) as float
'===================================
'return hypot( b.x-a.x, b.y-a.y)
'end function



function inrange(float m1,m2,p) as int
======================================
if p<=m2 
  if p>=m1
    return -1
  endif
endif
if p<=m1 
  if p>=m2
    return -1
  endif
endif
end function

function insiderange(float m1,m2,p) as int
==========================================
if p<m2 
  if p>m1
    return 1
  endif
endif
if p<m1 
  if p>m2
    return 1
  endif
endif
end function

function inbetween(point *p,*q,*r) as int
=========================================
if insiderange(p.x,q.x,r.x)
  if insiderange(p.y,q.y,r.y)
    return 1
  endif
endif
end function

function interpoint(point *a,*b,*p, float f)
============================================
p.x=a.x+(b.x-a.x)*f
p.y=a.y+(b.y-a.y)*f
end function

function midpoint(point *a,*b,*p)
=================================
static float half=0.5
p.x=a.x+(b.x-a.x)*half
p.y=a.y+(b.y-a.y)*half
end function

function rightangleBase(point *pa,*pb,*pc,*pd) as int
=====================================================
'anticlockwise triangle points:
'left  pa
'right pb
'apex  pc
'TriangleHt will be  between pc pd, forming 2 triangles
float a,b,c,d
a=sqdist pc,pa 'LEFT SIDE
b=sqdist pc,pb 'RIGHT SIDE
c=sqdist pa,pb 'BASE
if c<>0 then
  float rc=sqr(c)
  d=(a-b+c)/(2*rc)
  interpoint pa,pb,pd, d/rc
  return 1
end if
end function

function TriangleHt(point*pa,*pb,*pc,float*h) as int
====================================================
point pd
rightanglebase pa,pb,pc,pd
float ac,ad
ac=sqdist pa,pc
ad=sqdist pa,pd
if ac>ad
  h=sqr(ac-ad) 'height perpendicular
  return 1
endif
end function

function triangleArea(point a,b,c,float*ar) as int
==================================================
float ab,h
ab=dist a,b
TriangleHt a,b,c,h
ar=0.5*ab*h
return 1
end function

function approx(float a,b,pr) as int
====================================
'positive values
if a+pr>b
  if b+pr>a
    return 1
  endif
endif
'negative values
if a+pr<b
  if b+pr<a
    return 1
  endif
endif
end function

function insidetriangle(point*a,*b,*c,*p,float ap) as int
=========================================================
'using area matching
float ar1,ar2,ar3,art
trianglearea a,b,c,art
trianglearea a,b,p,ar1
trianglearea b,c,p,ar2
trianglearea c,a,p,ar3
'prox 1.0001 say
return approx art,ar1+ar2+ar3, ap
end function

function calcgradient(line *l, float *g,*i)
===========================================
float x = l.p2.x - l.p1.x
float y = l.p2.y - l.p1.y
if x=0 then
  g=0 : i = 1 'SET VERTICAL LINE FLAG
else
  g = y/x : i=0
end if
end function

function rot90(point *p,*q,*r)
==============================
'rotate q around p by 90 degrees anticlockwise
float ddx,ddy
ddx = q.x-p.x
ddy = q.y-p.y
r.x = p.x-ddy
r.y = p.y+ddx
end function

function rotm90(point *p,*q,*r)
===============================
'rotate q around p by 90 degrees clockwise
float ddx,ddy
ddx = q.x-p.x
ddy = q.y-p.y
r.x = p.x+ddy
r.y = p.y-ddx
end function

function linerot90(point *p,line *q,*r)
=======================================
rot90 p,q.p1,r.p1
rot90 p,q.p2,r.p2
end function

function linerotm90(point *p,line *q,*r)
========================================
rotm90 p,q.p1,r.p1
rotm90 p,q.p2,r.p2
end function




function topolar(point *p,*q,*r) as sys
=======================================
float dx,dy
dx=q.x-p.x
dy=q.y-p.y
r.x=hypot(dx,dy) 'as polar dist
r.y=atan(dy,dx)  'as polar angle
return 1
end function

function tocart(point *p,*q,*r) as sys
======================================
float dx,dy
dx=r.x*cos(r.y) 'r.x polar dist
dy=r.x*sin(r.y) 'r.y polar angle
q.x+=dx
q.y+=dy
return 1
end function


function exclusive(line *l1,*l2) as sys
=======================================
boundingrect b1,b2
copy @b1,@l1,sizeof boundingrect
copy @b2,@l2,sizeof boundingrect
'
'B1 X
if b1.p1.x>b1.p2.x
  swap b1.p1.x,b1.p2.x
endif
'B1 Y
if b1.p1.y>b1.p2.y
  swap b1.p1.y,b1.p2.y
endif
'B2 X
if b2.p1.x>b2.p2.x
  swap b2.p1.x,b2.p2.x
endif
'B2 Y
if b2.p1.y>b2.p2.y
  swap b2.p1.y,b2.p2.y
endif
'
if b1.p2.x<b2.p1.x then return 1
if b1.p1.x>b2.p2.x then return 1
if b1.p2.y<b2.p1.y then return 1
if b1.p1.y>b2.p2.y then return 1
'
end function

function boundstatus(line *l1,*l2, point *p, int mo=0) as int
=============================================================
int   a=0
float xp=p.x , yp=p.y
'
'COINCIDENT POINT TESTING
'
if samepos(p,l1.p1) then a or=1
if samepos(p,l1.p2) then a or=2
if samepos(p,l2.p1) then a or=4
if samepos(p,l2.p2) then a or=8
if a then return a
if mo=1
  if not insiderange(l1.p1.x,l1.p2.x,xp) then return 0
  if not insiderange(l1.p1.y,l1.p2.y,yp) then return 0
  if not insiderange(l2.p1.x,l2.p2.x,xp) then return 0
  if not insiderange(l2.p1.y,l2.p2.y,yp) then return 0
else
  if not inrange(l1.p1.x,l1.p2.x,xp) then return 0
  if not inrange(l1.p1.y,l1.p2.y,yp) then return 0
  if not inrange(l2.p1.x,l2.p2.x,xp) then return 0
  if not inrange(l2.p1.y,l2.p2.y,yp) then return 0
endif
return 16
end function

function intersect(line *l1,*l2, point *p) as int
=================================================
' 0 no intersection
' 1 same pos as l1.p1
' 2 same pos as l1.p2
' 4 same pos as l2.p1
' 8 same pos as l2.p2
'16 intersection: lines cross between points
'
float g1,g2,i1,i2
calcgradient l1,g1,i1
calcgradient l2,g2,i2
if g1=g2
  if i1=i2
    return 0 'PARALLEL LINES
  endif
endif
float xd,yd,y2v,xp,yp
xd=l2.p1.x - l1.p1.x
if i2 then
  'vertical line l2
  xp=l2.p1.x
  yp=l1.p1.y+g1*xd
elseif i1 then
  'vertical line l1
  xp=l1.p1.x
  yp=l2.p1.y-g2*xd
elseif g2=0
  'horizontal line l2
  yd  = l2.p1.y - l1.p1.y
  xp  = l1.p1.x+yd/g1
  yp  = l2.p1.y
elseif g1=0
  'horizontal line l1
  yd  = l1.p1.y - l2.p1.y
  xp  = l2.p1.x+yd/g2
  yp  = l1.p1.y
else
  y2v = l2.p1.y - xd*g2
  yd  = y2v - l1.p1.y
  xp  = yd/(g1-g2)
  yp  = l1.p1.y + xp*g1
  xp += l1.p1.x
end if
p = {xp,yp}
return 1
end function

function intersected(line *l1,*l2,point *p, int mo=0) as int
============================================================
if exclusive(l1,l2) then
  return 0
elseif intersect(l1,l2,p) then
  return boundstatus l1,l2,p, mo
end if
end function

function intersectedp(sys pl1,pl2,pp) as int
============================================
line l1 at pl1
line l2 at pl2
point p at pp
if exclusive(l1,l2) then
  return 0
elseif intersect(l1,l2,p) then
  return boundstatus l1,l2,p
end if
end function

function centrefrom3points(point *p,*q,*r,*c) as int
====================================================
point mpq,mqr
point npq,nqr
sys   a
midpoint p,q,mpq
midpoint q,r,mqr
rotm90 mpq,q,npq 
rotm90 mqr,r,nqr
line l1,l2
l1.p1.x = mpq.x : l1.p1.y=mpq.y
l1.p2.x = npq.x : l1.p2.y=npq.y
l2.p1.x = mqr.x : l2.p1.y=mqr.y
l2.p2.x = nqr.x : l2.p2.y=nqr.y
a=intersect l1,l2,c
return a
end function

function arcfrom3points(point *p,*q,*r,*c,*s,*e) as int
=======================================================
sys   a
a=centrefrom3points p,q,r,c
if a then
  topolar c,p,s
  topolar c,r,e
end if
return @r
end function

function CreatePolyPoints( point*p,float x=0,y=0,int n=4, float sx=1.0,sy=1.0,ia=0.0) as int
============================================================================================
int i
int e=n+1
float a
for i=1 to e 'closed
  if i=e
    a=0
  endif
  p.x=x+(sx*cos(rad(ia+a)))
  p.y=y+(sy*sin(rad(ia+a)))
  a=i*360/n
  stride p
next
end function

function CreatePolyLines( line*p,float x=0,y=0,int n=4, float sx=1.0,sy=1.0) as int
===================================================================================
int i
float a
for i=1 to n
  p.p1.x=x+(sx*cos(rad(a)))
  p.p1.y=y+(sy*sin(rad(a)))
  a=i*360/n
  p.p2.x=x+(sx*cos(rad(a)))
  p.p2.y=y+(sy*sin(rad(a)))
  stride p
next
end function

function PerturbPoints( sys pp, int n=4,float x=0.01,y=0.01 ) as int
====================================================================
point p at (pp)
int i
for i=1 to n
  p.x+=x*rnd
  p.y+=y*rnd
  stride p
next
end function
