/pliant/math/curve.pli
1	# Copyright  Hubert Tonneau  hubert.tonneau@pliant.cx
2	#
3	# This program is free software; you can redistribute it and/or
4	# modify it under the terms of the GNU General Public License version 2
5	# as published by the Free Software Foundation.
6	#
7	# This program is distributed in the hope that it will be useful,
8	# but WITHOUT ANY WARRANTY; without even the implied warranty of
9	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10	# GNU General Public License for more details.
11	#
12	# You should have received a copy of the GNU General Public License
13	# version 2 along with this program; if not, write to the Free Software
14	# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
15
16	module "/pliant/language/compiler.pli"
17	module "/pliant/language/context.pli"
18	module "/pliant/math/functions.pli"
19	module "/pliant/math/point.pli"
20	module "/pliant/math/transform.pli"
21
22	constant sqrt2 2^0.5
23	constant small 1/4
24	constant big 3/4
25
26
27
28	public
29
30	type CurvePoint
31	field Float in_x in_y
32	field Float x y
33	field Float out_x out_y
34	field uInt16 flags ; field uInt8 in_mode out_mode
35
36	constant curvepoint_automatic 0
37	constant curvepoint_directed 1
38	constant curvepoint_manual 2
39
40	constant curvepoint_selected 1
41	constant curvepoint_withinpoint 100h
42	constant curvepoint_withoutpoint 200h
43
44	constant automatic curvepoint_automatic
45	constant directed curvepoint_directed
46	constant manual curvepoint_manual
47
48
49
50	module "curve/mode.pli"
51	constant is_ready (cast 80h Int)
52
53	type Curve
54	field Array:CurvePoint pts
55	field Int status <- standard
56
57	function 'cast Status' c -> s
58	arg Curve c ; arg Status s
59	explicit
60	s := shunt (c:status .and. is_ready)<>0 success failure
61
62
63
64
65
66	function hsin x -> r
67	arg Float x r
68	r := sin x*pi/2
69
70	function hcos x -> r
71	arg Float x r
72	r := cos x*pi/2
73
74	function weigh x -> r
75	arg Float x r
76	r := 1-(cos x*pi/2)
77
78	function norme x y -> n
79	arg Float x y n
80	n := (x*x+y*y)^0.5
81
82	function scalarproduct x0 y0 x1 y1 -> p
83	arg Float x0 y0 x1 y1 p
84	p := x0*x1+y0*y1
85
86	function angle x y -> a
87	arg Float x y a
88	var Float xx yy
89	if x>=0
90	if y>=0
91	xx := x ; yy := y ; a := 0
92	else
93	xx := -y ; yy := x ; a := 1.5*pi
94	else
95	if y>=0
96	xx := y ; yy := -x ; a := 0.5*pi
97	else
98	xx := -x ; yy := -y ; a := pi
99	check xx>=0 and yy>=0
100	if yy<=xx
101	if xx=0
102	a := 0
103	else
104	a += atan yy/xx
105	else
106	a += 0.5*pi-(atan xx/yy)
107
108	function min a b c -> r
109	arg Float a b c r
110	r := min (min a b) c
111
112	function max a b c -> r
113	arg Float a b c r
114	r := max (max a b) c
115
116	function bound a m M -> r
117	arg Int a m M r
118	r := shunt a<=m m a>=M M a
119
120	function distance a b -> d
121	arg Point2 a b ; arg Float d
122	d := ( (a:x-b:x)*(a:x-b:x) + (a:y-b:y)*(a:y-b:y) )^0.5
123	if d=undefined
124	d := 0
125
126
127
128
129
130	function curve_point x y angle -> p
131	arg Float x y ; arg CBool angle ; arg CurvePoint p
132	p x := x ; p y := y
133	if angle
134	p in_x := 0 ; p in_y  := 0
135	p out_x := 0 ; p out_y := 0
136	p in_mode := manual ; p out_mode := manual
137	else
138	p in_mode := automatic ; p out_mode := automatic
139	p flags := 0
140
141	method p tg_in x y
142	arg_rw CurvePoint p ; arg Float x y
143	p in_x := x ; p in_y := y
144	p in_mode := directed
145
146	method p tg_out x y
147	arg_rw CurvePoint p ; arg Float x y
148	p out_x := x ; p out_y := y
149	p out_mode := directed
150
151	method p in x y
152	arg_rw CurvePoint p ; arg Float x y
153	p in_x := x ; p in_y := y
154	p in_mode := manual
155
156	method p out x y
157	arg_rw CurvePoint p ; arg Float x y
158	p out_x := x ; p out_y := y
159	p out_mode := manual
160
161	method p is_angle -> c
162	arg CurvePoint p ; arg CBool c
163	c := p:in_mode<>automatic or p:out_mode<>automatic
164
165	export curve_point '. tg_in' '. tg_out' '. in' '. out' '. is_angle'
166
167
168
169
170
171	method c reset
172	arg_rw Curve c
173	c:pts size := 0
174	c status := standard
175
176	method c through x y
177	arg_rw Curve c ; arg Float x y
178	c pts += curve_point x y false
179	if pliant_debugging_level>=2
180	c status := standard
181
182	method c angle x y
183	arg_rw Curve c ; arg Float x y
184	c pts += curve_point x y true
185	if pliant_debugging_level>=2
186	c status := standard
187
188	function '+=' c p
189	arg_rw Curve c ; arg CurvePoint p
190	c pts += p
191	if pliant_debugging_level>=2
192	c status := standard
193
194
195	method c compute mode
196	arg_rw Curve c ; arg Int mode
197	if c:pts:size<2
198	return
199	for (var Int k) 0 c:pts:size-1
200	var Pointer:CurvePoint p0 :> c:pts (shunt k<>0 k-1 c:pts:size-1)
201	var Pointer:CurvePoint p :> c:pts k
202	var Pointer:CurvePoint p1 :> c:pts (shunt k<>c:pts:size-1 k+1 0)
203	if p:in_mode<>manual or p:out_mode<>manual
204	var Float a0 := angle p:x-p0:x p:y-p0:y ; check a0=defined
205	var Float a1 := angle p1:x-p:x p1:y-p:y ; check a1=defined
206	var Float a := min (abs a1-a0-2.0*pi) (abs a1-a0) (abs a1-a0+2.0*pi)
207	var Float reduce := shunt a<=pi/2 1 sin:a
208	if p:in_mode=automatic or p:out_mode=automatic
209	var Float a := (a0+a1)/2
210	var Float tx := cos a
211	var Float ty := sin a
212	if (scalarproduct tx ty p1:x-p0:x p1:y-p0:y)<0 { tx := -tx ; ty := -ty }
213	if p:in_mode=automatic { p in_x := -tx ; p in_y := -ty }
214	if p:out_mode=automatic { p out_x := tx ; p out_y := ty }
215	if p:in_mode<>manual
216	var Float a0 := angle p0:x-p:x p0:y-p:y
217	var Float a1 := angle p:in_x p:in_y
218	var Float a := min (abs a1-a0-2*pi) (abs a1-a0) (abs a1-a0+2*pi)
219	var Float n := (1-cos:a)*(norme p0:x-p:x p0:y-p:y)
220	var Float d := (sqrt2-1) * sin:a * 2*sin:a * (norme p:in_x p:in_y)
221	var Float f := (shunt n<>0 and d<>0 n/d 0)*reduce;
222	p in_x *= f ; p in_y *=f
223	if p:out_mode<>manual
224	var Float a0 := angle p1:x-p:x p1:y-p:y
225	var Float a1 := angle p:out_x p:out_y
226	var Float a := min (abs a1-a0-2*pi) (abs a1-a0) (abs a1-a0+2*pi)
227	var Float n := (1-cos:a) * (norme p1:x-p:x p1:y-p:y)
228	var Float d := (sqrt2-1) * sin:a * 2*sin:a * (norme p:out_x p:out_y)
229	var Float f := (shunt n<>0 and d<>0 n/d 0)*reduce;
230	p out_x *= f ; p out_y *= f
231	c status := mode .or. is_ready
232	if pliant_debugging_level>=2
233	for (var Int k) 0 c:pts:size-1
234	var Pointer:CurvePoint p :> c:pts k
235	check p:x=defined and p:y=defined
236	check p:in_x=defined and p:in_y=defined
237	check p:out_x=defined and p:out_y=defined
238	(the_function '. compute' Curve Int) extra_module :> the_module "/pliant/math/curve/mode.pli"
239
240
241	method c size -> s
242	arg Curve c ; arg Int s
243	s := c:pts size
244
245	method c 'size :=' s
246	arg_rw Curve c ; arg Int s
247	c:pts size := s
248	c status := c:status .and. .not. is_ready
249
250	method c point i -> p
251	arg Curve c ; arg Int i ; arg_R CurvePoint p
252	p :> c:pts i
253
254	method c 'point :=' i p
255	arg_rw Curve c ; arg Int i ; arg CurvePoint p
256	c:pts i := p
257	c status := c:status .and. .not. is_ready
258
259	method c mode -> m
260	arg Curve c ; arg Int m
261	m := c:status .and. .not. (cast is_ready Int)
262
263	function compute_maximum -> m
264	arg Float m
265	m := 1
266	for (var Int step) 0 100
267	var Float f := 1+(exp -step)
268	while f<>1 and f*m=defined and -f*m=defined
269	m *= f
270	constant maximum compute_maximum
271
272	method c bbox x0 y0 x1 y1
273	arg Curve c ; arg_w Float x0 y0 x1 y1
274	check (c:status .and. is_ready)<>0
275	x0 := maximum ; y0 := maximum
276	x1 := -maximum ; y1 := -maximum
277	for (var Int i) 0 c:pts:size-1
278	var Pointer:CurvePoint p :> c:pts i
279	x0 := min x0 (min p:x p:x+p:in_x p:x+p:out_x)
280	y0 := min y0 (min p:y p:y+p:in_y p:y+p:out_y)
281	x1 := max x1 (max p:x p:x+p:in_x p:x+p:out_x)
282	y1 := max y1 (max p:y p:y+p:in_y p:y+p:out_y)
283	if x0=maximum
284	x0 := undefined ; y0 := undefined ; x1 := undefined ; y1 := undefined
285
286
287	method c pos param0 -> point
288	arg Curve c ; arg Float param0 ; arg Point2 point
289	check (c:status .and. is_ready)<>0
290	var CBool ol := (c:status .and. outline)<>0
291	var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
292	var Int i := bound (cast param-0.5 Int) 0 (shunt ol c:pts:size-1 c:pts:size-2)
293	param -= i
294	var Pointer:CurvePoint p0 :> c:pts i
295	var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
296	if (c:status .and. bezier)<>0
297	var Float x0 := p0 x
298	var Float x1 := p0:out_x+p0:x
299	var Float x2 := p1:in_x+p1:x
300	var Float x3 := p1 x
301	point x := (x3-3*x2+3*x1-x0)*param*param*param+3*(x2-2*x1+x0)*param*param+3*(x1-x0)*param+x0
302	var Float y0 := p0 y
303	var Float y1 := p0:out_y+p0:y
304	var Float y2 := p1:in_y+p1:y
305	var Float y3 := p1 y
306	point y := (y3-3*y2+3*y1-y0)*param*param*param+3*(y2-2*y1+y0)*param*param+3*(y1-y0)*param+y0
307	eif (c:status .and. conic)<>0
308	var Float x0 := p0 x
309	var Float x1 := ( p0:out_x+p0:x + p1:in_x+p1:x )/2
310	var Float x2 := p1 x
311	point x := (1-param)*(1-param)*x0 + 2*param*(1-param)*x1 + param*param*x2
312	var Float y0 := p0 y
313	var Float y1 := (p0:out_y+p0:y + p1:in_y+p1:y )/2
314	var Float y2 := p1 y
315	point y := (1-param)*(1-param)*y0 + 2*param*(1-param)*y1 + param*param*y2
316	else
317	var Float w0 w1
318	if c:pts:size=2 and not ol
319	w0 := 0 ; w1 := 0
320	eif i=0 and not ol and c:pts:0:out_mode=automatic
321	w0 := 0 ; w1 := param
322	eif i=c:pts:size-2 and not ol and c:pts:(c:pts:size-1):in_mode=automatic
323	w0 := 1-param ; w1 := 0
324	else
325	w0 := weigh 1-param ; w1 := weigh param
326	point x := w0*(p0:x+param*sqrt2*p0:out_x) + w1*(p1:x+(1-param)*sqrt2*p1:in_x) + (1-w0-w1)*((1-param)*p0:x+param*p1:x)
327	point y := w0*(p0:y+param*sqrt2*p0:out_y) + w1*(p1:y+(1-param)*sqrt2*p1:in_y) + (1-w0-w1)*((1-param)*p0:y+param*p1:y)
328	check point:x=defined and point:y=defined
329	check point:x=defined and point:y=defined
330
331
332
333
334
335	method c x_param param0 -> x
336	arg Curve c ; arg Float param0 x
337	check (c:status .and. is_ready)<>0
338	var CBool ol := (c:status .and. outline)<>0
339	var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
340	var Int i := bound (cast param-0.5 Int) 0 (shunt ol c:pts:size-1 c:pts:size-2)
341	param -= i
342	var Float w0 w1
343	if c:pts:size=2 and not ol
344	w0 := 0 ; w1 := 0
345	eif i=0 and not ol and c:pts:0:out_mode=automatic
346	w0 := 0 ; w1 := param
347	eif i=c:pts:size-2 and not ol and c:pts:(c:pts:size-1):in_mode=automatic
348	w0 := 1-param ; w1 := 0
349	else
350	w0 := weigh 1-param ; w1 := weigh param
351	var Pointer:CurvePoint p0 :> c:pts i
352	var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
353	x := w0*(p0:x+param*sqrt2*p0:out_x) + w1*(p1:x+(1-param)*sqrt2*p1:in_x) + (1-w0-w1)*((1-param)*p0:x+param*p1:x)
354	check x=defined
355
356	method c y_param param0 -> y
357	arg Curve c ; arg Float param0 y
358	check (c:status .and. is_ready)<>0
359	var CBool ol := (c:status .and. outline)<>0
360	var Float param := param0 * (shunt ol c:pts:size c:pts:size-1)
361	var Int i := bound (cast param-0.5 Int) 0 (shunt ol c:pts:size-1 c:pts:size-2)
362	param -= i
363	var Float w0 w1
364	if c:pts:size=2 and not ol
365	w0 := 0 ; w1 := 0
366	eif i=0 and not ol and c:pts:0:out_mode=automatic
367	w0 := 0 ; w1 := param
368	eif i=c:pts:size-2 and not ol and c:pts:(c:pts:size-1):in_mode=automatic
369	w0 := 1-param ; w1 := 0
370	else
371	w0 := weigh 1-param ; w1 := weigh param
372	var Pointer:CurvePoint p0 :> c:pts i
373	var Pointer:CurvePoint p1 :> c:pts (i+1)%c:pts:size
374	y := w0*(p0:y+param*sqrt2*p0:out_y) + w1*(p1:y+(1-param)*sqrt2*p1:in_y) + (1-w0-w1)*((1-param)*p0:y+param*p1:y)
375	check y=defined
376
377
378	method c param_x x deltax -> param
379	arg Curve c ; arg Float x deltax param
380	check (c:status .and. is_ready)<>0
381	var Float p0 := 0 ; var Float x0 := c x_param p0 ; if x0=x return:p0
382	var Float p1 := 1 ; var Float x1 := c x_param p1 ; if x1=x return:p1
383	if x0>x1  { swap x0 x1 ; swap p0 p1 }
384	if x<x0 return:p0
385	if x>x1 return:p1
386	var Int mode := 0
387	var Float fm pm xm
388	while true
389	if mode=0
390	fm := (x-x0)/(x1-x0)
391	pm := p0+fm*(p1-p0)
392	eif  mode=1
393	fm := (x-x0)/(x1-x0)
394	fm := shunt fm>=0.5 0.5+(fm-0.5)*(fm-0.5)*2 0.5-(0.5-fm)*(0.5-fm)*2
395	pm := p0+fm*(p1-p0)
396	else
397	pm := (p0+p1)/2
398	check pm=defined
399	if pm=p0 or pm=p1
400	return pm
401	check (p1>p0 and pm>p0 and pm<p1) or (p1<p0 and pm>p1 and pm<p0)
402	xm := c x_param pm
403	if (abs xm-x)<=deltax
404	return pm
405	eif xm<x
406	x0 := xm ; p0 := pm
407	mode := shunt fm<=small (min mode+1 2) (max mode-1 0)
408	eif xm>x
409	x1 := xm ; p1 := pm
410	mode := shunt fm>big (min mode+1 2) (max mode-1 0)
411	else
412	return pm
413
414	method c param_y y deltay -> param
415	arg Curve c ; arg Float y deltay param
416	check (c:status .and. is_ready)<>0
417	var Float p0 := 0 ; var Float y0 := c y_param p0 ; if y0=y return:p0
418	var Float p1 := 1 ; var Float y1 := c y_param p1 ; if y1=y return:p1
419	if y0>y1  { swap y0 y1 ; swap p0 p1 }
420	if y<y0 return:p0
421	if y>y1 return:p1
422	var Int mode := 0
423	var Float fm pm ym
424	while true
425	if mode=0
426	fm := (y-y0)/(y1-y0)
427	pm := p0+fm*(p1-p0)
428	eif  mode=1
429	fm := (y-y0)/(y1-y0)
430	fm := shunt fm>=0.5 0.5+(fm-0.5)*(fm-0.5)*2 0.5-(0.5-fm)*(0.5-fm)*2
431	pm := p0+fm*(p1-p0)
432	else
433	pm := (p0+p1)/2
434	check pm=defined
435	if pm=p0 or pm=p1
436	return pm
437	check (p1>p0 and pm>p0 and pm<p1) or (p1<p0 and pm>p1 and pm<p0)
438	ym := c y_param pm
439	if (abs ym-y)<=deltay
440	return pm
441	eif ym<y
442	y0 := ym ; p0 := pm
443	mode := shunt fm<=small (min mode+1 2) (max mode-1 0)
444	eif ym>y
445	y1 := ym ; p1 := pm
446	mode := shunt fm>big (min mode+1 2) (max mode-1 0)
447	else
448	return pm
449
450
451	method c y x deltax -> y
452	arg Curve c ; arg Float x deltax y
453	check (c:status .and. is_ready)<>0 and (c:status .and. y_from_x)<>0
454	var Float p := c param_x x deltax
455	y := c y_param p
456
457	method c x y deltay -> x
458	arg Curve c ; arg Float y deltay x
459	check (c:status .and. is_ready)<>0 and (c:status .and. x_from_y)<>0
460	var Float p := c param_y y deltay
461	x := c x_param p
462
463
464
465
466
467	method c polyline table t0 p0 t1 p1 m epsilon
468	arg Curve c ; arg_rw Array:Point2 table ; arg Float t0 t1 epsilon ; arg Transform2 m ; arg Point2 p0 p1
469	var Float t05 := (t0+t1)/2.0
470	var Point2 p05 := m (c pos t05)
471	if (distance (point (p0:x+p1:x)*0.5 (p0:y+p1:y)*0.5) p05)<=epsilon
472	table += p05 ; table += p1
473	else
474	c polyline table t0 p0 t05 p05 m epsilon
475	c polyline table t05 p05 t1 p1 m epsilon
476
477	method c polyline m epsilon -> table
478	arg Curve c ; arg Transform2 m ; arg Float epsilon ; arg Array:Point2 table
479	check (c:status .and. is_ready)<>0
480	check epsilon>0
481	table size := 1 ; table 0 := m (c pos 0)
482	var Int n := shunt (c:status .and. outline)<>0 c:pts:size c:pts:size-1
483	for (var Int i) 0 n-1
484	var Pointer:CurvePoint cp0 :> c:pts i
485	var Pointer:CurvePoint cp1 :> c:pts (i+1)%c:pts:size
486	if cp0:out_x=0 and cp0:out_y=0 and cp1:in_x=0 and cp1:in_y=0
487	table += m (point cp1:x cp1:y)
488	else
489	var Float t0 := i/n ; var Point2 p0 := table table:size-1
490	var Float t05 := (i+0.5)/n ; var Point2 p05 := m (c pos t05)
491	var Float t1 := (i+1)/n ; var Point2 p1 := m (c pos t1)
492	c polyline table  t0 p0  t05 p05 m epsilon
493	c polyline table  t05 p05 t1 p1 m epsilon
494
495
496	method c polyline table p0 p1 m epsilon
497	arg Curve c ; arg_rw Array:Point3 table ; arg Transform2 m ; arg Float epsilon ; arg Point3 p0 p1
498	var Point3 p05
499	p05 z := (p0:z+p1:z)/2.0
500	addressof:p05 map Point2 := m (c pos p05:z)
501	if (distance (point (p0:x+p1:x)*0.5 (p0:y+p1:y)*0.5) (addressof:p05 map Point2))<=epsilon
502	table += p05 ; table+=p1
503	else
504	c polyline table p0 p05 m epsilon
505	c polyline table p05 p1 m epsilon
506
507	method c polyline3 m epsilon -> table
508	arg Curve c ; arg Transform2 m ; arg Float epsilon ; arg Array:Point3 table
509	check (c:status .and. is_ready)<>0
510	check epsilon>0
511	table size := 1 ; table:0 z := 0 ; (addressof table:0) map Point2 := m (c pos 0)
512	var Int n := shunt (c:status .and. outline)<>0 c:pts:size c:pts:size-1
513	for (var Int i) 0 n-1
514	var Pointer:CurvePoint cp0 :> c:pts i
515	var Pointer:CurvePoint cp1 :> c:pts (i+1)%c:pts:size
516	if cp0:out_x=0 and cp0:out_y=0 and cp1:in_x=0 and cp1:in_y=0
517	var Point3 p1 ; p1 z := (i+1)/n ; addressof:p1 map Point2 := m (point cp1:x cp1:y)
518	table += p1
519	else
520	var Point3 p0 := table table:size-1
521	var Point3 p05 ; p05 z := (i+0.5)/n ; addressof:p05 map Point2 := m (c pos p05:z)
522	var Point3 p1 ; p1 z := (i+1)/n ; addressof:p1 map Point2 := m (c pos p1:z)
523	c polyline table p0 p05 m epsilon
524	c polyline table p05 p1 m epsilon
525
526
527	export Curve '. reset' '. through' '. angle' '+=' '. compute'
528	export '. size' '. size :=' '. point' '. point :=' '. mode' '. bbox' '. pos'
529	export '. x_param' '. y_param' '. param_x' '. param_y' '. x' '. y'
530	export '. polyline' '. polyline3'
531	export 'cast Status'
532
533