Reptorian G'MIC Filters

You can avoid looping over every value and use built-in image looping by replacing the end like this:

round.. 1
$1,1,1,{s#-1} $1,1,1,1
f[0] "i[#-1,i]+=1;I[#-2,i]+=I(#1,x,y,z);i"
+eq. 0 add[-2,-1] div[-2,-1] # avoid div by zero

The method there is to build two images during the loop rather than one - the second is a histogram to divide by. That also replaces the fill to round values with native round.

For your second question, we probably need to see an example of what might be in $cropinfo or how it is set.

Edit: just to demo how histogram can do the same thing

round.. 1
$1,1,1,{s#-1}
f[0] "I[#-1,i]+=I(#1,x,y,z);i"
+histogram[0] $1,0,255
+eq. 0 add[-2,-1] div[-2,-1]

Amazingly that works so well.

C:\Windows\System32>gmic sp cat +to_gray rv tic new_rep_regm 256 toc
[gmic]-0./ Start G'MIC interpreter.
[gmic]-1./ Input sample image 'cat' (1 image 600x550x1x3).
[gmic]-1./ Force image [0] to be in GRAY mode.
[gmic]-2./ Reverse positions of images [0,1].
[gmic]-2./ Initialize timer.
[gmic]-3./ Elapsed time: 0.022 s.

I think I will work from it as a base and see how I can fill in the missing value as well as looping along multiple image and convert them to gradient.

The second question, I learned that it doesn’t support == or !=, but it will support >, and <. I wonder if @David_Tschumperle can address that limitation of the math parser in gmic.

EDIT: I never realized that fill can access outside images and fill them in. Interesting. I must experiment with that actually.

EDIT:

new_rep_regm:
round.. 1
$1,1,1,{s#-1} $1,1,1,1
f[0] "i[#-1,i]+=1;I[#-2,i]+=I(#1,x,y,z);i"
+f. i!=0
+eq. 0 add[-3,-2] div[-3,-2] # avoid div by zero
inpaint_pde.. .
k..

This works to find gradient map of image.

1 Like

I don’t get it. What limitation ?

Here.

This works.

rep_test:
5,1,1,1,"begin(st=[0,1,2,3,4]);
m=st>1;
m[x];
"

Error!

rep_test:
5,1,1,1,"begin(st=[0,1,2,3,4]);
m=st==1;
m[x];
"

That’s because operators > and < are useless in a vector space, because a vector space is not ordered. So by default, these operators are returning also a boolean vector, when used with two vector-valued operands.

Operators == and != are (extremely) useful to test vector equality, so they always return a boolean value.

If you want to return a vector that test equality of each component, you can write this:

veq(A,B) = !(A-B);

C = veq(A,B);  // == between vector components
D = !veq(A,B); // != between vector components
1 Like

Thanks to @garagecoder’s help, I had upgraded Reverse Engineer Gradient Map/rep_regm filter. It is now completely usable in cli, and gui as well.


I added multi-threading support for Chirikov-Taylor GUI/CLI filter. 6 times faster than serial processing in my computer. Also, added Alpha Support for Grayscale Mode in Chirikov-Taylor GUI filter.


Rewritten the rep_sptbwgp again. It’s faster. However, I just don’t get why there’s a speed difference in horizontal and vertical processing regardless of using a rotated image or not. Nothing makes sense about it. Permute doesn’t really solve it seems.


Now, I had managed to optimize rep_axis_streak_color again. I think I am done for the new year upgrade to my filters. I still however would like to be able to solve the issue of why vertical processing is slower than horizontal processing in rep_sptbwgp, and never did, and this is the only case where this bug shows up in terms of gmic processing or I am missing something. It doesn’t show up in other gmic scripts that utilize processing along row/columns, just this one which points to “missing something”, and nothing seem to indicate what it is.


Roadmap for 2021:

  1. Refactor Naeddyar mapmaking script to more suitable version. He wrote scripts with little knowledge of gmic scripting. So, hence the need of refactoring. ← I’m here now.

  2. Refactor few of my script and adding new feature.

  3. New filters only after the refactor. I won’t be refactoring my script for optimization, etc.


Looks like I may work on Perspective Streak, and addressing problems. However, there is a issue I ran into with writing the Axis Streak Part. The axis streak editing part is meant to address a shortcoming with not being able to preserve image structure.

The problem: When setting _add2new to 1, the output is not as expected because if I were to blend the output with the old image, then areas where there are already non-zero transparent pixel would be off.

How do I address it with this code?

EDIT: I had addressed it. Realized I can just insert temp vector to new image before calculating newcol and newalp.

New Axis Streak
#@cli rep_axis_streak_color: orientation,direction,_alpha_exponential_factor,_maxval>0,_add2new={ 0=add2old | 1=add2new },_cmykmode={ 0=non-cmyk | 1=cmyka_mode }
#@cli : Streaks colored pixels taking into account of opacity.
#@cli : '_alpha_exponential_factor' is used to manipulate the alpha mixing within pixels. The more power that is assigned to the alpha, the more mixing there would be.
#@cli : '_maxval' divides the alpha channel internally to normalize ranges to 0-1. A error will appear if not normalized. By default, it used the max alpha channel value.
#@cli : '_add2new' defines whether to streak data on old image or to streak pixel data onto a new image. Used for Perspective Streak.
#@cli : '_cmykmode' is only used in case of using only cmyk images. Not needed in normal cases at all.
#@cli : Default values: '_alpha_exponential_factor=0','_maxalp=n/a','_cmykmode=0'
#@cli :
#@cli : Author: Reptorian.
skip ${3=0},${4=},${5=0},${6=0}

tcr=3

if $6 tcr+=1 fi

if narg($4) if $4==0 return fi fi

repeat $! l[$>]

 if s==1||s==$tcr break fi

 vv=0
 repeat s
  sh $>
  vv+={iv#-1}
  rm.
 done
 if !$vv break fi

 sh. 0,{s-2}
 sh.. {s}
 f.. i#-1?I
 if narg($4) alp={abs($4)}
 else alp={iM#-1}
 fi
 /. $alp
 if iM#-1>1||im#-1<0 error alpval(valid)==F fi
 if $3<=-1 exp_f=-{1-10^-8} else exp_f=$3 fi
 f. i^(1+$exp_f)
 
 if $5 {w#0},{h#0},{d#0},{s#0}
  col_info=-3
  alp_info=-2
  targ_info=-1
 else
  col_info=-2
  alp_info=-1
  targ_info=0
 fi

 if $1
  outdata_dim={w},1,{d},{s#0}
  outdata_coords=x,yy
  if $2 direction=yy=hh-1,yy>=0,yy--
  else direction=yy=0,yy<hh,yy++
  fi
 else
  outdata_dim=1,{h},{d},{s#0}
  outdata_coords=xx,y
  if $2 direction=xx=0,xx<ww,xx++
  else direction=xx=ww-1,xx>=0,xx--
  fi
 fi

 $outdata_dim,":begin(
  const ww=w#0;
  const hh=h#0;
 );
 start_val=1;
 for("$direction",
  start_val?(
   start_val=0;
   col=I(#"$col_info","$outdata_coords",z);
   alp=i(#"$alp_info","$outdata_coords",z,0);
   temp=[col,alp];
  ):(
   newcol=I(#"$col_info","$outdata_coords",z);
   newalp=i(#"$alp_info","$outdata_coords",z,0);
   newinfo=[newcol,newalp];
   !newalp?(I(#"$targ_info","$outdata_coords",z)=temp;):
   newalp==1?(temp=newinfo;
   ):(
      col=newcol*newalp+(1-newalp)*col;
      alp=alp+newalp*(1-alp);
      temp=[col,alp];
      I(#"$targ_info","$outdata_coords",z)=temp;
   );
  );
 );"
 if $5 rm[^-2]
 else rm[1-3]
 fi
 sh. {s-1} *. $alp rm.
endl done

Looks like upgrading Perspective Streak is going to take a while. That being said, I added 3 new formulas to Thorn Fractal:

  • Chaotic Hooks Unearthing
  • Sinusoidal Liquid
  • Cosinusoidal Liquid

Should be here after 1 hour.


Another 3 Formulas has been added to Thorn Fractal.

Side note: I realized with what @garagecoder shown me that filling outside images with fill is possible, it is actually feasible to do multithreaded version of Diffusion Limited Aggregation.

EDIT: I tried, didn’t really saved time. That means I’m done with current filters except for Perspective Streak.


1/11/2021:
Okay, now I have a new update. It seems that I am a bit closer to upgrading Perspective Streak. I’m working on transformation now. Now recpoltrans -3 to recpoltrans -2 takes .32 s as opposed to 1.5 s. More than 4.6 times faster. It shows promise at preserving more details too. After transformation, I would need to add more options to axis_streak_color just for the fact that it is crucial to make it work correctly with it.


1/14/2021:

Wow! Now I think I solved recpoltrans! The new changes is by far more faster. From 1.5 s to .194 s. It is also more accurate. I only worked on if -3, and 2 part of rep_recpoltrans since they matter a lot more than the others as they’re both designed for precision. The others if in rep_recpoltrans don’t really matter too much for me to fix them. In addition to these changes, I’m closer to fixing a bug in perspective streak!

See here.

Left: Original Middle: recpoltrans -3 to recpoltrans 2 Right: Xor Analysis

[gmic]-2./ Elapsed time: 0.19 s.

Also, saving code here:

Code of new_reprecpoltrans
new_rep_recpoltrans:
skip ${1=0},${2=0},${3=0}
if $1<-1||$1>1 error "($1>=-1&&$1<=1)=0" fi
if $2<-1||$2>1 error "($2>=-1&&$2<=1)=0" fi
repeat $! l[$>]
 ov=${-rep_2dcr}
 if $3==-3
  maxlength={max(w,h)*1.5}
  perimeter={(w+h)*2}
  {$perimeter},{$maxlength},{d},{s},":begin(
    surface=expr('begin(const dpi=2*pi;);(x/w)*dpi;',w);
    const dpi=pi*2;
    const ww=w#0-1;
    const hh=h#0-1;
    const ihh=h-1;
    const point_x=(($1*-1)*.5+.5)*ww;
    const point_y=($2*.5+.5)*hh;
    const inv_point_x=ww-point_x;
    const inv_point_y=hh-point_y;
    const cut_ang_s0=abs(atan2(inv_point_y,inv_point_x));
    const cut_ang_s1=pi-abs(atan2(inv_point_y,point_x));
    const cut_ang_s2=pi+abs(atan2(point_y,point_x));
    const cut_ang_s3=dpi-abs(atan2(point_y,inv_point_x));
    distanceaway=[ww-point_x,hh-point_y,point_x,point_y];
   );
   surface_angle=surface[x];
   surface_angle>cut_ang_s0&&surface_angle<=cut_ang_s1?side=1:
   surface_angle>cut_ang_s1&&surface_angle<=cut_ang_s2?side=2:
   surface_angle>cut_ang_s2&&surface_angle<=cut_ang_s3?side=3:
   side=0;
   mdist=abs(side%2?1/sin(surface_angle):1/cos(surface_angle));
   dix=(point_x+cos(surface_angle)*distanceaway[side]*mdist*y/ihh);
   diy=(point_y+sin(surface_angle)*distanceaway[side]*mdist*y/ihh);
   I(#0,round(ww-dix),round(hh-diy),z);"
   
  k.
 elif $3==-1||$3==-2
  f ":begin(
   const point_x=(($1*-1)*.5+.5)*w;
   const point_y=($2*.5+.5)*h;
   const inv_point_x=w-point_x;
   const inv_point_y=h-point_y;
   const cut_ang_s0=abs(atan2(inv_point_y,inv_point_x)*180/pi);
   const cut_ang_s1=180-abs(atan2(inv_point_y,point_x)*180/pi);
   const cut_ang_s2=180+abs(atan2(point_y,point_x)*180/pi);
   const cut_ang_s3=360-abs(atan2(point_y,inv_point_x)*180/pi);
   distanceaway(value)=(
    value==0?ww-point_x:
    value==1?hh-point_y:
    value==2?point_x:
    point_y;
    );
  );
  surface_angle=(x/w)*360;
  surface_angle>cut_ang_s0&&surface_angle<=cut_ang_s1?side=1:
  surface_angle>cut_ang_s1&&surface_angle<=cut_ang_s2?side=2:
  surface_angle>cut_ang_s2&&surface_angle<=cut_ang_s3?side=3:
  side=0;
  mdist=abs(side%2?1/sin(surface_angle):1/cos(surface_angle));
  dix=(point_x+cos(surface_angle)*distanceaway(side)*mdist*y/h)*((w-1)/w);
  diy=(point_y+sin(surface_angle)*distanceaway(side)*mdist*y/h)*((h-1)/h);
  I(w-(dix+1),h-(diy+1),z,2);
  "
  if $3<-1 r2dx 50%,6 fi
 elif $3==0||$3==1
  if $3>0 r2dx 200%,6 fi
  f ":begin(
   const ww=w-1;
   const hh=h-1;
   const sd=max(w,h)/min(w,h);
   const sx=w>h?sd:1;
   const sy=w>h?1:sd;
   const cx=.5+$1*.5;
   const cy=.5+$2*.5;
   const px=cx*w;
   const py=(1-cy)*h;
   const sxl=(w/2)/px;
   const sxr=(w/2)/(w-px);
   const syt=(h/2)/py;
   const syb=(h/2)/(h-py);
  );
  atx=(x/ww-cx)*sx;
  aty=(y/hh-(1-cy))*sy;
  sur_atan=(atan2(aty,atx)+pi)/(2*pi);
  xl=-1+(x/ww)*2*sxl;
  xr=1-(1-x/ww)*2*sxr;
  yt=-1+(y/hh)*2*syt;
  yb=1-(1-y/hh)*2*syb;
  xx=x>=px?xr:xl;
  yy=y>=py?yb:yt;
  sur_max=max(abs(xx),abs(yy));
  I(sur_atan*w,sur_max*h,z,2,2);
  "
  if $3>0 sharpen 2 fi
 elif $3==2
  $=val
  orientation=${val{$>+4}}
  perimeter={w}
  length_1={$perimeter-h*(4/3)}
  length_2={$perimeter-$length_1}
  if $orientation
   width={min($length_1,$length_2)}
   height={max($length_1,$length_2)}
  else
   width={max($length_1,$length_2)}
   height={min($length_1,$length_2)}
  fi
  v + echo {w},{h} v -
  {$width/2},{$height/2},{d},{s},":begin(
   const ww=w#0-1;
   const hh=h#0-1;
   const nw=w-1;
   const nh=h-1;
   const ox="$width";
   const oy="$height";
   const sd=max(ox,oy)/min(ox,oy);
   if(w>h,
    const sxf=ox>oy?sd:1;
    const syf=ox>oy?1:sd;
   ,
    const sxf=ox<oy?1:sd;
    const syf=ox<oy?sd:1;
   );
   const cx=.5+$1*.5;
   const cy=.5+$2*.5;
   const px=cx*(ww-1);
   const py=(1-cy)*(hh-1);
   const sxl=(ww/2)/px;
   const sxr=(ww/2)/(ww-px);
   const syt=(hh/2)/py;
   const syb=(hh/2)/(hh-py);
  );
  xx=(x/nw)*(ww);
  yy=(y/nh)*(hh);
  xl=-1+(xx/ww)*2*sxl;
  xr=1-(1-xx/ww)*2*sxr;
  yt=-1+(yy/hh)*2*syt;
  yb=1-(1-yy/hh)*2*syb;
  nxx=xx>px?xr:xl;
  nyy=yy>py?yb:yt;
  ay=max(abs(nxx),abs(nyy));
  ax=(atan2((yy/hh-(1-cy))*syf,(xx/ww-cx)*sxf)+pi)/(2*pi);
  I(#0,abs(ax*ww),abs(ay*hh),z,2,1);
  "
  k.
 else error "$3|"$"3!=intnum[-3,2]"
 fi
endl done

EDIT: Okay, there seem to be difference now that I analyzed a different way. Now, I’m having the micro-shift problem caused by atan2. I don’t know how to fix that though. :confused:

Here’s an analysis of difference. I decided to try to map coordinate into image, and then find the difference between x,y. This is the result. The red is x, and the y is green. Max difference is 3 pixels off. That’s not good.


1/17/2021:

Looks like I have a much better version of rep_recpoltran now. My next step is to make changes to axis_streak to address the shortcomings with the new rep_recpoltrans, and finally a even better version of Perspective Streak. One that is not that slow and preserve details better.

1 Like

Does anyone find this useful? This one works for volumetric image or videos. The algorithm could be better though, I’ll admit.

rep_mode_find:
repeat $! l[$>]
 {w},{h},1,{s},":begin(
   const dd=d#-1;
  );
  v_m=vectordd(0);
  for(m=0,m<dd,m++,
  ref_pixel=I(#-1,x,y,m);
  for(n=m,n<dd,n++,
   if(I(#-1,x,y,n)==ref_pixel,
    v_m[m]++;
   );
  );
 );
 pos=max(v_m)-1;
 pos?I(#-1,x,y,pos):vectors(nan);
 "
 rm..
endl done

This is able to extract features which are split over the Z axis…?

It extract the mode over the z-axis. It does not detect ambiguous cases like {3,3,3,4,1,5,5,5}. So yes.

I’m thinking of a certain algorithm for this:

  1. Create two vectors, one for the terms found and one for the frequencies.
  2. Find the position(s) of the maximum value(s) in that frequency vector.
  3. Look up the numbers from the vector containing the terms and store them in a vector which gets sent to an image whose depth is the maximum number of modes. You may need run() for that.

There are some algorithm that does 1, and 2. I just don’t know how to translate those algorithm I found.

EDIT: Outdated post is outdated.

New Mode CLI Filter (WIP)
  • Case IMGS>2 : Find the mode of all images within the same pixel coordinates. Outputs nan if multiple modes are found or no modes are found.

  • Case IMGS==2: Checks if pixels are the same. If they are, output the current pixel, otherwise nan.

  • Case IMGS==1: Find all the mode within a image. If there are no mode, do nothing. If there are multiple modes, then you can choose which modes to choose or output a nan if no input is detected

$1 == Chosen id of multiple modes. If less than 0, then if there is more than 1 mode, does nothing. Otherwise, fill the image with one of those mode. This only works with case IMGS==1.
rep_mode:
skip ${1=-1}
if $!>2
 +a z colormap. 0
 a[^-1] c
 {w#0},{h#0},{d#0},{s#1},"begin(
  const vs=w#1;
  const ss=s#1;
  const css=s#0/ss;
  endnan=vectors(nan);
 );
 nmaxfreq=0;
 vfreq=vectorvs(-1);
 vp=I#0;
 for(pos=0,pos<css,pos++,
  curcol=vp[pos*ss,ss];
  for(posfreq=0,posfreq<vs,posfreq++,
   coltestpos=I(#1,posfreq,0,0);
   if(curcol==coltestpos,vfreq[posfreq]++;);
  );
 );
 maxfreqcount=max(vfreq);
 if(maxfreqcount+1,
  posmaxfreqcol=find(vfreq,maxfreqcount,0,1);
  for(fs=0,fs<vs,fs++,if(vfreq[fs]==maxfreqcount,nmaxfreq++;););
  if(nmaxfreq==1,I(#1,posmaxfreqcol,0,0);,endnan);
 ,endnan;
 );"
 k.
elif $!==2
 f. "begin(endnan=vectors(nan););I==I#0?I:endnan;" rm..
else
 minv={im}
 - $minv
 
 repeat s
  sh. $>
  bsize$>={iM+1}
  rm.
 done
 
 if s==1
  $bsize0,1,1,1 eval.. ++i(#-1,i,0,0,0)
 elif s==2
  $bsize0,$bsize1,1,1 eval.. ++i(#-1,I,0)
  1,1,1,2
  eval. "
  "
 elif s==3
  $bsize0,$bsize1,$bsize2,1 eval.. ++i(#-1,I)
  
  1,1,1,3 1,1,1,1
  
  eval... :${-math_lib}"if(i,
    dar_insert(#-2,[x,y,z]);
    dar_insert(#-1,i);
   );
   end(
    resize(#-2,1,dar_size(#-2),1,3);
    resize(#-1,1,dar_size(#-1),1,1);
   );"
   
   +.. $minv
   
   
   pixelsort.. -,y,[-1]
   sort. -
   
   u {"freq=0;
    for(n=0,n<h#-1,n++,
     if(i(#-1,0,n)==iM,freq++,break(););
     );
     (freq)"}
     
    maxfreq={${}}
    
    if $1==-1
     if $maxfreq==1
      f[-4] I(#-2,0,0,0)
     fi
     k[-4]
    else
     pos={$1%$maxfreq}
     f[-4] I(#-2,0,$pos,0)
     k[-4]
    fi
 else
  $bsize0,$bsize1,$bsize2,$bsize3 eval.. ++i(#-1,I)
 fi
fi

Hi @Reptorian,

Something which may interest you… thanks to lots of inspiration and guidance from @afre , I’ve made a separable “summed area table” based box filter which allows per-pixel kernel sizes (could be extended to allow arbitrary rectangular kernels per-pixel).

It performs quite fast. I haven’t added this to my own gmic filters yet, but probably will soon - I suggest keeping your own version if you want to use it though:

gcd_satboxfilter_local : check ${"is_image_arg $1"}
  repeat $! pass$1 0 max. 1 l[$>,-1]
    f.. ">i+j(-1)"
    f.. "
    begin(
      const boundary=1; const interpolation=1; const W=w-1;
      L(X)=X>1?i(X-1):X*i(0);
      R(X)=X<W?i(X):(1+X-W)*i(W)-(X-W)*i(W-1);
    );
    S=i#1; K=(S-1)/2;
    (R(x+K) - L(x-K))/S"
    f.. ">i+j(0,-1)"
    f.. "
    begin(
      const boundary=1; const interpolation=1; const H=h-1;
      T(Y)=Y>1?i(x,Y-1):Y*i(x,0);
      B(Y)=Y<H?i(x,Y):(1+Y-H)*i(x,H)-(Y-H)*i(x,H-1);
    );
    S=i#1; K=(S-1)/2;
    (B(y+K) - T(y-K))/S"
  endl rm. done

Example use:

gmic sp tiger sp inside ri. ..,3 n. 0,101 gcd_satboxfilter_local.. .

(Edited to handle a small boundaries glitch)

1 Like

I guess our conversations yielded real life results. :face_with_hand_over_mouth: This is different from what I envisioned. Maybe I am thinking of

1 Like

Yes indeed! I honestly would not even have considered this being possible were it not for your ideas/questions :slight_smile:
There are plenty interesting applications, so I’m very happy to have it.

I think the utility of my involvement in things is that I always give people new ideas to think about.

That forward and backward processing reminds me of some other filtering. I forget what type. You might be on to more than just box filtering.

The part I’m excited about is a nice area blur. The g’mic tutorials mention the “cheesy” result you get when simulating it with a mask. But a couple iters of this and it’s all smooth:

1 Like

Has a fading into unconsciousness feel to it. :male_detective:

@garagecoder Noted. I do have plans to add these filters, but they’re still back in my head for months now.

That being said.

New progress on mode command. I’ll simplify the code soon, it shouldn’t take too long to finish this mode filter.

Reptorian Mode CLI Filter
$1 == Mode to use if there are 1 or more modes. If it less than 0, output will be generated automatically from one mode, and if there are multiple/none, output is original image.
Note that the above is not applicable to 2 or more images. In the case of multiple images, the mode of all images by pixels by pixels is generated instead, and if no modes are found within that pixel coordinate, then output nan within the pixel coordinate.
rep_mode_color:
skip ${1=-1}
if $!>2
 +a z colormap. 0
 a[^-1] c
 {w#0},{h#0},{d#0},{s#1},"begin(
  const vs=w#1;
  const ss=s#1;
  const css=s#0/ss;
  endnan=vectors(nan);
 );
 nmaxfreq=0;
 vfreq=vectorvs(-1);
 vp=I#0;
 for(pos=0,pos<css,pos++,
  curcol=vp[pos*ss,ss];
  for(posfreq=0,posfreq<vs,posfreq++,
   coltestpos=I(#1,posfreq,0,0);
   if(curcol==coltestpos,vfreq[posfreq]++;);
  );
 );
 maxfreqcount=max(vfreq);
 if(maxfreqcount+1,
  posmaxfreqcol=find(vfreq,maxfreqcount,0,1);
  for(fs=0,fs<vs,fs++,if(vfreq[fs]==maxfreqcount,nmaxfreq++;););
  if(nmaxfreq==1,I(#1,posmaxfreqcol,0,0);,endnan);
 ,endnan;
 );"
 k.
elif $!==2
 f. "begin(endnan=vectors(nan););I==I#0?I:endnan;" rm..
else
 minv={im}
 - $minv
 
 repeat s
  sh. $>
  bsize$>={iM+1}
  rm.
 done
 
 if s==1
 
  $bsize0,1,1,1 eval.. ++i(#-1,i,0,0,0) {w},1,1,1,x 
  
  u {"freq=0;
   for(n=0,n<w#-2,n++,
    if(i(#-2,n)==iM#-2,freq++;);
    );
   (freq)"}
  
  maxfreq={${}}
  
  echo $maxfreq
  
  +.. $minv
  
  pixelsort. -,x,[-2]
    
  if $1<0
   if $maxfreq==1
    f... i(#-1,0,0,0)
   fi
   k...
  else
   pos={$1%$maxfreq}
   f... i(#-1,$pos,0,0)
   k...
  fi
  
 elif s==2
 
  $bsize0,$bsize1,1,1 eval.. ++i(#-1,I,0)
  
  1,1,1,2 1,1,1,1
  
  eval... :${-math_lib}"if(i,
    dar_insert(#-2,[x,y]);
    dar_insert(#-1,i);
   );
   end(
    resize(#-2,1,dar_size(#-2),1,2);
    resize(#-1,1,dar_size(#-1),1,1);
   );
  "
  
  +.. $minv
  
  pixelsort.. -,y,[-1]
  sort. -

  u {"freq=0;
   for(n=0,n<h#-1,n++,
    if(i(#-1,0,n)==iM,freq++,break(););
    );
    (freq)"}
     
   maxfreq={${}}
    
   if $1<0
    if $maxfreq==1
     f[-4] I(#-2,0,0,0)
    fi
    k[-4]
   else
    pos={$1%$maxfreq}
    f[-4] I(#-2,0,$pos,0)
    k[-4]
   fi
    
 elif s==3
  $bsize0,$bsize1,$bsize2,1 eval.. ++i(#-1,I)
  
  1,1,1,3 1,1,1,1
  
  eval... :${-math_lib}"if(i,
    dar_insert(#-2,[x,y,z]);
    dar_insert(#-1,i);
   );
   end(
    resize(#-2,1,dar_size(#-2),1,3);
    resize(#-1,1,dar_size(#-1),1,1);
   );"
   
   +.. $minv
   
   
   pixelsort.. -,y,[-1]
   sort. -
   
   u {"freq=0;
    for(n=0,n<h#-1,n++,
     if(i(#-1,0,n)==iM,freq++,break(););
     );
     (freq)"}
     
    maxfreq={${}}
    
    if $1<0
     if $maxfreq==1
      f[-4] I(#-2,0,0,0)
     fi
     k[-4]
    else
     pos={$1%$maxfreq}
     f[-4] I(#-2,0,$pos,0)
     k[-4]
    fi
 else
  $bsize0,$bsize1,$bsize2,$bsize3 eval.. ++i(#-1,I)
  
   1,1,1,4 1,1,1,1
  
  eval... :${-math_lib}"if(i,
    dar_insert(#-2,[x,y,z,c]);
    dar_insert(#-1,i);
   );
   end(
    resize(#-2,1,dar_size(#-2),1,4);
    resize(#-1,1,dar_size(#-1),1,1);
   );"

   +.. $minv
   
   
   pixelsort.. -,y,[-1]
   sort. -
   
   u {"freq=0;
    for(n=0,n<h#-1,n++,
     if(i(#-1,0,n)==iM,freq++,break(););
     );
     (freq)"}
     
    maxfreq={${}}
    
    if $1<0
     if $maxfreq==1
      f[-4] I(#-2,0,0,0)
     fi
     k[-4]
    else
     pos={$1%$maxfreq}
     f[-4] I(#-2,0,$pos,0)
     k[-4]
    fi  
 fi
fi

Added rep_mode_color command. Next step is to figure out shape version of it, but don’t know how to do that yet. It is slow with blending unfortunately, but oh well.

Note: It doesn’t support CMYKA in case of mode per images rather than mode of multiple images per pixels. Tried supporting it, but turns out dynamic arrays were of two different sizes, so therefore cancelled it. Not to mention the huge amount of data to be processed in some cases.

Here’s a test:

$ sp cat rep_dupsdaxy 12,50,20,1,1,0 rep_mode_color ,, {w},{h},{d},1,!isnan(i#0)?255 replace_nan.. 0 a c drgba

image

The result does seem very interesting actually. Could in theory be used as a new artistic filter.

@David_Tschumperle Since I moved the discussion about func(), and func(a,b,c) to here. This is my reply:

I tested:

foo :
3,3,1,1
  f. "
  begin(
    $1?(
      func() = 3*x;
    ):(
      func() = pi + x;
    );
   );
   func();
  "

The above works.

However. this doesn’t seem to work:

#Error Output:1( fnx()=x+j; fny()=y; ):( fnx()=x+j; fny()=y+k; )' in expression '...): 0==1( fnx()=x+j; fny()=y; ):( fnx()=x+j; fny()=y+k; )...'.
foo2:
skip ${2=0},${3=0},${4=0},${5=0}
f "
begin(
const ang=pi*($4/180);
const ang2=pi*($5/180);
const cos_ang=cos(ang);
const sin_ang=sin(ang);
const cos_ang2=cos(ang2);
const sin_ang2=sin(ang2);
rot_x(a,b)=a*cos_ang-b*sin_ang;
rot_y(a,b)=a*sin_ang+b*cos_ang;
rot_y_alt(a,b)=a*sin_ang2+b*cos_ang2;
const blur=(round(abs($2)+1)-1)*2+3;
const qblur=sqr(blur);
const eblur=floor(blur/2);
const sblur=floor(blur/-2);
const coeff=$1/100;
const mult=qblur*coeff+eblur/round(abs($2)+1);
const choice=abs($3)<=2?abs($3):0;
j=k=0;
$3==5?(
 fnx()=x+rot_x(j,k);  
 fny()=y+rot_y_alt(j,k);
):
$3==4?(
 fnx()=x+rot_x(j,k);
 fny()=y+rot_y(j,k);
):
$3==3?(
 fnx()=x+rot_x(j,0);
 fny()=y+rot_y(j,0);
):
$3==2?(
 fnx()=x;
 fny()=y+k;
):
$3==1(
 fnx()=x+j;
 fny()=y;
):(
 fnx()=x+j;
 fny()=y+k;
);
);
sum=i*mult;
for(k=int(sblur),k<=eblur,k++,
    for(j=int(sblur),j<=eblur,j++,
        if(k!=0&j!=0,
        nx=fnx();
        ny=fny();
        nfcalc=i(nx,ny,z,c,0,1);
        sum-=nfcalc*coeff;
        );
    );
);
sum;
"

Hence, why I use if elif … fi and did fnx()="$fnx" to make it work.

EDIT: Just realized the error I did for so long. I forgot the ?. Sorry.