Reptorian G'MIC Filters

[Reptorian]

@garagecoder

Use usual optimisation tricks… make your loops contain as little as possible etc.

That’s exactly what I did.

Here’s the current base for Fragment Blur (The code in rep_frblur uses every of those below, but that is not the base):

#@cli rep_shift_distance_angle_xy: _offset>0,_angle,_boundary_condition={ 0=None | 1=Neumann | 2=Periodic | 3=Mirror },_interpolation={ 0=nearest_neighbor | 1=linear }
#@cli : Offset images at angle and distance.
#@cli : (eq. to 'rep_sdaxy').\n
#@cli : _offset refers to the displacement of image.\n
#@cli : Default value: 'boundary_condition=3','interpolation=1'\n
#@cli : Author: Reptorian.
rep_shift_distance_angle_xy:
skip ${3=3},${4=1}
rad_ang={$2/180*pi}
shift {cos($rad_ang)*$1},{sin($rad_ang)*$1},0,0,$3,$4
#@cli rep_dupsdaxy: eq. to 'rep_duplicate_by_shift_distance_angle_xy'. : (+)
rep_dupsdaxy: rep_duplicate_by_shift_distance_angle_xy $*
#@cli rep_duplicate_by_shift_distance_angle_xy: duplicates_count>=2,_radius_offset>0,_offset_angle,_keep_original={ 0=remove_original | 1=keep_original },_boundary_condition={ 0=None | 1=Neumann | 2=Periodic | 3=Mirror },_interpolation={ 0=nearest_neighbor | 1=linear }
#@cli : Creates copies of image at a distance from original point.
#@cli : (eq. to 'rep_dupsdaxy').\n
#@cli : _radius_offset refers to the distance the duplicates are from the original image. Each circular duplicates are of the same distance. radius_offset can be in percentage form, or in integer form. When it is in percentage form, the length of each duplicates offset is equal to the percentage of diagonal of image.
#@cli : _offset_angle refers to the starting angle the original duplicate is at. The primary duplicate starts from the right.\n
#@cli : Default value: '_radius_offset={sqrt(w^2+h^2)*.05}','_offset_angle=0','_keep_original=0','_boundary_condition=0','_interpolation=1'\n
#@cli : Author: Reptorian.
rep_duplicate_by_shift_distance_angle_xy:
skip ${2={norm(w,h)*.05}},${3=0},${4=0},${5=1},${6=1}
if $1<2 v + error "Invalid duplicate numbers!" fi
aang={360/$1}
repeat $1 +rep_shift_distance_angle_xy[0] $2,{$>*$aang+$3},$5,$6 done if !$4 rm[0] fi

rep_frblur_1:
skip ${2={sqrt(w^2+h^2)*.05}},${3=0},${4=0},${5=1},${6=1}
sh 0,{s-2} sh.. {s} maxa={iM#-1} /. $maxa *.. . rm[-2,-1]
rep_duplicate_by_shift_distance_angle_xy $1,$2,$3,$4,{abs($5)},$6
add
sh 0,{s-2} sh.. {s} f.. i/(i0#-1?i0#-1:1) *. {$maxa/iM#-1} rm[-2,-1]

rep_frblur_2:
skip ${2={sqrt(w^2+h^2)*.05}},${3=0},${4=0},${5=1},${6=1}
rep_duplicate_by_shift_distance_angle_xy $1,$2,$3,$4,{$5!=0?abs($5):1},$6
ti=$!
add / $ti

frblur_1 is for images with alpha, and frblur_2 is for image without alpha.

The code is faster than your j approach, and the result is more in line with the PDN Fragment Blur. The problem becomes apparent when using images with much bigger dimension. Like as in 5000x3000 and then it becomes clear that it’s 10 times slower than the one in PDN.

[garagecoder]

@Reptorian

Like as in 5000x3000 and then it becomes clear that it’s 10 times slower

I think that’s because of memory allocation. I notice rep_duplicate_by_shift_distance_angle_xy makes copies of the image. If you can avoid creating copies, that probably won’t be a problem. I find it strange image would be slower than duplication, it’s one of the most efficient commands in g’mic.

[garagecoder]

An example:

gcd_fragblur : skip ${1=3},${2=0.05}
  100%,100%,100%,100%
  repeat $1
    j. ..,{T=$>*2*pi/$1;[round($2*[cos(T)*w,sin(T)*h]),0,0,1/$1]}
  done

I notice you use neumann boundary, if you must use that then indeed shift could be better - just try not to create more than one image at a time. This could also be done with a convolution kernel, I’ll maybe try that too.

Edit: here it is

gcd_fragblur_convolve : skip ${1=5},${2=1},${3=21}
  $3,$3 eval "C=floor(($3-1)/2);repeat($1,k,T=k*2*pi/$1;P=round(C*[cos(T),sin(T)]+C);i(P[0],P[1])=1)"
  normalize_sum.
  repeat $!-1 l[$>,-1]
    c={ceil(([w,h,d]-1)/2)} whd={0,[w,h,d]-1}
    convolve.. .,1,0,1,$c,0,0,0,$whd,1,1,1,$2,$2,$2
  endl done rm.

gcd_fragblur params: count, distance
gcd_fragblur_convolve params: count, scale, kernel width

[Reptorian]

I changed convolve to convolve_fft and gotten decent times with it. The current gcd_fragblur_convolve wasn’t much faster than original implementation. The first one is the fastest, but way off the ground truth picture.

[garagecoder]

I did some experiments with larger images. It can be surprising which things are faster at larger scales, probably due to multi-threading. Anyway, this was about as fast as I could get using shift:

gcd_fragblur_shift : skip ${1=5},${2=0.05}
  repeat $! l[$>]
  100%,100%,100%,100% 100%,100%,100%,100%
    repeat $1
      j. [0]
      shift. {T=$>*2*pi/$1;round($2*[cos(T)*w,sin(T)*h])},0,0,1
      add.. .
    done
    k.. div $1
  endl done

The basic idea is to avoid memory allocations and just keep copying and shifting within one image.

[Reptorian]

I don’t get how come that one is faster than this one I made which avoids copying new image. But, I think I will be using that code. Just gotta modify it a bit to make it suitable for my filter.

rep_duplicate_by_shift_distance_angle_xy_add:
100%,100%,100%,100%
ang,shift_ang={(2*pi)/$1},{deg2rad($3)}
repeat $1
 +shift[0] {cos($>*$ang+$shift_ang)*$2},{sin($>*$ang+$shift_ang)*$2},0,0,$4,$5
 add[-2,-1]
done

I guess it’s the plus in this modification of my own is why it’s slower.

EDIT:
This is my final code for use. I would use the convolve_fft instead if it allows for float numbers instead of integer as that was 3 s instead of 6 s.

rep_gcd_duplicate_by_shift_distance_angle_xy_add:
skip ${2={norm(w,h)*.05}},${3=0},${4=0},${5=1},${6=1}

100%,100%,100%,100% .

$1,1,1,2,"begin(
  const iter_ang=(2*pi)/$1;
  const shift_ang=deg2rad($3);
 );
 [cos(x*iter_ang+shift_ang)*$2,sin(x*iter_ang+shift_ang)*$2]"

repeat $1
 j.. [0]
 shift.. {I(#-1,$>)},0,0,$5,$6
 add... ..
done

if $4 k[0,-3] add else k... fi

[Reptorian]

In the G’MIC 3.0 thread, it turns out that image/j actually accepts linear opacity. That means it would be possible to create a much faster version of Fragment Blur. There’s a catch though, subpixel-processing wouldn’t be a thing here which is why I have chosen shift. The processing speed should be nearly the one and the same as in PDN.

[Reptorian]

Here’s my three base code to use for Fragment Blur:

The top two is for nearest neighbor approach. The latter is for case of boundary condition.
The bottom one is for bilinear condition regardless of boundary. All of them are faster than the current implementation.

rep_duplicate_by_shift_distance_angle_xy_add_round:
100%,100%,100%,100%

$1,1,1,2,"begin(
  const iter_ang=(2*pi)/$1;
  const shift_ang=deg2rad($3);
 );
 [round(cos(x*iter_ang+shift_ang)*$2),round(sin(x*iter_ang+shift_ang)*$2)];"
 
repeat $1
 j[1] [0],{I(#-1,$>)},0,0,-1
done

if $4 k[0,1] add else k[1] fi
rep_duplicate_by_shift_distance_angle_xy_add_boundary_round:
distoff={ceil($2)}

100%,100%,100%,100% 

+expand_xy[0] $distoff,$5

$1,1,1,2,"begin(
  const distoff=$distoff;
  const iter_ang=(2*pi)/$1;
  const shift_ang=deg2rad($3);
 );
 [round(cos(x*iter_ang+shift_ang)*$2-distoff),round(sin(x*iter_ang+shift_ang)*$2-distoff)];"
 
repeat $1
 j[1] [-2],{I(#-1,$>)},0,0,-1
done

if $4 k[0,1] add else k[1] fi
rep_gcd_duplicate_by_shift_distance_angle_xy_add_bilinear:
skip ${2={norm(w,h)*.05}},${3=0},${4=0},${5=1},${6=1}

100%,100%,100%,100% .

$1,1,1,2,"begin(
  const iter_ang=(2*pi)/$1;
  const shift_ang=deg2rad($3);
 );
 [cos(x*iter_ang+shift_ang)*$2,sin(x*iter_ang+shift_ang)*$2]"

repeat $1
 j.. [0]
 shift.. {I(#-1,$>)},0,0,$5,1
 add... ..
done

if $4 k[0,-3] add else k... fi

[Reptorian]

Optimized Fragment Blur has been released. I don’t know if it works though. At the cli level, I think it does. For GUI, I don’t know as I only have access to 2.9.9 for the GUI version. The changes can be from 120% to 715% faster. Dependent on setup.

Optimization to rep_fragment_blur

[Reptorian]

Earlier post is deleted. Now, I had upgraded Emboss-Relief! I also adjusted Construction Material Texture in light of the update though I don’t know if it’ll work since I don’t have 3.0.0. If you could test it, that’d be wonderful.

@samj Now, it has been upgraded. Consider replacing the one you have with this one.

Here’s the changes:

• Faster Processing. 1D convolution, and the plus prepend operator. Also, several command command tricks as they’re a fancy goto. Yes, I know gotos aren’t recommended.
• Slightly more flexible in output. Not just color is changed.
• Color Output is enabled.

[Reptorian]

I have something interesting here:

$ 512,512 rep_mt_test
rep_mt_test:
seed={u(0,1)}

{w*2},{h*2},$_cpus,4

half_diag={norm(w,h)/2}

pal 150

$_cpus,1,1,1,:"begin(
  srand($seed);
  diagonal(a,b)=sqrt(a^2+b^2);
  const ss=s#-2;
  const nc=w#-1;
  const nt=$_cpus;
  const tau=2*pi;
  const nt_mi=nt-1;
  const ww=w#-2;
  const hh=h#-2;
  const mx=ww-1;
  const my=hh-1;
  const max_dist=diagonal(mx,my);
  const min_dist=max_dist*.1;
  const steps=max(2,50);
  const div_step=steps-1;
  const radius=15;
  const length=sqr(radius)*ss;
  const off=ceil(radius/2);
  md_arr=expr('x',nt); # Matrix to define where to draw circle / Matriz para dibujar circulos
  fy_shuffle()=(    # Fisher-Yates Shuffle md_arr
   for(p=nt_mi,p>0,p--,
    r_ind=int(u(0,nt));
    swap(md_arr[p],md_arr[r_ind]);
   );
  );
  opacity=expr('
   begin(
    const half_dist=(w-1)/2;
    remap(a)=(a-half_dist)/half_dist;
   );
   xx=remap(x);
   yy=remap(y);
   norm(xx,yy)<=1;'
  ,radius,radius,1,1);
  paint_circles()=(
   paint_color=[I(#-1,int(u(nc)),0,0),255];
   sqr_col=resize(paint_color,length,1,0);
   ang=u(tau);
   dist=u(min_dist,max_dist);
   step_dist=dist/div_step;
   shift_x=cos(ang)*step_dist;
   shift_y=sin(ang)*step_dist;
   pos_x=int(u(ww));
   pos_y=int(u(hh));
   repeat(steps,
    draw(#-2,sqr_col,pos_x-off,pos_y-off,md_arr[x],0,radius,radius,1,ss,1,opacity,1);
    pos_x+=shift_x;
    pos_y+=shift_y;
   );
  );
  fy_shuffle();
 );
 repeat(20,
  critical(fy_shuffle(););
  paint_circles();
 );"

l... s z blend alpha endl

The big block of code inside rep_mt_test draws lines in parallel at random z using an array with unique integer number. Is this thread-safe? The reason I asked is because I want to draw lines of circles at different zs, and no threads should be drawing on the same z. Print only outputs md_arr once, so I can’t tell what it is on different threads.

To explain some more. fy_shuffle() shuffles an array which is used to define where to draw the lines of circle. And paint_circles() draws lines of circles.

And for those who’s wondering what is this test code for, I’m basing it off the principles found in this plugin for PDN - TR's iPollock (March 3rd 2015) - Plugins - Publishing ONLY! - paint.net Forum

Output:

[Reptorian]

I added this to the end of the code.

end_t(
  print(md_arr);
 );

Results shows that it isn’t thread-safe. I wanted to do multi-threaded drawing always at different z at different threads.

md_arr is supposed to be the same in all threads.

---

EDIT: Never mind, now I got it.

New Test Code

rep_mt_test:
nt,num_rep,pal_id,min_perc,steps,radius={$_cpus},90,150,{cut(.1,0,1)},{max(2,50)},15

$nt,1,1,1,x nm. pz
$nt,1,1,1 nm. swapper

pal $pal_id

m "rpz : eval. >r_ind=int(u(0,"$_cpus"));swap(i(#-2,x),i(#-2,r_ind));"
m "flatten : l. s z blend alpha endl r. {w#0},{h#0},100%,100%,3 sh. 0,2 sh.. 3 /. 255 /.. . *. 255 rm[-2,-1]"

repeat $!-3 l[$>,-3,-2,-1]

 {w#0*2},{h#0*2},$nt,4
 
 repeat $num_rep
  rpz[pz,swapper]
  
  eval[pz] :"begin(
    diagonal(a,b)=sqrt(a^2+b^2);
    const ss=s#-1;
    const nc=w#-2;
    const nt=$nt;
    const tau=2*pi;
    const nt_mi=nt-1;
    const ww=w#-1;
    const hh=h#-1;
    const mx=ww-1;
    const my=hh-1;
    const max_dist=diagonal(mx,my);
    const min_dist=max_dist*$min_perc;
    const steps=$steps;
    const div_step=steps-1;
    const radius=$radius;
    const length=sqr(radius)*ss;
    const off=ceil(radius/2);
    opacity=expr('
     begin(
       const half_dist=(w-1)/2;
       remap(a)=(a-half_dist)/half_dist;
      );
      xx=remap(x);
      yy=remap(y);
      norm(xx,yy)<=1;'
      ,radius,radius,1,1);
    paint_circles()=(
     paint_color=[I(#-2,int(u(nc)),0,0),255];
     sqr_col=resize(paint_color,length,1,0);
     ang=u(tau);
     dist=u(min_dist,max_dist);
     step_dist=dist/div_step;
     shift_x=cos(ang)*step_dist;
     shift_y=sin(ang)*step_dist;
     pos_x=int(u(ww));
     pos_y=int(u(hh));
     repeat(steps,
      draw(#-1,sqr_col,pos_x-off,pos_y-off,i,0,radius,radius,1,ss,1,opacity,1);
      pos_x+=shift_x;
      pos_y+=shift_y;
     );
    );
   );
   paint_circles();"
 done
 flatten
endl done

um rpz

[Reptorian]

Looks like another pause to another script of mine. This one, well, I wish I knew how to calculate length of cubic spline. There is 5 point langrage formula, but I don’t know how to use it.

Here’s another:

Based on this idea: Hitomezashi Stitch Patterns - Numberphile - YouTube

And here’s the original spreadsheet code - Mobile Numbers: Hitomezashi Stitching | The Aperiodical (Found in bottom of blog post)

+rep_binary_hitomezashi:
skip "${3=}" #If there is no 3rd argument, then $3 is treated as null. So therefore, 2 arguments.

if narg($3) pass$3 1 # If there is a 3rd argument
else
 # Vowels as 1 is part of the rule. 
 # Whether number is even is also part of rule. 
 # Emoticons as 1 is part of rule. 
 # There's some few more arbitrary rule.
 (0,1,1,1,1,1,0,0,0,0,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,0,0,1,1,0,1,1,0,1,1,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,1,1,1,1,1,1,1,0,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,1,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,0,1,0,0,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,0)
fi

('"$1"') # Get array of unicode value of input 1
('"$2"') # Get array of unicode value of input 2

f[-2,-1] i(#-3,i) # Convert unicode image/array to binary with rules.

{w#-1},1,1,1,>begin(v=0;);v=x?(j(#-1,-1)+v)%2:0; rm.. # Process first column. Creates auxilatory surface and then remove reference surface.

{w#-2},{w#-1},1,1 # Create a 2D surface to write on

eval.. "
 v=i(#-1,0,x)=i;
 row=x+1;
 repeat(w#-3-1,p,
  v=i(#-1,p+1,x)=sum(v,i(#-3,p,0),row,1)%2;
 );" #Generate value from spreadsheet rule.

rm[-4--2] # Keep only the Hitomezashi Pattern

---

[grosgood (Garry R. Osgood)]

@Reptorian

I wish I knew how to calculate length of cubic spline.

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.

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:

Author: Phil Tregoning, Wikimedia Commons

Well, um, OK. This is all just a linear interpolation between a and b. (P₀ and P₁ 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:

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?

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:

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 4^th 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):

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

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

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

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.

[Reptorian]

I’ll keep that in note.

More work on the Hitomezashi project:

That’s a beauty.

Leaving this as backup for fix later

gmic bin=111000011111100011000001110100110000001110011000010010001111101001110000011100010010011011100010 +rep_binary_hitomezashi $bin,$bin pal steamlords,wcmyk tic +_rep_colorize_binary_hitomezashi[0] [-2],[-1] toc r. 200%,200%,100%,100%,0,3,0,0

[afre (Alan)]

Unfortunately, the pattern is repetitive and less distinct about 75% to the edges from the centre.

[Reptorian]

You can have distinct output:

But, yeah you’re right that it looks more uniform when the number of colors are less than the number of shapes with areas greater than 1.

Here’s something else:

[Reptorian]

Just dropping this here. I made it animate-able. Only thing left to do is to allow users to have a upscaled version, and whether borders are there.

[Reptorian]

Is there already a command that does this? I’m aware of the resize and ax/ay/az/ac options , but that creates 2 pixel border on the center. This command I created solves it and these are the four options. I’ll push it though I"ll remove it if this is already implemented. My hitomezashi will be using this command.

---

With the new rep_symmetrize_xy, I used it on the hitomezashi.

[afre (Alan)]

Where is the 2 px border? Do you mean between the square patterns in the second image? That you could get rid of if you do some preprocessing before appending.