/* * Recommended reading: * * Wavelet Rasterization *                 (2011) J. Manson, S.Schaefer * * A Simple and Fast Line-Clipping Method as a Scratch Extension for Computer Graphics Education *                 (2019) Dimitrios Matthes, Vasileios Drakopoulos * * Clipping of Bézier Curves *                 (1985) Fuhua Cheng, Chi-Cheng Lin * * On the numerical condition of polynomials in Bernstein form *                 (1987) R. T. Farouki, V. T. Rajan * * Numerical Condition of Polynomials in Different Forms *                 (2001) Hong Zhang * * Computer Aided Geometric Design *                 (2017) Thomas W. Sederberg */#include "otfsys.h"#include "otf.h"#include "otfpriv.h" #pragma clang float_control(precise, on) typedef struct SegQ SegQ;typedef struct Spt Spt;typedef double Sval; struct Spt {	Sval x;	Sval y;}; struct SegQ {	union {		struct {			union {				Spt p0;				Sval v0[2];			};			union {				Spt p1;				Sval v1[2];			};			union {				Spt p2;				Sval v2[2];			};		};		Sval v[3*2];	};}; /* the divisions when calculating coefficients are reduced, * thus pixel value range is [0.0-QBZR_PIX_SCALE] instead of [0.0-1.0] */#define QBZR_PIX_SCALE 4.0 #define MAXCOMPONENTS 16#define ε 1.0e-9#define is₀(v) ((v) > -ε && (v) < ε) static intclosest²(int v){	v--;	v |= v>>1;	v |= v>>2;	v |= v>>4;	v |= v>>8;	v |= v>>16;	return v + 1;} static Sval *qslv(Sval *qs, Sval a, Sval b, Sval c){	Sval p, q, d, t;	int i; #define ins(q) \	if(q > 0 && q < 1){ \		i = -1; \		if((t = qs[i]) < q){ \			*qs++ = q; \		}else if(q < t){ \			qs++; \			do{ \				qs[i--] = t; \				qs[i] = q; \				t = qs[i-1]; \			}while(q < t); \		} \	} 	if(unlikely(is₀(a))){		if(!is₀(b)){			q = -c/b;			ins(q);		}	}else{		p = b/(2.0*a);		d = p*p - c/a;		if(d > ε){			d = sqrt(d);			q = -d-p;			ins(q);			q = d-p;			ins(q);		}	}#undef ins	return qs;} static intqKL(SegQ *s₀, int jj, Sval px, Sval py, Sval *K, Sval *L){	Sval *q, *e, qs[1+1+4*2+1];	Sval a, b, c, w, v;	int i, n, r;	SegQ s, z; 	/* transform */	s.p0.x = s₀->p0.x*jj - px;	s.p0.y = s₀->p0.y*jj - py;	s.p1.x = s₀->p1.x*jj - px;	s.p1.y = s₀->p1.y*jj - py;	s.p2.x = s₀->p2.x*jj - px;	s.p2.y = s₀->p2.y*jj - py; #define e(t,a) (s.p0.a*(1-t)*(1-t) + 2*s.p1.a*(1-t)*t + s.p2.a*t*t)#define within(v) ((w = e(v, x)) >= -ε && w <= 1+ε && (w = e(v, y)) >= -ε && w <= 1+ε) 	/* clip to quadtree cell */	q = qs;	*q++ = -1;	if(s.p0.x >= 0 && s.p0.x <= 1 && s.p0.y >= 0 && s.p0.y <= 1)		*q++ = 0;	c = s.p0.x;	a = c - 2*s.p1.x + s.p2.x;	b = 2*(s.p1.x - c);	q = qslv(q, a, b, c);	q = qslv(q, a, b, c-1);	c = s.p0.y;	a = c - 2*s.p1.y + s.p2.y;	b = 2*(s.p1.y - c);	q = qslv(q, a, b, c);	q = qslv(q, a, b, c-1);	if(s.p2.x >= 0 && s.p2.x <= 1 && s.p2.y >= 0 && s.p2.y <= 1)		*q++ = 1; 	n = 1;	for(e = qs+1; e < q; e++){		v = *e;		if(within(v))			qs[n++] = v;	}#undef within 	/* calculate K & L on subsections */	r = 0;	for(i = 1; i < n-1; i++){		a = qs[i];		b = qs[i+1]; 		z.p0.y = e(a, y);		z.p2.y = e(b, y);		if(is₀(z.p0.y) && is₀(z.p2.y))			continue;		z.p0.x = e(a, x);		z.p2.x = e(b, x);		if(is₀(z.p0.x) && is₀(z.p2.x))			continue;#undef e 		c = 1 - a;		v = (b - a)/c;		z.p1.x = (1-v)*z.p0.x + v*(c*s.p1.x + a*s.p2.x);		z.p1.y = (1-v)*z.p0.y + v*(c*s.p1.y + a*s.p2.y); 		K[0] += z.p2.y - z.p0.y;		K[1] += z.p0.x - z.p2.x;		L[0] += 3*z.p2.x*z.p2.y		      + 2*z.p2.y*z.p1.x		      - 2*z.p2.x*z.p1.y		      +   z.p2.y*z.p0.x		      + 2*z.p1.y*z.p0.x		      -   z.p0.y*(  z.p2.x + 2*z.p1.x + 3*z.p0.x);		L[1] += 2*z.p1.y*z.p0.x		      +   z.p2.y*(2*z.p1.x +   z.p0.x)		      - 2*z.p1.x*z.p0.y		      + 3*z.p0.x*z.p0.y		      -   z.p2.x*(3*z.p2.y + 2*z.p1.y +   z.p0.y);		r = 1;	}	return r;} static intlKL(SegQ *s₀, int jj, Sval px, Sval py, Sval *K, Sval *L){	Sval x₁, x₂, y₁, y₂, α₁, α₂, β₁, β₂; 	/* transform */	α₁ = x₁ = s₀->p0.x*jj - px;	β₁ = y₁ = s₀->p0.y*jj - py;	α₂ = x₂ = s₀->p1.x*jj - px;	β₂ = y₂ = s₀->p1.y*jj - py; 	if(α₁ < 0){		if(α₂ < 0) return 0;		α₁ = 0;		β₁ = (y₂-y₁)/(x₂-x₁)*(0-x₁)+y₁;	}else if(α₁ >= 1){		if(α₂ >= 1) return 0;		α₁ = 1;		β₁ = (y₂-y₁)/(x₂-x₁)*(1-x₁)+y₁;	}	if(α₂ < 0){		α₂ = 0;		β₂ = (y₂-y₁)/(x₂-x₁)*(0-x₁)+y₁;	}else if(α₂ > 1){		α₂ = 1;		β₂ = (y₂-y₁)/(x₂-x₁)*(1-x₁)+y₁;	}	if(β₁ < 0){		if(β₂ < 0) return 0;		β₁ = 0;		α₁ = (x₂-x₁)/(y₂-y₁)*(0-y₁)+x₁;	}else if(β₁ >= 1){		if(β₂ >= 1) return 0;		β₁ = 1;		α₁ = (x₂-x₁)/(y₂-y₁)*(1-y₁)+x₁;	}	if(β₂ < 0){		β₂ = 0;		α₂ = (x₂-x₁)/(y₂-y₁)*(0-y₁)+x₁;	}else if(β₂ > 1){		β₂ = 1;		α₂ = (x₂-x₁)/(y₂-y₁)*(1-y₁)+x₁;	} 	K[0] += β₂ - β₁;	K[1] += α₁ - α₂;	L[0] += 3.0*(β₂ - β₁)*(α₁ + α₂);	L[1] += 3.0*(α₁ - α₂)*(β₁ + β₂);	return 1;} static u64intCxy(SegQ *s, int ns, int jj, int px, int py, Sval *c, u64int *m, u64int tm){	int (*f)(SegQ*, int, Sval, Sval, Sval*, Sval*);	Sval K[4][2], L[4][2], x₀, x₁, x₂, y₀, y₁, y₂, psz;	u64int tx₀, tx₁, z, all;	u8int w;	int i; 	jj *= 2;	px *= 2;	py *= 2;	if(jj > 8){		tx₀ = 0;		tx₁ = 0;	}else{		tx₀ = 1ULL<<(px*jj + py);		tx₁ = 1ULL<<((px+1)*jj + py);	}	psz = 1.0/(Sval)jj;	x₀ = psz*px;	y₀ = psz*py;	x₁ = psz+x₀;	y₁ = psz+y₀;	x₂ = psz+x₁;	y₂ = psz+y₁;	all = 0;	for(i = 0; i < ns; i++, s++){		if((m[i] & tm) == 0 ||		   s->p0.x >= x₂ && s->p1.x >= x₂ && s->p2.x >= x₂ ||		   s->p0.x <= x₀ && s->p1.x <= x₀ && s->p2.x <= x₀ ||		   s->p0.y >= y₂ && s->p1.y >= y₂ && s->p2.y >= y₂ ||		   s->p0.y <= y₀ && s->p1.y <= y₀ && s->p2.y <= y₀)			continue; 		K[0][0] = K[0][1] = 0;		K[1][0] = K[1][1] = 0;		K[2][0] = K[2][1] = 0;		K[3][0] = K[3][1] = 0;		L[0][0] = L[0][1] = 0;		L[1][0] = L[1][1] = 0;		L[2][0] = L[2][1] = 0;		L[3][0] = L[3][1] = 0;		z = 0;		f = s->p1.x == s->p2.x && s->p1.y == s->p2.y ? lKL : qKL;		w =			(s->p0.x <= x₁ || s->p1.x <= x₁ || s->p2.x <= x₁)<<0 |			(s->p0.x >= x₁ || s->p1.x >= x₁ || s->p2.x >= x₁)<<1 |			(s->p0.y <= y₁ || s->p1.y <= y₁ || s->p2.y <= y₁)<<2 |			(s->p0.y >= y₁ || s->p1.y >= y₁ || s->p2.y >= y₁)<<3; 		if(tx₀ == 0){			if((w & 5) == 5) z |= f(s, jj, px+0, py+0, K[0], L[0]);			if((w & 6) == 6) z |= f(s, jj, px+1, py+0, K[2], L[2]);			if((w & 9) == 9) z |= f(s, jj, px+0, py+1, K[1], L[1]);			if((w & 10) == 10) z |= f(s, jj, px+1, py+1, K[3], L[3]);		}else{			if((w & 5) == 5 && f(s, jj, px+0, py+0, K[0], L[0])) z |= tx₀;			if((w & 6) == 6 && f(s, jj, px+1, py+0, K[2], L[2])) z |= tx₁;			if((w & 9) == 9 && f(s, jj, px+0, py+1, K[1], L[1])) z |= tx₀<<1;			if((w & 10) == 10 && f(s, jj, px+1, py+1, K[3], L[3])) z |= tx₁<<1;			m[ns+i] |= z;			all |= z;		} 		if(z != 0){			c[0] += (L[0][1] - L[1][1] + L[2][1] - L[3][1])/6.0 + K[1][1] + K[3][1];			c[1] += (L[0][0] + L[1][0] - L[2][0] - L[3][0])/6.0 + K[2][0] + K[3][0];			c[2] += (L[0][0] - L[1][0] - L[2][0] + L[3][0])/6.0 + K[2][0] - K[3][0];		}	}	return all;} static intqbzr(SegQ *seg, int ns, int dw, int dh, int pitch, u8int *p){	int ju², j², w, h, i, j, ju, px, py, v;	Sval rju², s₀, s, sl, sq, *c, *c₀x, *c₀y, *cs, zx, zy;	u64int *ma, *mb, all, nall; 	ju² = closest²(dh);	if(ju² < (j = closest²(dw)))		ju² = j;	rju² = 1.0 / ju²;	for(ju = 1; (1<<ju) < ju²; ju++); 	/* calculate area */	sl = sq = 0;	for(i = 0; i < ns; i++){#define det(a,b) (a.x*b.y - a.y*b.x)#define cL(s) det(s.p0, s.p2)#define cQ(s) (2*det(s.p0, s.p1) + 2*det(s.p1, s.p2) + det(s.p0, s.p2))		if(seg[i].p1.x == seg[i].p2.x && seg[i].p1.y == seg[i].p2.y)			sl -= cL(seg[i]);		else			sq -= cQ(seg[i]);#undef det#undef cL#undef cQ		/* scale */		seg[i].v[0] *= rju²;		seg[i].v[1] *= rju²;		seg[i].v[2] *= rju²;		seg[i].v[3] *= rju²;		seg[i].v[4] *= rju²;		seg[i].v[5] *= rju²;	}	s₀ = (sl + sq/3)*QBZR_PIX_SCALE/(2*ju²*ju²); 	/* working space:	 * quick is-segment-in-the-cell tests (two 64-bit nums per segment)	 * wavelet coefficients (count by geometric series, 3 per "pixel" of every zoom level)	 */	j = 2*ns*sizeof(*ma) + (1-(4<<2*(ju-1)))/(1-4)*3*sizeof(*c);	/* current zoom level bitmasks for quick tests */	if((ma = calloc(1, j)) == nil){		werrstr("no memory");		return -1;	}	/* next zoom level bitmasks */	mb = ma + ns;	/* first (1x1) zoom level - all segments are there */	for(i = 0; i < ns; i++)		ma[i] = 1; 	/* coefficients */	c = cs = (Sval*)(ma + 2*ns); 	/* bitmasks used per-segment are united (bitwise OR) to skip entire empty cells */	all = 1;	nall = 0;	/* first few loops build the test bitmasks: 1x1 -> 4x4 -> 8x8 */	for(j = 0, j² = 1; j < ju && j <= 3; j++, j² <<= 1){		for(px = 0; px < j²; px++){			for(py = 0; py < j²; py++, c += 3){				u64int tm = 1ULL<<(px*j² + py);				if(all & tm)					nall |= Cxy(seg, ns, j², px, py, c, ma, tm);			}		}		if(j != 3){			for(i = 0; i < ns; i++){				ma[i] = mb[i];				mb[i] = 0;			}			all = nall;			nall = 0;		}	} 	/* no more bitmasks, just calculate coefficients */	for(; j < ju; j++, j² <<= 1){		for(px = 0; px < j²; px++){			/* bitmasks are expanded only up to 8x8 (64 bits)			 * so use that for more coarse testing			 */			u64int tm₀ = 1ULL<<((px>>(j-3))*8);			for(py = 0; py < j²; py++, c += 3){				u64int tm = tm₀ << (py>>(j-3));				if(all & tm)					Cxy(seg, ns, j², px, py, c, ma, tm);			}		}	} 	/* reverse Y axis (opentype coord vs plan 9) */	for(h = dh-1; h >= 0; h--){		for(w = 0; w < dw; w++){			c = cs;			s = s₀;			for(j² = 1; j² < ju²; j² *= 2){				zx = rju² * j² * w;				c₀x = c;				px = (int)zx;				zx -= px;				c += px*3*j²;				for(; zx >= 0; px++, zx--, c = c₀y + 3*j²){					c₀y = c;					zy = rju² * j² * h;					py = (int)zy;					zy -= py;					c += 3*py;					for(; zy >= 0; py++, zy--, c += 3){						switch((zy < 0.5)<<1 | (zx < 0.5)){						case 0:  s += +c[0] + c[1] - c[2]; break;						case 1:  s += +c[0] - c[1] + c[2]; break;						case 2:  s += -c[0] + c[1] + c[2]; break;						default: s += -c[0] - c[1] - c[2]; break;						}					}				}				c = c₀x + 3*j²*j²;			}			v = s/QBZR_PIX_SCALE*255;			v *= v > 0;			*p++ = v < 255 ? v : 255;		}		p += pitch - dw;	}	free(ma);	return 0;} static SegQ *addpoints(SegQ *s, Spt *pts₀, Glyf *g, Sval dx, Sval dy, Sval sx, Sval sy, Sval *bb){	int k, npts, xMin, yMin, xMax, yMax;	Point *p₀, *p, *prev;	Spt *e, *pts;	Sval v; 	xMin = yMin = 99999;	xMax = yMax = -99999; 	for(k = 0; k < g->numberOfContours; k++){		npts = g->simple->endPtsOfContours[k]+1;		if(k > 0)			npts -= g->simple->endPtsOfContours[k-1]+1;		if(npts < 2)			continue; 		e = pts = pts₀;		p₀ = p = g->simple->points + (k > 0 ? g->simple->endPtsOfContours[k-1]+1 : 0);		prev = p₀+npts-1;		if(!p->onCurve){			/* fun stuff */			if(prev->onCurve)				*e = (Spt){dx+sx*prev->x, dy+sy*prev->y};			else				*e = (Spt){dx+sx/2*(p->x+prev->x), dy+sy/2*(p->y+prev->y)};			e++;		}		for(prev = nil; p < p₀+npts; prev = p, p++, e++){			if(prev != nil && p->onCurve == prev->onCurve){				/* more fun stuff */				if(p->onCurve) /* straight line */					*e = (Spt){dx+sx*p->x, dy+sy*p->y};				else /* a point in the middle */					*e = (Spt){dx+sx/2*(p->x+prev->x), dy+sy/2*(p->y+prev->y)};				e++;			}			*e = (Spt){dx+sx*p->x, dy+sy*p->y}; 			if(bb != nil){				if(xMin > p->x) xMin = p->x;				if(yMin > p->y) yMin = p->y;

Lines 1–500 of 898. Load more