1 /* Copyright (C) 2019 by Skef Iterum */
2 /*
3  * Redistribution and use in source and binary forms, with or without
4  * modification, are permitted provided that the following conditions are met:
5 
6  * Redistributions of source code must retain the above copyright notice, this
7  * list of conditions and the following disclaimer.
8 
9  * Redistributions in binary form must reproduce the above copyright notice,
10  * this list of conditions and the following disclaimer in the documentation
11  * and/or other materials provided with the distribution.
12 
13  * The name of the author may not be used to endorse or promote products
14  * derived from this software without specific prior written permission.
15 
16  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
17  * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
18  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
19  * EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
20  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
21  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
22  * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
23  * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
24  * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
25  * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  */
27 
28 #include <fontforge-config.h>
29 
30 #include "utanvec.h"
31 #include "splineutil.h"
32 #include "splineutil2.h"
33 
34 #include <math.h>
35 #include <assert.h>
36 
37 #define UTMARGIN (1e-7)     // Arrived at through testing
38 #define UTSOFTMARGIN (1e-5) // Fallback margin for some cases
39 #define BPNEAR(bp1, bp2) BPWithin(bp1, bp2, UTMARGIN)
40 #define UTRETRY (1e-9)      // Amount of "t" to walk back from the
41                             // ends of degenerate splines to find
42                             // a slope.
43 
MakeUTanVec(bigreal x,bigreal y)44 BasePoint MakeUTanVec(bigreal x, bigreal y) {
45     BasePoint ret = { 0.0, 0.0 };
46     bigreal len = x*x + y*y;
47     if ( len!=0 ) {
48 	len = sqrt(len);
49 	ret.x = x/len;
50 	ret.y = y/len;
51     }
52     return ret;
53 }
54 
NormVec(BasePoint v)55 BasePoint NormVec(BasePoint v) {
56     return MakeUTanVec(v.x, v.y);
57 }
58 
59 /* Orders unit tangent vectors on an absolute -PI+e to PI
60  * equivalent basis
61  */
UTanVecGreater(BasePoint uta,BasePoint utb)62 int UTanVecGreater(BasePoint uta, BasePoint utb) {
63     assert(    RealNear(BPLenSq(uta),1)
64             && RealNear(BPLenSq(utb),1) );
65 
66     if (uta.y >= 0) {
67 	if (utb.y < 0)
68 	    return true;
69 	return uta.x < utb.x && !BPNEAR(uta, utb);
70     }
71     if (utb.y >= 0)
72 	return false;
73     return uta.x > utb.x && !BPNEAR(uta, utb);
74 }
75 
76 /* True if rotating from ut1 to ut3 in the specified direction
77  * the angle passes through ut2.
78  *
79  * Note: Returns true if ut1 is near ut2 but false if ut2 is
80  * near ut3
81  */
UTanVecsSequent(BasePoint ut1,BasePoint ut2,BasePoint ut3,int ccw)82 int UTanVecsSequent(BasePoint ut1, BasePoint ut2, BasePoint ut3,
83                            int ccw) {
84     BasePoint tmp;
85 
86     if ( BPNEAR(ut1, ut2) )
87 	return true;
88 
89     if ( BPNEAR(ut2, ut3) || BPNEAR(ut1, ut3) )
90 	return false;
91 
92     if (ccw) {
93 	tmp = ut1;
94 	ut1 = ut3;
95 	ut3 = tmp;
96     }
97 
98     if ( UTanVecGreater(ut1, ut3) ) {
99 	return UTanVecGreater(ut1, ut2) && UTanVecGreater(ut2, ut3);
100     } else {
101 	return    (UTanVecGreater(ut1, ut2) && UTanVecGreater(ut3, ut2))
102 	       || (UTanVecGreater(ut2, ut1) && UTanVecGreater(ut2, ut3));
103     }
104 }
105 
JointBendsCW(BasePoint ut_ref,BasePoint ut_vec)106 int JointBendsCW(BasePoint ut_ref, BasePoint ut_vec) {
107     // magnitude of cross product
108     bigreal r = ut_ref.x*ut_vec.y - ut_ref.y*ut_vec.x;
109     if ( RealWithin(r, 0.0, UTMARGIN) )
110 	return false;
111     return r > 0;
112 }
113 
SplineUTanVecAt(Spline * s,bigreal t)114 BasePoint SplineUTanVecAt(Spline *s, bigreal t) {
115     BasePoint raw;
116 
117     if ( SplineIsLinearish(s) ) {
118 	raw.x = s->to->me.x - s->from->me.x;
119 	raw.y = s->to->me.y - s->from->me.y;
120     } else {
121 	// If one control point is colocated with its on-curve point the
122 	// slope will be undefined at that end, so walk back a bit for
123 	// consistency
124 	if (   RealWithin(t, 0, UTRETRY)
125 	    && BPWithin(s->from->me, s->from->nextcp, 1e-13) )
126 	    t = UTRETRY;
127 	else if (   RealWithin(t, 1, UTRETRY)
128 	    && BPWithin(s->to->me, s->to->prevcp, 1e-13) )
129 	    t = 1-UTRETRY;
130 	raw = SPLINEPTANVAL(s, t);
131 	// For missing slopes take a very small step away and try again.
132 	if ( raw.x==0 && raw.y==0 )
133 	    raw = SPLINEPTANVAL(s, (t+UTRETRY <= 1.0) ? t+UTRETRY : t-UTRETRY);
134     }
135     return MakeUTanVec(raw.x, raw.y);
136 }
137 
138 /* Return the lowest t value greater than min_t such that the
139  * tangent at the point is parallel to ut, or -1 if the tangent
140  * never has that slope.
141  */
142 #define ROTY(p, ut) (-(ut).y*(p).x + (ut).x*(p).y)
143 #define ZCHECK(v) (RealWithin((v),0,1e-10) ? 0 : (v))
144 
SplineSolveForUTanVec(Spline * spl,BasePoint ut,bigreal min_t,bool picky)145 bigreal SplineSolveForUTanVec(Spline *spl, BasePoint ut, bigreal min_t,
146                               bool picky) {
147     bigreal tmp, yto;
148     extended te1, te2;
149     Spline1D ys1d;
150 
151     // Nothing to "solve" for; lines should be handled by different means
152     // because you can't "advance" along them based on slope changes
153     if ( SplineIsLinear(spl) )
154         return -1;
155 
156     // Only check t==0 if min_t is negative
157     if ( min_t<0 && BPNEAR(ut, SplineUTanVecAt(spl, 0.0)) )
158        return 0.0;
159 
160     // This is probably over-optimized
161     ys1d.d = ROTY(spl->from->me, ut);
162     yto = ROTY(spl->to->me, ut);
163     if ( spl->order2 ) {
164         ys1d.c = ZCHECK(2*(ROTY(spl->from->nextcp, ut)-ys1d.d));
165         ys1d.b = ZCHECK(yto-ys1d.d-ys1d.c);
166         ys1d.a = 0;
167     } else {
168 	tmp = ROTY(spl->from->nextcp, ut);
169 	ys1d.c = ZCHECK(3*(tmp-ys1d.d));
170 	ys1d.b = ZCHECK(3*(ROTY(spl->to->prevcp, ut)-tmp)-ys1d.c);
171 	ys1d.a = ZCHECK(yto-ys1d.d-ys1d.c-ys1d.b);
172     }
173     if ( ys1d.a!=0 && (   Within16RoundingErrors(ys1d.a+ys1d.d,ys1d.d)
174                        || Within16RoundingErrors(ys1d.a+yto,yto)))
175 	ys1d.a = 0;
176 
177     // After rotating by theta the desired angle will be at a y extrema
178     SplineFindExtrema(&ys1d, &te1, &te2);
179 
180     if (   te1 != -1 && te1 > (min_t+1e-9)
181         && BPNEAR(ut, SplineUTanVecAt(spl, te1)) )
182 	return te1;
183     if (   te2 != -1 && te2 > (min_t+1e-9)
184         && BPNEAR(ut, SplineUTanVecAt(spl, te2)) )
185 	return te2;
186 
187     // Check t==1
188     if ( (min_t+UTRETRY) < 1 && BPNEAR(ut, SplineUTanVecAt(spl, 1.0)) )
189         return 1.0;
190 
191     if ( picky )
192 	return -1.0;
193 
194     if (   te1 != -1 && te1 > (min_t+1e-9)
195         && BPWithin(ut, SplineUTanVecAt(spl, te1), UTSOFTMARGIN) )
196 	return te1;
197     if (   te2 != -1 && te2 > (min_t+1e-9)
198         && BPWithin(ut, SplineUTanVecAt(spl, te2), UTSOFTMARGIN) )
199 	return te2;
200 
201     // Check t==1
202     if (    (min_t+UTRETRY) < 1
203         && BPWithin(ut, SplineUTanVecAt(spl, 1.0), UTSOFTMARGIN) )
204         return 1.0;
205 
206     // Nothing found
207     return -1.0;
208 }
209 
210 #if 0
211 void UTanVecTests() {
212     int i, j, k, a, b;
213     bigreal r, x, y, z;
214     BasePoint ut[361], utmax = { -1, 0 }, utmin = UTMIN;
215 
216     for (i = 0; i<=360; i++) {
217 	r = ((180-i) * FF_PI / 180);
218 	ut[i] = MakeUTanVec(cos(r), sin(r));
219 	// printf("i:%d, i.x:%lf, i.y:%lf\n", i, ut[i].x, ut[i].y);
220     }
221 
222     for (i=0; i<360; ++i) {
223 	for (j = 0; j<360; ++j) {
224 	    if ( i==j )
225 		continue;
226 	    if ( i<j )
227 		assert(    UTanVecGreater(ut[i], ut[j])
228 		       && !UTanVecGreater(ut[j], ut[i]));
229 	    else
230 		assert(    UTanVecGreater(ut[j], ut[i])
231 		       	&& !UTanVecGreater(ut[i], ut[j]));
232 	    assert(!BPNEAR(ut[i], ut[j]));
233 	    x = atan2(ut[i].x, ut[i].y);
234 	    y = atan2(ut[j].x, ut[j].y);
235 	    z = x-y;
236 	    z = atan2(sin(z), cos(z));
237 	    if ( RealNear(z, FF_PI) || RealNear(z, -FF_PI) )
238 		continue;
239 	    assert( (z>0) == JointBendsCW(ut[i], ut[j]) );
240 	}
241 	assert( i==0 || UTanVecGreater(utmax, ut[i]) );
242 	assert( UTanVecGreater(ut[i], utmin) );
243     }
244 
245     for (i=0; i<360; ++i)
246 	for (j = 0; j<360; ++j)
247 	    for (k = 0; k<360; ++k) {
248 		a = UTanVecsSequent(ut[i], ut[j], ut[k], false);
249 		b = UTanVecsSequent(ut[i], ut[j], ut[k], true);
250 		if ( i==j )
251 		    assert(a && b);
252 		else if ( j==k || i==k )
253 		    assert(!a && !b);
254 		else {
255 		    x = atan2(ut[i].x, ut[i].y);
256 		    y = atan2(ut[j].x, ut[j].y);
257 		    z = atan2(ut[k].x, ut[k].y);
258 		    y = fmod(4*FF_PI+y-x,2*FF_PI);
259 		    z = fmod(4*FF_PI+z-x,2*FF_PI);
260 		    if ( y<z )
261 			assert( a && !b );
262 		    else
263 			assert( !a && b );
264 		}
265 	    }
266 }
267 #endif
268