Reptorian G'MIC Filters

You could approximate curve length with a series of line segments. Your (approximate) length is the sum of line segments drawn between curve plots.

Not exact, of course, but the approximation improves with increasing plot counts.

The rabbit to catch is a G’MIC-cal way to plot a Bézier curve. Here's how I would lay the trap.

TL;DR

  1. Usual Math Functions: lerp(a,b,t)→(1-t)*a+t*b, 15th bullet point down.
  2. De Casteljau's algorithm

Details
Apologies for the slovenly diagram.
post_01
A Lerp is a unit pyramid
The key idea is a linear interpolation between two points, samples, data items — at this level of generality, take your pick. This is a “lerp.” A lerp takes a weighted average of two given data items, a and b, at a parametric value of, say, t=0.6. A lerp of a and b is just ‘a*(1-0.6) + b*0.6’, or, using G’MIC's math function:

$ gmic a=27.6,34.2 b=9.875,4.75 param=0.6 aprime='{lerp([$a],[$b],$param)}'
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable 'a=27.6,34.2'.
[gmic]-0./ Set local variable 'b=9.875,4.75'.
[gmic]-0./ Set local variable 'param=0.6'.
[gmic]-0./ Set local variable 'aprime=16.965,16.530000000000001'.
[gmic]-0./ End G'MIC interpreter.

Diagrammatically, we might depict this as a an enactment of the “1st degree pyramid algorithm”, (see slovenly diagram), aka a lerp between points a and b and producing a weighted average point, a'. In numbers, the point a' is 3/5 of the way from point a, (x=27.6, y=34.2) to point b, (x=9.875,y=4.75), (x=16.965, y=16.53), point a'. If we do this for a lot of distinct parametric samples, we could plot the path of a' and make an animation of it:

Bezier_1_big
Author: Phil Tregoning, Wikimedia Commons

Well, um, OK. This is all just a linear interpolation between a and b. (P0 and P1 in Phil’s animation). But here is the kicker: Who’s to say that a and b aren't motionless? Could it be that a is actually an a' transiting in a lerpish kind of way between another a and b, and that b is also, actually, a b' transiting in a lerpish kind of way between b and c, and that our resultant is really a''? That is, we’re stacking a lerp on a pair of lerps: enacting a 2nd degree pyramid algorithm. Time for another slovenly diagram:
post_02
A stack of lerps

This added layer of lerps finds a' traversing between a and b, a linear path, as well as b' traversing between b and c, another linear path. But we also have a'' lerping betwixt a' and b', both endpoints undergoing linear traversals as a' traverses between them. So what does _a''_s traversal path look like?
Bezier_2_big
Author: Phil Tregoning, Wikimedia Commons

Algebraically collapsing the pyramidical mesh of lerps finds the straight polynomial relation in t and control points a, b, and c, a quadratic:
post_05

We can pile lerps arbitrarily high for any degree Bézier curve, not just linear, quadratic or cubic flavors. Now, since I haven't run out of slovenly diagrams yet, here's the pyramid for a fourth degree quartic Bézier, which takes five control points, a, b, c, d and e and bootstraps it's way up to a'''', the plot of a point on the quartic Bézier curve at some parameter t, drawn from the closed unit interval [0…1]:


Pyramid for a 4th degree Bézier curve

I've labeled this somewhat differently, marking all the places where a G'MIC ‘lerp(a,b,t)’ enters into play with various data items and their primes.

Now, if we look at this diagram sideways and squint a little bit, we can maybe see a way to dress it up G’MICally.

  1. Data structure: an array of arrays. Outer array length (width) equals the number of control points we care to deal with. Say, four for cubic Béziers, five for quartics, three for quadratics. Inner array length (depth) equals the spatial dimensions of the space where we plot points, two for surfaces, three for volumes and so forth. Let me label the elements of the outer array a, b, c, d, e instead of 0, 1, 2, 3 and 4.

  2. Process: For a given parametric value, t, taken from the closed interval [0…1]
    a. Bottom rank:
    a.1 ‘lerp(a,b,t)’, put the results, a' in array position a; we won’t use a again and can overwrite it.
    a.2 ‘lerp(b,c,t)’, put the results, b' in array position b; we won’t use b again and can overwrite it.
    a.3 ‘lerp(c,d,t)’, put the results, c' in array position c; we won’t use c again and can overwrite it.
    a.4 ‘lerp(d,e,t)’, put the results, d' in array position d; we won’t use d again and can overwrite it.
    b. Next rank up and each rank following to the top of the pyramid: Ditto the bottom rank, but, for the second rank, only use array elements a, b, c, and d, which now contain lerps a', b', c', and d'. We ignore array element e. We don’t need it any more.

Overall, the game is: pair-wise lerp array elements from left to right, putting the results in the left hand slot of the current pair, and then bootstrap up to the next rank, where we do it again, removing the rightmost array element from the computation and pair-wise lerping the remainder. At the top of the pyramid, we ‘lerp(a''',b''',t)’ to produce the point a'''' on the quartic curve corresponding to parametric value t. For the visually-oriented among you, all of this is just following the network of computations depicted in my previous slovenly diagram.

Here is an implementation:

bezlength.gmic
bezlength:
   check "isint(${1=10}) && ${1}>0 && isint(${2=512}) && ${2}>128"
   plotcount=$1
   imgsz=$2
   if size([$[]])==0
      # Default control point set; format: # xxxxx…^yyyyy…
      (0.0;-0.75;-0.9;0.1;0.9;0.95;0.6^0.0;-0.675;-0.9;-1.0;-0.9;-0.5;0.4)
   fi
   nm. cpoints
   permute[cpoints] xczy
   (0;1)
   nm. param
   r. 1,$plotcount,1,1,3
   $imgsz,$imgsz,1,1
   nm. canvas
   bzcl={"
           CP=crop(#$cpoints);
           T=crop(#$param);
           const vl=h#$cpoints;
           const ccnt=s#$cpoints;
           const pcnt=h#$param;
           const psz=pcnt*vl;
           const sspsz=pcnt*(vl+1);
           const sscpsz=ccnt*(vl+1);
           const slsz=pcnt-1;
           const rcnt=ccnt-1;
           P=vectorpsz(0);

           # Traverse pyramid; find plot points residing on Bézier
           repeat(
                     size(T),i,
                     LCP=CP;
                     t=T[i];
                     repeat(
                              rcnt,j,
                              repeat(
                                       rcnt-j,k,
                                       PT=lerp(LCP[vl*k,vl,1],LCP[vl*(k+1),vl,1],t);
                                       repeat(vl,m,LCP[vl*k+m]=PT[m])
                                    )
                           );
                     repeat(vl,m,P[vl*i+m]=LCP[m])
                 );

           # 1., 2., not relevant to pyramid alogorithms; these draw diagnostic
           # curves on the canvas and may be omitted, but retain 3.  

           # 1. Homogeneous Screen xform: 2 × Unit square, centered origin → canvas image
           specw=h#$canvas/2;
           id=eye(3);
           id[0]=specw/1.01;
           id[1]=0;
           id[2]=specw;
           id[3]=0;
           id[4]=-specw/1.01;
           id[5]=specw;

           # 2. Draw Bézier curve on canvas 
           SSP=vectorsspsz(0);
           repeat(pcnt,k,SSPPT=id*[P[vl*k,vl,1],1];repeat(vl+1,m,SSP[k*(vl+1)+m]=SSPPT[m]));
           SSCP=vectorsscpsz(0);
           repeat(ccnt,k,SSCPT=id*[CP[vl*k,vl,1],1];repeat(vl+1,m,SSCP[k*(vl+1)+m]=SSCPT[m]));
           repeat(pcnt-1,k,pt0=SSP[k*(vl+1),2,1];pt1=SSP[(k+1)*(vl+1),2,1];polygon(#$canvas,-2,[pt0,pt1],1,0xffffffff,255));
           repeat(ccnt-1,k,pt0=SSCP[k*(vl+1),2,1];pt1=SSCP[(k+1)*(vl+1),2,1];polygon(#$canvas,-2,[pt0,pt1],1,0xffffffff,127));
           repeat(ccnt,k,pt0=SSCP[k*(vl+1),2,1];ellipse(#$canvas,pt0,2,2,0,1,190));

           # 3. Find and export curve length: diff plots; find lengths; sum lengths.
           DIFF=P[vl,(pcnt-1)*vl,1]-P[0,(pcnt-1)*vl,1];
           SLEN=vectorslsz(0);
           fill(SLEN,k,norm2(DIFF[vl*k,vl,1]));
           sum(SLEN)
        "}
   text[canvas] "Est. curve length: "$bzcl", Samples: "$plotcount".",5%,5%,15,1,255
   rm[^-1]

And here’s some examples. The approximate length of a Bézier curve, at various resolutions (first argument - plot count, second argument: plotting image size):

bezlen_5
gmic bezlength.gmic (-1;-0.5;0;0.5;1^0;0.8;0;-0.625;0.7) bezlength. 5,512

bezlen_10

gmic bezlength.gmic (-1;-0.5;0;0.5;1^0;0.8;0;-0.625;0.7) bezlength. 10,512

bezlen_20

gmic bezlength.gmic (-1;-0.5;0;0.5;1^0;0.8;0;-0.625;0.7) bezlength. 20,512

bezlen_500

gmic bezlength.gmic (-1;-0.5;0;0.5;1^0;0.8;0;-0.625;0.7) bezlength. 500,512

Convergence to a useful results does not take too many samples for uncomplicated curves such as this. Convergence slows with more control points. Spline length follows from the approximate lengths of the constituent curve segments, added together.

Hope this helps.

3 Likes