🚀 go-pugleaf

Talking about functional bezier curves, the
de Casteljau Algorithm provides a very efficient
and elegant way to calculate points on a bezier curve
(I suppose that Adobe implemented it's bezier curves
using this algorithm).
The de Casteljau scheme uses repeated linear interpolation.
If you have 4 points A,B,C and D with coordinates
(Ax,Ay) through (Dx,Dy) (that's why it's a functional curve,
usually one has parametric curves with the parameter interval
[0;1], i.e. the function b(t) return a point (x,y) for every
value t in [0;1] with the end conditions b(0) = A, b(1) = D;
but I guess you want to handle functional curves so let's see
how they can be handled...)

Suppose your bezier points look like this:


^ y   
|           B                       C 
|           o-----------------------o
|          /                         \
|         /                           \      
|        /                             \
|       /                               o D
|      /                                 
|     /
|    o A
|    
+-----*----------X-----------------------*----> x
    
You'd like to know the y value for a special X within the interval [Ax,Dx]
(If you compute the y value outside this range the algorithm gets numerically
unstable).
To do this, you simply compute the value X/(Dx-Ax), i.e. the relative position
of your desired x-value within the given interval. Then you create three new
bezier points from the given four, i.e. for every line segment (AB,BC and CD)
you create a new point that has the same relative position within the line
segment the X value has within the regarded x interval.

As an example, suppose that X is approximately like in the 'drawing', i.e.

X/(Dx-Ax) = 1/3 =: b

Then you compute the new points A1, B1 and C1 as barycentric combinations
of the original points:

A1 = (1 - b) * A   + b * B    in vector notation : A1 = A + 1/3 ( B - A )
B1 = (1 - b) * B   + b * C
C1 = (1 - b) * C   + b * D


result (ASCII characters are not that handy for bezier curve drawings) :

^ y   
|                 B1                    
|                 o                    
|                                    
|                                     o C1
|                                     
|                                       
|       o A1                                 
|    
|    
|    
+-----*----------X-----------------------*----> x


If you repeat this so called de Casteljau step again, you get two points
A2 and B2, located approximately as follows:




^ y   
|                                       
|                                      
|                       o B2                     
|                                          
|           o A2                                  
|                                       
|                                            
|    
|    
|    
+-----*----------X-----------------------*----> x



and doing it again you get what you wanted :


^ y   
|                                       
|                                      
|                                                
|                o A3 = b(X)             
|                                                 
|                                       
|                                            
|    
|    
|    
+-----*----------X-----------------------*----> x



I supply a small postscript fragment that expects an array of
bezier points and an x value and then computes the point on 
the bezier curve for you using the de Casteljau scheme:

hope that helps,

		Patsch

--------------------------------------------------------------------
   "We are sorry, but the number you have dialed is imaginary."
       "Please rotate your phone 90 degrees and try again."
--------------------------------------------------------------------


%%----------------------------- CUT HERE -----------------------------------
%!

%
% [ x0 y0  ...  xN yN  ] x-value calcbez y-value
%
% calcbez : computes a point on a functional bezier curve
%

/calcbez
{
    exch	    	    	    	% get array
    dup length /anzbez exch def		% length of array
    dup 0 get		    		% x value of 'leftmost' point
    exch dup anzbez 2 sub get		% and of rightmost point
    3 -1 roll sub 
    3 -1 roll
    exch div   	    	    	    	% compute weight
    casteljau				% call de Casteljau algo
} def

%
% [ x0 y0  ...  xN yN  ] weight casteljau  x y 
%
% casteljau : computes a point on a bezier curve using repeated linear int.
%

/s 100 string def

/casteljau
{
    /bary exch def	    	    	% get weight
    /fld  exch def			% and array of points
    fld length /anzbez exch def		% #  of points

    %% (\012\012) print

    anzbez 2 eq
    {
	fld aload pop	    	    	% recursion ended
    }
    {
	[   	    	    	    	% push mark
	    0 2 anzbez 3 sub    	% loop
	    {
		/indx exch def	    	% current array index
		fld indx get
		1 bary sub mul 	    	% compute barycentric combination
		fld indx 2 add get
		bary mul add   	    	% new x value
		fld indx 1 add get
		1 bary sub mul
		fld indx 3 add get
		bary mul add   	    	% new y value
		%% (\012INDX ) print indx s cvs print 
		%% (  X = ) print exch dup s cvs print
		%% ( , Y = ) print exch dup s cvs print
	    } for
	]     	    	    	    	% push end of array mark
	bary casteljau			% do recursion
    } ifelse
} def

/showresults
{
    (\012\012Result : X = ) print
    exch 100 string cvs print
    ( , Y = ) print 
    100 string cvs print (\012) print
} def

%%
%% MAIN
%%

[ 0 0  1 4  2 4  3 0  ] 0 calcbez showresults %% change your parameters here
[ 0 0  1 4  2 4  3 0  ] 1 calcbez showresults %% change your parameters here
[ 0 0  1 4  2 4  3 0  ] 2 calcbez showresults %% change your parameters here
[ 0 0  1 4  2 4  3 0  ] 3 calcbez showresults %% change your parameters here

%%----------------------------- CUT HERE -----------------------------------