Skip to content

Instantly share code, notes, and snippets.

@rkolka
Last active July 4, 2021 20:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save rkolka/edd177f01d24e5f72044461b10bc6e19 to your computer and use it in GitHub Desktop.
Save rkolka/edd177f01d24e5f72044461b10bc6e19 to your computer and use it in GitHub Desktop.
Manifold 9 SQL functions for vector and geometry manipulation
-- $manifold$
-- $include$ [Geom_Functions]
VALUE @p00 GEOM = p2(0,0) ;
VALUE @p01 GEOM = p2(0,1) ;
VALUE @p10 GEOM = p2(1,0) ;
VALUE @p11 GEOM = p2(1,1) ;
VALUE @p000 GEOM = p3(0,0,0) ;
VALUE @p100 GEOM = p3(1,0,0) ;
VALUE @p010 GEOM = p3(0,1,0) ;
VALUE @p001 GEOM = p3(0,0,1) ;
VALUE @p234 GEOM = p3(2,3,4) ;
VALUE @p34 GEOM = p2(3,4) ;
VALUE @ihat FLOAT64X3 = v3(1,0,0) ;
VALUE @jhat FLOAT64X3 = v3(0,1,0) ;
VALUE @khat FLOAT64X3 = v3(0,0,1) ;
-- $manifold$
-- $include$ [Vector_Functions]
-- Short accessors of individual coord values
FUNCTION x(@g GEOM) FLOAT64 AS VectorValue( GeomCoordXY(@g, 0), 0) END ;
FUNCTION y(@g GEOM) FLOAT64 AS VectorValue( GeomCoordXY(@g, 0), 1) END ;
FUNCTION z(@g GEOM) FLOAT64 AS VectorValue( GeomCoordXYZ(@g, 0), 2) END ;
-- Analogous to M8 CentroidX/CentroidX
FUNCTION cx(@g GEOM) FLOAT64 AS VectorValue( GeomCenter(@g, 0), 0) END ;
FUNCTION cy(@g GEOM) FLOAT64 AS VectorValue( GeomCenter(@g, 0), 1) END ;
FUNCTION cxy(@g GEOM) FLOAT64X2 AS VectorMakeX2( VectorValue(GeomCenter(@g, 0), 0), VectorValue(GeomCenter(@g, 0), 1) ) END ;
-- Short version for creating points
FUNCTION p2(@x FLOAT64, @y FLOAT64) GEOM AS GeomMakePoint(VectorMakeX2(@x, @y)) END
FUNCTION p3(@x FLOAT64, @y FLOAT64, @z FLOAT64) GEOM AS GeomMakePoint3(VectorMakeX3(@x, @y, @z)) END
FUNCTION g2(@xy FLOAT64X2) GEOM AS GeomMakePoint(@xy) END
FUNCTION g3(@xyz FLOAT64X3) GEOM AS GeomMakePoint3(@xyz) END
--
FUNCTION StartPoint2(@g GEOM) GEOM AS GeomMakePoint(GeomCoordXY(@g, 0)) END ;
FUNCTION EndPoint2(@g GEOM) GEOM AS GeomMakePoint(GeomCoordXY(@g, GeomCoordCount(@g)-1)) END ;
FUNCTION StartPoint3(@g GEOM) GEOM AS GeomMakePoint3(GeomCoordXYZ(@g, 0)) END ;
FUNCTION EndPoint3(@g GEOM) GEOM AS GeomMakePoint3(GeomCoordXYZ(@g, GeomCoordCount(@g)-1)) END ;
-- From geom to vectors
FUNCTION xy(@g GEOM) FLOAT64X2 AS GeomCoordXY(@g, 0) END ;
FUNCTION xyz(@g GEOM) FLOAT64X3 AS GeomCoordXYZ(@g, 0) END ;
FUNCTION StartPointXY(@g GEOM) FLOAT64X2 AS GeomCoordXY(@g, 0) END ;
FUNCTION EndPointXY(@g GEOM) FLOAT64X2 AS GeomCoordXY(@g, GeomCoordCount(@g)-1) END ;
FUNCTION StartPointXYZ(@g GEOM) FLOAT64X3 AS GeomCoordXYZ(@g, 0) END ;
FUNCTION EndPointXYZ(@g GEOM) FLOAT64X3 AS GeomCoordXYZ(@g, GeomCoordCount(@g)-1) END ;
FUNCTION abXY(@g1 GEOM, @g2 GEOM) FLOAT64X2 AS ab2(StartPointXY(@g1), EndPointXY(@g2)) END ;
FUNCTION abXYZ(@g1 GEOM, @g2 GEOM) FLOAT64X3 AS ab3(StartPointXYZ(@g1), EndPointXYZ(@g2)) END ;
FUNCTION LineToXY(@g GEOM) FLOAT64X2 AS abXY(@g, @g) END ;
FUNCTION LineToXYZ(@g GEOM) FLOAT64X3 AS abXYZ(@g, @g) END ;
FUNCTION GeomIsRing(@g GEOM) BOOLEAN AS GeomCoordXY(@g,0) = GeomCoordXY(@g, GeomCoordCount(@g)-1) END ;
FUNCTION GeomInterpolate(@g GEOM, @q FLOAT64) GEOM AS
(
GeomMakePoint(Interpolate2(StartPointXY(@g), EndPointXY(@g), @q))
)
END
;
FUNCTION GeomMidPoint(@g GEOM) GEOM AS
(
GeomInterpolate(@g, 0.5)
)
END
;
FUNCTION GeomToGrid(@g GEOM, @step FLOAT64) GEOM AS
(
GeomSnapToGrid(@g, VectorMakeX2(@step, @step))
)
END
;
FUNCTION GeomSmoothGrid(@g GEOM, @SmoothValue FLOAT64, @GridStep FLOAT64, @NormalizeValue FLOAT64) GEOM AS
(
GeomNormalize(GeomSnapToGrid(GeomSmooth(@g, @SmoothValue), VectorMakeX2(@GridStep, @GridStep)), @NormalizeValue)
)
END
;
FUNCTION AzimuthSegment(@g GEOM) FLOAT64 AS
(
57.2958 * Atan2((VectorValue(GeomCoordXY(@g, 1),0) - VectorValue(GeomCoordXY(@g, 0),0)), (VectorValue(GeomCoordXY(@g, 1),1) - VectorValue(GeomCoordXY(@g, 0),1)))
)
END
;
FUNCTION GeomScaleAt2(@g GEOM, @scale FLOAT64X2, @at FLOAT64X2) GEOM AS
(
GeomScaleShift(@g, @scale, ab2(scaleComponents2(@at, @scale), @at))
)
END
;
FUNCTION GeomScaleAtCenter(@g GEOM, @scale FLOAT64X2) GEOM AS
(
GeomScaleAt2(@g, @scale, cxy(@g))
)
END
;
FUNCTION GeomProjectZUp(@geom GEOM, @orientation FLOAT64x2) GEOM AS
p2( x2(@orientation)*x(@geom) + y2(@orientation)*y(@geom), z(@geom) )
END
;
FUNCTION GeomToPixelDepth(@dims INT32X3, @origin GEOM, @geom GEOM ) GEOM AS
g3(SphericalToPixelDepth(@dims, SphericalCoordsCentered(xyz(@origin), xyz(@geom))))
END
;
-- $include$ [Geom_Functions]
-- $include$ [Constants]
-- select one expression and use ALT-SHIFT-ENTER to evaluate.
v2(3,3)
v2(3,-3)
v3(4,3,1)
v3(-4,3,2)
v4(4,3,1,1)
v4(-4,3,2,-2)
dot2(v2(3,3), v2(3,-1))
dot3(v3(3,3,1), v3(3,-3,-1))
dot4(v4(4,3,1,1), v4(-4,3,2,-2))
norm3(hat3(cross3(v3(3,3,1), v3(3,-3,-1))))
add3(v3(5,7,8), v3(3,2,1))
scale3(v3(6,4,2), 0.21)
x2(v2(34,35))
xy(p2(34,35))
xy(@p00)
xyz(@p000)
-- for VALUES ... use ALT-ENTER to evaluate
VALUES
( 'dummy', 1),
( 'x2 x', x2(v2(4545.9999999, 8.9999)) ) ,
( 'x2 y', y2(v2(4545.9999999, 8.9999)) ) ,
( 'x3 x', x3(v3(4545.9999999, 8.439, 4334.43434334)) ) ,
( 'x3 y', y3(v3(4545.9999999, 8.439, 4334.43434334)) ) ,
( 'x3 z', z3(v3(4545.9999999, 8.439, 4334.43434334)) ) ,
( 'x4 x', x4(v4(4545.9999999, 8.439, 63.43434334, 45.3)) ) ,
( 'x4 y', y4(v4(4545.9999999, 8.439, 63.43434334, 45.3)) ) ,
( 'x4 z', z4(v4(4545.9999999, 8.439, 63.43434334, 45.3)) ) ,
( 'x4 w', w4(v4(4545.9999999, 8.439, 63.43434334, 45.3)) )
AS
([desc], [Value])
;
VALUES
( V2(4545.9999999, 8.439) ,
V3(4545.9999999, 8.439, 63.43434334) ,
V4(4545.9999999, 8.439, 63.43434334, 45.3)
)
as ([V2], [V3], [V4] )
;
xy(StringWktGeom('POINT(45 23)'))
-- find slope in degrees
Rad2Deg(acos(dot3(
v3(0,0,1),
hat3(cross3(
ab3(v3(7290.41, 116.5, 483.52), v3(7291.87, 115.47, 390.59))
,
ab3(v3(7290.41, 116.5, 483.52), v3(7289.92, 113.62, 389.58))
))
)))
TriangleInclineRatio(v3(0,0,35), v3(4,2,0), v3(4,-2,0))
Rad2Deg(TriangleAspect(v3(0,0,35), v3(4,-2,0), v3(4,2,0)))
cross3( ab3(v3(0,0,35), v3(4,2,0)) , ab3(v3(0,0,35), v3(4,-2,0)) )
TriangleInclineRatio(v3(90.41, 6.5, 483.52), v3(91.87, 5.47, 390.59), v3(90, 3.62, 389.58))
Rad2Deg(TriangleAspect(v3(90.41, 6.5, 483.52), v3(91.87, 5.47, 390.59), v3(90, 3.62, 389.58)))
Rad2Deg(TriangleAspect(v3(1, 1, 1), v3(2, 2, -2), v3(-2, 2, -2)))
Rad2Deg(AngleBetween3(v3(1,2,3), v3(-2,-1,-3)))
norm3(RotateAroundAxis3(v3(1000,-2000,3000), v3(-12040060606,-888,-60346234463), Deg2Rad(-475754)))
norm3(v3(1000,-2000,3000))
Rad2Deg(AngleBetween3(v3(1,1,1), v3(-1,1,1)))
RotateTowards3(v3(1,1,1), v3(-1,1,1), Deg2Rad(90))
Rad2Deg(AngleBetween3(v3(1,1,1), v3(-1,1,1)))
+
Rad2Deg(AngleBetween3(v3(-1,1,1),
RotateTowards3(v3(1,1,1), v3(-1,1,1), Deg2Rad(90)) ))
AffineTransform1(6, 0.5, 2)
AffineTransform1(-6, 0.5, 2)
AffineTransform2(v2(-1,-1), v2(2, -1), v2(2, 1), v2(2, 2))
setz3(v2(5,6),7)
-- $manifold$
-- Shorter aliases for creating Vectors
FUNCTION v2(@x FLOAT64, @y FLOAT64) FLOAT64X2 AS VectorMakeX2(@x, @y) END ;
FUNCTION v3(@x FLOAT64, @y FLOAT64, @z FLOAT64) FLOAT64X3 AS VectorMakeX3(@x, @y, @z) END ;
FUNCTION v4(@x FLOAT64, @y FLOAT64, @z FLOAT64, @w FLOAT64) FLOAT64X4 AS VectorMakeX4(@x, @y, @z, @w) END ;
-- Short accessors of vector components
FUNCTION x2(@v FLOAT64X2) FLOAT64 AS VectorValue(@v, 0) END ;
FUNCTION y2(@v FLOAT64X2) FLOAT64 AS VectorValue(@v, 1) END ;
FUNCTION x3(@v FLOAT64X3) FLOAT64 AS VectorValue(@v, 0) END ;
FUNCTION y3(@v FLOAT64X3) FLOAT64 AS VectorValue(@v, 1) END ;
FUNCTION z3(@v FLOAT64X3) FLOAT64 AS VectorValue(@v, 2) END ;
FUNCTION x4(@v FLOAT64X4) FLOAT64 AS VectorValue(@v, 0) END ;
FUNCTION y4(@v FLOAT64X4) FLOAT64 AS VectorValue(@v, 1) END ;
FUNCTION z4(@v FLOAT64X4) FLOAT64 AS VectorValue(@v, 2) END ;
FUNCTION w4(@v FLOAT64X4) FLOAT64 AS VectorValue(@v, 3) END ;
-- Short setters of vector components
FUNCTION setx2(@v FLOAT64X2, @x FLOAT64) FLOAT64X2 AS v2(@x, y2(@v)) END ;
FUNCTION sety2(@v FLOAT64X2, @y FLOAT64) FLOAT64X2 AS v2(x2(@v), @y) END ;
FUNCTION setx3(@v FLOAT64X3, @x FLOAT64) FLOAT64X3 AS v3(@x, y4(@v), z4(@v) ) END ;
FUNCTION sety3(@v FLOAT64X3, @y FLOAT64) FLOAT64X3 AS v3(x4(@v), @y, z4(@v) ) END ;
FUNCTION setz3(@v FLOAT64X3, @z FLOAT64) FLOAT64X3 AS v3(x4(@v), y4(@v), @z ) END ;
FUNCTION setx4(@v FLOAT64X4, @x FLOAT64) FLOAT64X4 AS v4(@x, y4(@v), z4(@v), w4(@v) ) END ;
FUNCTION sety4(@v FLOAT64X4, @y FLOAT64) FLOAT64X4 AS v4(x4(@v), @y, z4(@v), w4(@v) ) END ;
FUNCTION setz4(@v FLOAT64X4, @z FLOAT64) FLOAT64X4 AS v4(x4(@v), y4(@v), @z, w4(@v) ) END ;
FUNCTION setw4(@v FLOAT64X4, @w FLOAT64) FLOAT64X4 AS v4(x4(@v), y4(@v), z4(@v), @w ) END ;
-- drop 3rd and/or 4th dimension
FUNCTION v2f3(@v FLOAT64X3) FLOAT64X2 AS v2( x3(@v), y3(@v) ) END ;
FUNCTION v2f4(@v FLOAT64X4) FLOAT64X2 AS v2( x4(@v), y4(@v) ) END ;
FUNCTION v3f4(@v FLOAT64X4) FLOAT64X3 AS v3( x4(@v), y4(@v), z4(@v) ) END ;
-- regular vector to homogeneous coordinates
FUNCTION v2thom(@v FLOAT64X2) FLOAT64X3 AS v3( x2(@v), y2(@v), 1 ) END ;
FUNCTION v3thom(@v FLOAT64X3) FLOAT64X4 AS v4( x3(@v), y3(@v), z3(@v), 1) END ;
-- regular vector from homogeneous coordinates
FUNCTION v2fhom(@v FLOAT64X3) FLOAT64X2 AS v2( x3(@v)/z3(@v), y3(@v)/z3(@v) ) END ;
FUNCTION v3fhom(@v FLOAT64X4) FLOAT64X3 AS v3( x4(@v)/w4(@v), y4(@v)/w4(@v), z4(@v)/w4(@v) ) END ;
FUNCTION perp2(@v FLOAT64X2) FLOAT64X2 AS v2( -y2(@v), x2(@v) ) END ;
FUNCTION Atan2xCCW(@v FLOAT64X2) FLOAT64 AS Atan2( y2(@v), x2(@v) ) END ;
FUNCTION Atan2xCW(@v FLOAT64X2) FLOAT64 AS Atan2( x2(@v), y2(@v) ) END ;
--
FUNCTION EqualsEpsilon1(@v FLOAT64, @u FLOAT64, @epsilon FLOAT64) BOOLEAN AS
(
Abs( @v - @u ) <= @epsilon
)
END ;
FUNCTION EqualsEpsilon2(@v FLOAT64X2, @u FLOAT64X2, @epsilon FLOAT64) BOOLEAN AS
(
Abs( x2(@v) - x2(@u) ) <= @epsilon
AND Abs( y2(@v) - y2(@u) ) <= @epsilon
)
END ;
FUNCTION EqualsEpsilon3(@v FLOAT64X3, @u FLOAT64X3, @epsilon FLOAT64) BOOLEAN AS
(
Abs( x3(@v) - x3(@u) ) <= @epsilon
AND Abs( y3(@v) - y3(@u) ) <= @epsilon
AND Abs( z3(@v) - z3(@u) ) <= @epsilon
)
END ;
FUNCTION EqualsEpsilon4(@v FLOAT64X4, @u FLOAT64X4, @epsilon FLOAT64) BOOLEAN AS
(
Abs( x4(@v) - x4(@u) ) <= @epsilon
AND Abs( y4(@v) - y4(@u) ) <= @epsilon
AND Abs( z4(@v) - z4(@u) ) <= @epsilon
AND Abs( w4(@v) - w4(@u) ) <= @epsilon
)
END ;
-- Add vectors / also ShiftComponents
FUNCTION add2(@a FLOAT64X2, @b FLOAT64X2) FLOAT64X2 AS v2(x2(@a) + x2(@b), y2(@a) + y2(@b) ) END ;
FUNCTION add3(@a FLOAT64X3, @b FLOAT64X3) FLOAT64X3 AS v3(x3(@a) + x3(@b), y3(@a) + y3(@b), z3(@a) + z3(@b) ) END ;
FUNCTION add4(@a FLOAT64X4, @b FLOAT64X4) FLOAT64X4 AS v4(x4(@a) + x4(@b), y4(@a) + y4(@b), z4(@a) + z4(@b), w4(@a) + w4(@b) ) END ;
-- Makes Vector ab === Subtracts @a from @b
FUNCTION ab2(@a FLOAT64X2, @b FLOAT64X2) FLOAT64X2 AS v2(x2(@b) - x2(@a), y2(@b) - y2(@a) ) END ;
FUNCTION ab3(@a FLOAT64X3, @b FLOAT64X3) FLOAT64X3 AS v3(x3(@b) - x3(@a), y3(@b) - y3(@a), z3(@b) - z3(@a) ) END ;
FUNCTION ab4(@a FLOAT64X4, @b FLOAT64X4) FLOAT64X4 AS v4(x4(@b) - x4(@a), y4(@b) - y4(@a), z4(@b) - z4(@a), w4(@b) - w4(@a) ) END ;
-- Scale vectors
FUNCTION scale2(@a FLOAT64X2, @b FLOAT64) FLOAT64X2 AS v2(@b*x2(@a), @b*y2(@a)) END ;
FUNCTION scale3(@a FLOAT64X3, @b FLOAT64) FLOAT64X3 AS v3(@b*x3(@a), @b*y3(@a), @b*z3(@a)) END ;
FUNCTION scale4(@a FLOAT64X4, @b FLOAT64) FLOAT64X4 AS v4(@b*x4(@a), @b*y4(@a), @b*z4(@a), @b*w4(@a)) END ;
-- Scale vector components by another vector
FUNCTION scaleComponents2(@a FLOAT64X2, @b FLOAT64X2) FLOAT64X2 AS v2( x2(@a) * x2(@b), y2(@a) * y2(@b) ) END ;
FUNCTION scaleComponents3(@a FLOAT64X3, @b FLOAT64X3) FLOAT64X3 AS v3( x3(@a) * x3(@b), y3(@a) * y3(@b), z3(@a) * z3(@b) ) END ;
FUNCTION scaleComponents4(@a FLOAT64X4, @b FLOAT64X4) FLOAT64X4 AS v4( x4(@a) * x4(@b), y4(@a) * y4(@b), z4(@a) * z4(@b), w4(@a) * w4(@b) ) END ;
-- MostOrthogonalAxis
FUNCTION MostOrthogonalAxis2(@v FLOAT64X2) FLOAT64X2 AS
(
CASE
WHEN Abs(x2(@v)) < Abs(y2(@v)) THEN v2(1,0)
ELSE v2(0,1)
END
)
END ;
FUNCTION MostOrthogonalAxis3(@v FLOAT64X3) FLOAT64X3 AS
(
CASE
WHEN Abs(x3(@v)) < Abs(y3(@v)) AND Abs(x3(@v)) < Abs(z3(@v)) THEN v3(1,0,0)
WHEN Abs(y3(@v)) < Abs(x3(@v)) AND Abs(y3(@v)) < Abs(z3(@v)) THEN v3(0,1,0)
ELSE v3(0,0,1)
END
)
END ;
FUNCTION MostOrthogonalAxis4(@v FLOAT64X4) FLOAT64X4 AS
(
CASE
WHEN Abs(x4(@v)) < Abs(y4(@v)) AND Abs(x4(@v)) < Abs(z4(@v)) AND Abs(x4(@v)) < Abs(w4(@v)) THEN v4(1,0,0,0)
WHEN Abs(y4(@v)) < Abs(x4(@v)) AND Abs(y4(@v)) < Abs(z4(@v)) AND Abs(y4(@v)) < Abs(w4(@v)) THEN v4(0,1,0,0)
WHEN Abs(z4(@v)) < Abs(x4(@v)) AND Abs(z4(@v)) < Abs(y4(@v)) AND Abs(z4(@v)) < Abs(w4(@v)) THEN v4(0,0,1,0)
ELSE v4(0,0,0,1)
END
)
END ;
---
--MaximumComponent
--MinimumComponent
FUNCTION MinimumComponent2(@v FLOAT64X4) FLOAT64 AS
(
CASE
WHEN x4(@v) < y4(@v) THEN x4(@v)
ELSE y4(@v)
END
)
END
;
FUNCTION MinimumComponent3(@v FLOAT64X4) FLOAT64 AS
(
CASE
WHEN x4(@v) < y4(@v) AND x4(@v) < z4(@v) THEN x4(@v)
WHEN y4(@v) < x4(@v) AND y4(@v) < z4(@v) THEN y4(@v)
ELSE z4(@v)
END
)
END
;
FUNCTION MinimumComponent4(@v FLOAT64X4) FLOAT64 AS
(
CASE
WHEN x4(@v) < y4(@v) AND x4(@v) < z4(@v) AND x4(@v) < w4(@v) THEN x4(@v)
WHEN y4(@v) < x4(@v) AND y4(@v) < z4(@v) AND y4(@v) < w4(@v) THEN y4(@v)
WHEN z4(@v) < x4(@v) AND z4(@v) < y4(@v) AND z4(@v) < w4(@v) THEN z4(@v)
ELSE w4(@v)
END
)
END
;
FUNCTION MaximumComponent2(@v FLOAT64X4) FLOAT64 AS
(
CASE
WHEN x4(@v) > y4(@v) THEN x4(@v)
ELSE y4(@v)
END
)
END
;
FUNCTION MaximumComponent3(@v FLOAT64X4) FLOAT64 AS
(
CASE
WHEN x4(@v) > y4(@v) AND x4(@v) > z4(@v) THEN x4(@v)
WHEN y4(@v) > x4(@v) AND y4(@v) > z4(@v) THEN y4(@v)
ELSE z4(@v)
END
)
END
;
FUNCTION MaximumComponent4(@v FLOAT64X4) FLOAT64 AS
(
CASE
WHEN x4(@v) > y4(@v) AND x4(@v) > z4(@v) AND x4(@v) > w4(@v) THEN x4(@v)
WHEN y4(@v) > x4(@v) AND y4(@v) > z4(@v) AND y4(@v) > w4(@v) THEN y4(@v)
WHEN z4(@v) > x4(@v) AND z4(@v) > y4(@v) AND z4(@v) > w4(@v) THEN z4(@v)
ELSE w4(@v)
END
)
END
;
-- dot products VectorDot(<valuexN>, <valuexN>)
-- VectorDot(<valuexN>, <valuexN>)
-- Cross product in 3D. Returns vector perpendicular to both @a and @b.
-- VectorCross(<valuex3>, <valuex3>)
-- norm is the length of vector
FUNCTION norm2(@a FLOAT64X2) FLOAT64 AS
(
Hypot( x2(@a), y2(@a) )
)
END
;
FUNCTION norm3(@a FLOAT64X3) FLOAT64 AS
(
sqrt( pow(x3(@a), 2) + pow(y3(@a), 2) + pow(z3(@a), 2) )
)
END
;
FUNCTION norm4(@a FLOAT64X4) FLOAT64 AS
(
sqrt( pow(x3(@a), 2) + pow(y3(@a), 2) + pow(z3(@a), 2) + pow(z3(@a), 2) )
)
END
;
-- 2d vec as complex number
FUNCTION ComplexFromPolar(@a FLOAT64X2) FLOAT64X2 AS v2( x2(@a)*Cos(y2(@a)), x2(@a)*Sin(y2(@a)) ) END ;
FUNCTION ComplexToPolar(@a FLOAT64X2) FLOAT64X2 AS v2( norm2(@a), Atan2(y2(@a), x2(@a)) ) END ;
-- hat(@a) is a unitvector in the same direction as @a
FUNCTION hat2(@a FLOAT64X2) FLOAT64X2 AS
(
v2(
x2(@a) / norm2(@a),
y2(@a) / norm2(@a)
)
)
END
;
FUNCTION hat3(@a FLOAT64X3) FLOAT64X3 AS
(
v3(
x3(@a) / norm3(@a),
y3(@a) / norm3(@a),
z3(@a) / norm3(@a)
)
)
END
;
FUNCTION hat4(@a FLOAT64X3) FLOAT64X3 AS
(
v4(
x4(@a) / norm4(@a),
y4(@a) / norm4(@a),
z4(@a) / norm4(@a),
w4(@a) / norm4(@a)
)
)
END
;
-- Interpolate or mix a and b
FUNCTION interpolate1(@a FLOAT64, @b FLOAT64, @q FLOAT64) FLOAT64X2 AS
(
@a * (1 - @q) + @b * @q
)
END
;
FUNCTION interpolate2(@a FLOAT64X2, @b FLOAT64X2, @q FLOAT64) FLOAT64X2 AS
(
v2(
x2(@a) * (1 - @q) + x2(@b) * @q ,
y2(@a) * (1 - @q) + y2(@b) * @q
)
)
END
;
FUNCTION interpolate3(@a FLOAT64X3, @b FLOAT64X3, @q FLOAT64) FLOAT64X3 AS
(
v3(
x3(@a) * (1 - @q) + x3(@b) * @q ,
y3(@a) * (1 - @q) + y3(@b) * @q ,
z3(@a) * (1 - @q) + z3(@b) * @q
)
)
END
;
FUNCTION interpolate4(@a FLOAT64X4, @b FLOAT64X4, @q FLOAT64) FLOAT64X4 AS
(
v4(
x4(@a) * (1 - @q) + x4(@b) * @q ,
y4(@a) * (1 - @q) + y4(@b) * @q ,
z4(@a) * (1 - @q) + z4(@b) * @q ,
w4(@a) * (1 - @q) + w4(@b) * @q
)
)
END
;
FUNCTION AngleBetween2(@a FLOAT64X2, @b FLOAT64X2) FLOAT64 AS
(
Acos( VectorDot( hat2(@a), hat2(@b) ) )
)
END
;
FUNCTION AngleBetween3(@a FLOAT64X3, @b FLOAT64X3) FLOAT64 AS
(
Acos( VectorDot( hat3(@a), hat3(@b) ) )
)
END
;
-- use Rodrigues' formula.
FUNCTION RotateAroundAxis3(@v FLOAT64X3, @axis FLOAT64X3, @theta FLOAT64) FLOAT64X3 AS
(
add3(add3(
scale3(hat3(@axis), (1 - Cos(@theta)) * VectorDot(@v, hat3(@axis)))
,
scale3(@v, Cos(@theta)))
,
scale3(VectorCross(@v, hat3(@axis)), Sin(@theta)))
)
END
;
FUNCTION RotateTowards3(@v FLOAT64X3, @b FLOAT64X3, @theta FLOAT64) FLOAT64X3 AS
(
RotateAroundAxis3(@v, VectorCross(@b, @v), @theta)
)
END
;
-- @igt i-hat goes to
-- @jgt j-hat goes to
-- @kgt k-hat goes to
-- @lgt w-hat goes to
-- @Ogt Origin goes to
FUNCTION LinearTransform1(@v FLOAT64, @igt FLOAT64) FLOAT64 AS
(
@igt * @v
)
END
;
FUNCTION LinearTransform2(@v FLOAT64X2, @igt FLOAT64X2, @jgt FLOAT64X2) FLOAT64X2 AS
(
v2(
x2(@igt) * x2(@v) + x2(@jgt) * y2(@v),
y2(@igt) * x2(@v) + y2(@jgt) * y2(@v)
)
)
END
;
FUNCTION LinearTransform3(@v FLOAT64X3, @igt FLOAT64X3, @jgt FLOAT64X3, @kgt FLOAT64X3) FLOAT64X3 AS
(
v3(
x3(@igt) * x3(@v) + x3(@jgt) * y3(@v) + x3(@kgt) * z3(@v) ,
y3(@igt) * x3(@v) + y3(@jgt) * y3(@v) + y3(@kgt) * z3(@v) ,
z3(@igt) * x3(@v) + z3(@jgt) * y3(@v) + z3(@kgt) * z3(@v)
)
)
END
;
FUNCTION LinearTransform4(@v FLOAT64X4, @igt FLOAT64X4, @jgt FLOAT64X4, @kgt FLOAT64X4, @lgt FLOAT64X4) FLOAT64X4 AS
(
v4(
x4(@igt) * x4(@v) + x4(@jgt) * y4(@v) + x4(@kgt) * z4(@v) + x4(@lgt) * w4(@v) ,
y4(@igt) * x4(@v) + y4(@jgt) * y4(@v) + y4(@kgt) * z4(@v) + y4(@lgt) * w4(@v) ,
z4(@igt) * x4(@v) + z4(@jgt) * y4(@v) + z4(@kgt) * z4(@v) + z4(@lgt) * w4(@v) ,
w4(@igt) * x4(@v) + w4(@jgt) * y4(@v) + w4(@kgt) * z4(@v) + w4(@lgt) * w4(@v)
)
)
END
;
FUNCTION AffineTransform1(@v FLOAT64, @igt FLOAT64, @Ogt FLOAT64) FLOAT64 AS
(
x2(LinearTransform2(
v2(@v, 1), v2(@igt, 0), v2(@Ogt, 1)
))
)
END
;
FUNCTION AffineTransform2(@v FLOAT64X2, @igt FLOAT64X2, @jgt FLOAT64X2, @Ogt FLOAT64X2) FLOAT64X2 AS
(
v2f3(LinearTransform3(
setz3(@v, 1), setz3(@igt, 0), setz3(@jgt, 0), setz3(@Ogt, 1)
))
)
END
;
FUNCTION AffineTransform3(@v FLOAT64X3, @igt FLOAT64X3, @jgt FLOAT64X3, @kgt FLOAT64X3, @Ogt FLOAT64X3) FLOAT64X3 AS
(
v3f4(LinearTransform4(
setw4(@v, 1), setw4(@igt, 0), setw4(@jgt, 0), setw4(@kgt, 0), setw4(@Ogt, 1)
))
)
END
;
FUNCTION PerspectiveDivision1(@v FLOAT64X2) FLOAT64 AS
(
x2(@v) / y2(@v)
)
END
;
FUNCTION PerspectiveDivision2(@v FLOAT64X3) FLOAT64X2 AS
(
v2f3(scale3(@v, 1/z3(@v)))
)
END
;
FUNCTION PerspectiveDivision3(@v FLOAT64X4) FLOAT64X3 AS
(
v3f4(scale4(@v, 1/w4(@v)))
)
END
;
FUNCTION PerspectiveTransform1(@v FLOAT64, @igt FLOAT64, @Ogt FLOAT64, @p FLOAT64) FLOAT64 AS
(
PerspectiveDivision1(LinearTransform2(
v2(@v, 1), v2(@igt, @p), v2(@Ogt, 1)
))
)
END
;
FUNCTION PerspectiveAffineTransform1(@v FLOAT64, @igt FLOAT64, @Ogt FLOAT64, @p FLOAT64) FLOAT64 AS
(
PerspectiveDivision1(LinearTransform2(
v2(@v, 1), v2(@igt, @p), v2(@Ogt, 1)
))
)
END
;
FUNCTION PerspectiveTransform2(@v FLOAT64X2, @igt FLOAT64X2, @jgt FLOAT64X2, @Ogt FLOAT64X2, @p FLOAT64X2) FLOAT64X2 AS
(
PerspectiveDivision2(LinearTransform3(
setz3(@v, 1), setz3(@igt, x2(@p)), setz3(@jgt, y2(@p)), setz3(@Ogt, 1)
))
)
END
;
FUNCTION PerspectiveTransform3(@v FLOAT64X3, @igt FLOAT64X3, @jgt FLOAT64X3, @kgt FLOAT64X3, @Ogt FLOAT64X3, @p FLOAT64X3) FLOAT64X3 AS
(
PerspectiveDivision3(LinearTransform4(
setw4(@v, 1), setw4(@igt, x3(@p)), setw4(@jgt, y3(@p)), setw4(@kgt, z3(@p)), setw4(@Ogt, 1)
))
)
END
;
FUNCTION Deg2Rad(@d FLOAT64) FLOAT64 AS @d/180*PI END ;
FUNCTION Rad2Deg(@r FLOAT64) FLOAT64 AS @r/PI*180 END ;
FUNCTION Slope3(@v FLOAT64X3) FLOAT64 AS
(
z3(@v) / Hypot(y3(@v), x3(@v))
)
END
;
--- inclination from the zenith
FUNCTION PolarAngle3(@v FLOAT64X3) FLOAT64 AS
(
Atan2( Hypot(y3(@v), x3(@v)), z3(@v) )
)
END
;
-- elevation from the xy plane
FUNCTION ElevationAngle3(@v FLOAT64X3) FLOAT64 AS
(
Atan2( z3(@v), Hypot(y3(@v), x3(@v)) )
)
END
;
-- counter clock-wise (like in math)
FUNCTION CCWAzimuthAngle3(@v FLOAT64X3) FLOAT64 AS
(
Atan2( y3(@v), x3(@v) )
)
END
;
-- clock-wise (like in geo)
FUNCTION CWAzimuthAngle3(@v FLOAT64X3) FLOAT64 AS
(
Atan2( x3(@v), y3(@v) )
)
END
;
-- clock-wise (like in geo) with offset orientation -
FUNCTION CWOffsetAzimuthAngle3(@v FLOAT64X3, @offset FLOAT64X2) FLOAT64 AS
(
Atan2xCW( LinearTransform2( v2f3(@v), @offset, perp2(@offset) ) )
)
END
;
FUNCTION TriangleInclineRatio(@a FLOAT64X3, @b FLOAT64X3, @c FLOAT64X3) FLOAT64 AS
(
Abs( 1 / Slope3(VectorCross( ab3(@a, @b), ab3(@a, @c) )))
)
END
;
FUNCTION TriangleAspect(@a FLOAT64X3, @b FLOAT64X3, @c FLOAT64X3) FLOAT64 AS
(
CCWAzimuthAngle3( VectorCross(ab3(@a, @b), ab3(@a, @c)) )
)
END
;
FUNCTION SphericalCoords(@v FLOAT64X3) FLOAT64x3 AS
v3( CWAzimuthAngle3(@v), ElevationAngle3(@v), norm3(@v) )
END
;
FUNCTION SphericalCoordsCentered(@origin FLOAT64X3, @v FLOAT64X3) FLOAT64x3 AS
SphericalCoords( ab3(@origin, @v) )
END
;
FUNCTION SphericalCoordsOriented(@orientation FLOAT64X2, @v FLOAT64X3) FLOAT64x3 AS
v3( CWOffsetAzimuthAngle3(@v, @orientation), ElevationAngle3(@v), norm3(@v) )
END
;
FUNCTION SphericalCoordsCenteredOriented(@origin FLOAT64X3, @orientation FLOAT64X2, @v FLOAT64X3) FLOAT64x3 AS
SphericalCoordsOriented(@orientation, ab3(@origin, @v) )
END
;
FUNCTION SphericalToPixelDepth(@dims INT32X3, @v FLOAT64X3) INT32X3 AS
add3(
v3(x3(@dims)/2, y3(@dims)/2, 0),
v3(
x3(@dims)*(1/(2*Pi))*x3(@v),
y3(@dims)*(1/Pi)*y3(@v),
z3(@dims)*z3(@v)
)
)
END
;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment