1// geomx adds some geometry functions needed by rasterx
2// Copyright 2017 by the rasterx Authors. All rights reserved.
3// Created: 2/12/2017 by S.R.Wiley
4
5package rasterx
6
7import (
8 "fmt"
9 "math"
10
11 "golang.org/x/image/math/fixed"
12)
13
14// Invert returns the point inverted around the origin
15func Invert(v fixed.Point26_6) fixed.Point26_6 {
16 return fixed.Point26_6{X: -v.X, Y: -v.Y}
17}
18
19// turnStarboard90 returns the vector 90 degrees starboard (right in direction heading)
20func turnStarboard90(v fixed.Point26_6) fixed.Point26_6 {
21 return fixed.Point26_6{X: -v.Y, Y: v.X}
22}
23
24// turnPort90 returns the vector 90 degrees port (left in direction heading)
25func turnPort90(v fixed.Point26_6) fixed.Point26_6 {
26 return fixed.Point26_6{X: v.Y, Y: -v.X}
27}
28
29// DotProd returns the inner product of p and q
30func DotProd(p fixed.Point26_6, q fixed.Point26_6) fixed.Int52_12 {
31 return fixed.Int52_12(int64(p.X)*int64(q.X) + int64(p.Y)*int64(q.Y))
32}
33
34// Length is the distance from the origin of the point
35func Length(v fixed.Point26_6) fixed.Int26_6 {
36 vx, vy := float64(v.X), float64(v.Y)
37 return fixed.Int26_6(math.Sqrt(vx*vx + vy*vy))
38}
39
40//PathCommand is the type for the path command token
41type PathCommand fixed.Int26_6
42
43// Human readable path constants
44const (
45 PathMoveTo PathCommand = iota
46 PathLineTo
47 PathQuadTo
48 PathCubicTo
49 PathClose
50)
51
52// A Path starts with a PathCommand value followed by zero to three fixed
53// int points.
54type Path []fixed.Int26_6
55
56// ToSVGPath returns a string representation of the path
57func (p Path) ToSVGPath() string {
58 s := ""
59 for i := 0; i < len(p); {
60 if i != 0 {
61 s += " "
62 }
63 switch PathCommand(p[i]) {
64 case PathMoveTo:
65 s += fmt.Sprintf("M%4.3f,%4.3f", float32(p[i+1])/64, float32(p[i+2])/64)
66 i += 3
67 case PathLineTo:
68 s += fmt.Sprintf("L%4.3f,%4.3f", float32(p[i+1])/64, float32(p[i+2])/64)
69 i += 3
70 case PathQuadTo:
71 s += fmt.Sprintf("Q%4.3f,%4.3f,%4.3f,%4.3f", float32(p[i+1])/64, float32(p[i+2])/64,
72 float32(p[i+3])/64, float32(p[i+4])/64)
73 i += 5
74 case PathCubicTo:
75 s += "C" + fmt.Sprintf("C%4.3f,%4.3f,%4.3f,%4.3f,%4.3f,%4.3f", float32(p[i+1])/64, float32(p[i+2])/64,
76 float32(p[i+3])/64, float32(p[i+4])/64, float32(p[i+5])/64, float32(p[i+6])/64)
77 i += 7
78 case PathClose:
79 s += "Z"
80 i++
81 default:
82 panic("freetype/rasterx: bad pather")
83 }
84 }
85 return s
86}
87
88// String returns a readable representation of a Path.
89func (p Path) String() string {
90 return p.ToSVGPath()
91}
92
93// Clear zeros the path slice
94func (p *Path) Clear() {
95 *p = (*p)[:0]
96}
97
98// Start starts a new curve at the given point.
99func (p *Path) Start(a fixed.Point26_6) {
100 *p = append(*p, fixed.Int26_6(PathMoveTo), a.X, a.Y)
101}
102
103// Line adds a linear segment to the current curve.
104func (p *Path) Line(b fixed.Point26_6) {
105 *p = append(*p, fixed.Int26_6(PathLineTo), b.X, b.Y)
106}
107
108// QuadBezier adds a quadratic segment to the current curve.
109func (p *Path) QuadBezier(b, c fixed.Point26_6) {
110 *p = append(*p, fixed.Int26_6(PathQuadTo), b.X, b.Y, c.X, c.Y)
111}
112
113// CubeBezier adds a cubic segment to the current curve.
114func (p *Path) CubeBezier(b, c, d fixed.Point26_6) {
115 *p = append(*p, fixed.Int26_6(PathCubicTo), b.X, b.Y, c.X, c.Y, d.X, d.Y)
116}
117
118// Stop joins the ends of the path
119func (p *Path) Stop(closeLoop bool) {
120 if closeLoop {
121 *p = append(*p, fixed.Int26_6(PathClose))
122 }
123}
124
125// AddTo adds the Path p to q.
126func (p Path) AddTo(q Adder) {
127 for i := 0; i < len(p); {
128 switch PathCommand(p[i]) {
129 case PathMoveTo:
130 q.Stop(false) // Fixes issues #1 by described by Djadala; implicit close if currently in path.
131 q.Start(fixed.Point26_6{X: p[i+1], Y: p[i+2]})
132 i += 3
133 case PathLineTo:
134 q.Line(fixed.Point26_6{X: p[i+1], Y: p[i+2]})
135 i += 3
136 case PathQuadTo:
137 q.QuadBezier(fixed.Point26_6{X: p[i+1], Y: p[i+2]}, fixed.Point26_6{X: p[i+3], Y: p[i+4]})
138 i += 5
139 case PathCubicTo:
140 q.CubeBezier(fixed.Point26_6{X: p[i+1], Y: p[i+2]},
141 fixed.Point26_6{X: p[i+3], Y: p[i+4]}, fixed.Point26_6{X: p[i+5], Y: p[i+6]})
142 i += 7
143 case PathClose:
144 q.Stop(true)
145 i++
146 default:
147 panic("AddTo: bad path")
148 }
149 }
150 q.Stop(false)
151}
152
153// ToLength scales the point to the length indicated by ln
154func ToLength(p fixed.Point26_6, ln fixed.Int26_6) (q fixed.Point26_6) {
155 if ln == 0 || (p.X == 0 && p.Y == 0) {
156 return
157 }
158
159 pX, pY := float64(p.X), float64(p.Y)
160 lnF := float64(ln)
161 pLen := math.Sqrt(pX*pX + pY*pY)
162
163 qX, qY := pX*lnF/pLen, pY*lnF/pLen
164 q.X, q.Y = fixed.Int26_6(qX), fixed.Int26_6(qY)
165 return
166}
167
168// ClosestPortside returns the closest of p1 or p2 on the port side of the
169// line from the bow to the stern. (port means left side of the direction you are heading)
170// isIntersecting is just convienice to reduce code, and if false returns false, because p1 and p2 are not valid
171func ClosestPortside(bow, stern, p1, p2 fixed.Point26_6, isIntersecting bool) (xt fixed.Point26_6, intersects bool) {
172 if isIntersecting == false {
173 return
174 }
175 dir := bow.Sub(stern)
176 dp1 := p1.Sub(stern)
177 dp2 := p2.Sub(stern)
178 cp1 := dir.X*dp1.Y - dp1.X*dir.Y
179 cp2 := dir.X*dp2.Y - dp2.X*dir.Y
180 switch {
181 case cp1 < 0 && cp2 < 0:
182 return
183 case cp1 < 0 && cp2 >= 0:
184 return p2, true
185 case cp1 >= 0 && cp2 < 0:
186 return p1, true
187 default: // both points on port side
188 dirdot := DotProd(dir, dir)
189 // calculate vector rejections of dp1 and dp2 onto dir
190 h1 := dp1.Sub(dir.Mul(fixed.Int26_6((DotProd(dp1, dir) << 6) / dirdot)))
191 h2 := dp2.Sub(dir.Mul(fixed.Int26_6((DotProd(dp2, dir) << 6) / dirdot)))
192 // return point with smallest vector rejection; i.e. closest to dir line
193 if (h1.X*h1.X + h1.Y*h1.Y) > (h2.X*h2.X + h2.Y*h2.Y) {
194 return p2, true
195 }
196 return p1, true
197 }
198}
199
200// RadCurvature returns the curvature of a Bezier curve end point,
201// given an end point, the two adjacent control points and the degree.
202// The sign of the value indicates if the center of the osculating circle
203// is left or right (port or starboard) of the curve in the forward direction.
204func RadCurvature(p0, p1, p2 fixed.Point26_6, dm fixed.Int52_12) fixed.Int26_6 {
205 a, b := p2.Sub(p1), p1.Sub(p0)
206 abdot, bbdot := DotProd(a, b), DotProd(b, b)
207 h := a.Sub(b.Mul(fixed.Int26_6((abdot << 6) / bbdot))) // h is the vector rejection of a onto b
208 if h.X == 0 && h.Y == 0 { // points are co-linear
209 return 0
210 }
211 radCurve := fixed.Int26_6((fixed.Int52_12(a.X*a.X+a.Y*a.Y) * dm / fixed.Int52_12(Length(h)<<6)) >> 6)
212 if a.X*b.Y > b.X*a.Y { // xprod sign
213 return radCurve
214 }
215 return -radCurve
216}
217
218// CircleCircleIntersection calculates the points of intersection of
219// two circles or returns with intersects == false if no such points exist.
220func CircleCircleIntersection(ct, cl fixed.Point26_6, rt, rl fixed.Int26_6) (xt1, xt2 fixed.Point26_6, intersects bool) {
221 dc := cl.Sub(ct)
222 d := Length(dc)
223
224 // Check for solvability.
225 if d > (rt + rl) {
226 return // No solution. Circles do not intersect.
227 }
228 // check if d < abs(rt-rl)
229 if da := rt - rl; (da > 0 && d < da) || (da < 0 && d < -da) {
230 return // No solution. One circle is contained by the other.
231 }
232
233 rlf, rtf, df := float64(rl), float64(rt), float64(d)
234 af := (rtf*rtf - rlf*rlf + df*df) / df / 2.0
235 hfd := math.Sqrt(rtf*rtf-af*af) / df
236 afd := af / df
237
238 rOffx, rOffy := float64(-dc.Y)*hfd, float64(dc.X)*hfd
239 p2x := float64(ct.X) + float64(dc.X)*afd
240 p2y := float64(ct.Y) + float64(dc.Y)*afd
241 xt1x, xt1y := p2x+rOffx, p2y+rOffy
242 xt2x, xt2y := p2x-rOffx, p2y-rOffy
243 return fixed.Point26_6{X: fixed.Int26_6(xt1x), Y: fixed.Int26_6(xt1y)},
244 fixed.Point26_6{X: fixed.Int26_6(xt2x), Y: fixed.Int26_6(xt2y)}, true
245}
246
247// CalcIntersect calculates the points of intersection of two fixed point lines
248// and panics if the determinate is zero. You have been warned.
249func CalcIntersect(a1, a2, b1, b2 fixed.Point26_6) (x fixed.Point26_6) {
250 da, db, ds := a2.Sub(a1), b2.Sub(b1), a1.Sub(b1)
251 det := float32(da.X*db.Y - db.X*da.Y) // Determinate
252 t := float32(ds.Y*db.X-ds.X*db.Y) / det
253 x = a1.Add(fixed.Point26_6{X: fixed.Int26_6(float32(da.X) * t), Y: fixed.Int26_6(float32(da.Y) * t)})
254 return
255}
256
257// RayCircleIntersection calculates the points of intersection of
258// a ray starting at s2 passing through s1 and a circle in fixed point.
259// Returns intersects == false if no solution is possible. If two
260// solutions are possible, the point closest to s2 is returned
261func RayCircleIntersection(s1, s2, c fixed.Point26_6, r fixed.Int26_6) (x fixed.Point26_6, intersects bool) {
262 fx, fy, intersects := RayCircleIntersectionF(float64(s1.X), float64(s1.Y),
263 float64(s2.X), float64(s2.Y), float64(c.X), float64(c.Y), float64(r))
264 return fixed.Point26_6{X: fixed.Int26_6(fx),
265 Y: fixed.Int26_6(fy)}, intersects
266
267}
268
269// RayCircleIntersectionF calculates in floating point the points of intersection of
270// a ray starting at s2 passing through s1 and a circle in fixed point.
271// Returns intersects == false if no solution is possible. If two
272// solutions are possible, the point closest to s2 is returned
273func RayCircleIntersectionF(s1X, s1Y, s2X, s2Y, cX, cY, r float64) (x, y float64, intersects bool) {
274 n := s2X - cX // Calculating using 64* rather than divide
275 m := s2Y - cY
276
277 e := s2X - s1X
278 d := s2Y - s1Y
279
280 // Quadratic normal form coefficients
281 A, B, C := e*e+d*d, -2*(e*n+m*d), n*n+m*m-r*r
282
283 D := B*B - 4*A*C
284
285 if D <= 0 {
286 return // No intersection or is tangent
287 }
288
289 D = math.Sqrt(D)
290 t1, t2 := (-B+D)/(2*A), (-B-D)/(2*A)
291 p1OnSide := t1 > 0
292 p2OnSide := t2 > 0
293
294 switch {
295 case p1OnSide && p2OnSide:
296 if t2 < t1 { // both on ray, use closest to s2
297 t1 = t2
298 }
299 case p2OnSide: // Only p2 on ray
300 t1 = t2
301 case p1OnSide: // only p1 on ray
302 default: // Neither solution is on the ray
303 return
304 }
305 return (n - e*t1) + cX, (m - d*t1) + cY, true
306}