G'MIC exercises

Software G'MIC
668

Here is the exercise. Have fun!

Selectively filter an image based on a noise mask using the math interpreter. The masking and filtering should be done with a minimal number interpreter instances possible so as to limit switching back and forth between it and the regular interpreter, which seems to slow down the command in my experience.

The way the mask is generated is up to you. Let’s say we use the median filter; i.e. N=crop(…); med(N) but only with neighbourhoods that fulfill the mask criterion and not for every pixel. Doing this will save computational time because we are filtering less pixels and taking less neighbourhoods into consideration, which is beneficial for larger window sizes.

Bonus challenge

This can be a separate exercise but if you can write and use the vector median filter (VMF) instead of MF I would love it since VMF although processing intensive is quite often used to deal with noise. It is a perfect match for the so-called switching method described above.

669

Something like this perhaps:

gcd_filterbymask : skip ${1=5}
  rm sp cat +channels 0 f. "u>0.9" # dummy image and mask
  f.. "begin(const boundary=1;const S=$1;const N=2*S+1);
  ref(crop(x-S,y-S,z,c,N,N,1,1),K);
  i0#1 > 0 ? med(K) : i"

Edit: do you have an example of what you mean by “vector median filter”?

670

I think vmf means all of the result in crop per channels as a vector.

671

Hmm… sounds a bit like geometric median, which we actually have an example of further up in this very thread. That’s still present as the filter “Smooth [geometric-median]” and could easily be altered to use a mask. That is indeed a comparitively slow filter though.

672

u is a random number, right?

I prefer it all to be inside one fill or eval. I guess eval would be suitable for that if I prepared empty images in advance to be populated by the calculations.

VMF is an oldie, so everyone assumes everyone else knows what it is. :sweat_smile: There have been better faster implementations since. I found these examples to be helpful. (Others are too optimized or abstracted for me to understand what is going on.)

673

That does look like geometric median then, but where the resulting vector is one of the input vectors (a “medoid” perhaps?). If you find a good algorithm I could have a look to see if it’s likely to be quicker. That one looks like brute force / sorting based so probably even slower than geometric example above though.

The extra fill is indeed random - a fake noise mask. I think there’s a good chance creating your mask separately is still optimal, because keeping loops small is often crucial.

Edit: ah yes, good old k-means algorithms etc. - gmic has those amongst the palette commands

674

New year :partying_face: exercise: https://ietresearch.onlinelibrary.wiley.com/doi/pdf/10.1049/ipr2.12023.

Let’s add another wrinkle of complexity: find median of masked pixels. :slight_smile:

676

This is actually quite easy. Just put inf in the vector you get with crop(), at the locations of non-masked pixels, and sort the vectors.

V = crop(#0,x - 2,y - 2,z,c,5,5,1,1);  # 5x5 image neighborhood.
M = crop(#1,x - 2,y - 2,z,c,5,5,1,1);  # 5x5 mask (with values 0 or 1).
N = sum(M); # Number of masked pixels.
# for (m = 0, m<size(M), ++m, M[m]==0?(V[m] = inf)); # Put inf before sorting.
V*=(1e20*!M + 1);  # Do something similar
V2 = sort(V);
med = N?V2[int(N/2)]:i;  # Masked median value.
677

Elegant solution (after the edits).

Shouldn’t this be i since none of the pixels are selected? As is, we would retrieve the value from V[0], which is the top left corner of the neighbourhood.

678

Yes, you’re right. I’ll edit the post.

679

It seems to suggest that it’s also quite a slow algorithm! In the meantime, you might be interested in a median alternative: gcd_l2median_patch. It’s pushed to git, but I don’t know when that appears in filters these days. Only per-channel, not vector though…

It behaves like a slightly smoother median and allows non-integer kernel sizes. Perhaps not many uses, but it works well for the filter Upscale [Edge] (when it arrives):
gummyscaled

Edit: if you decide to use it, probably best to copy paste the code - it could change later!

680

Code looks familiar. Speaking of upscale…

I have been looking for robust round trip scaling methods for multi-scale processing. Scaling destroys the image so much that the rest of the processing isn’t worth it. Machine learning is in vogue but I don’t have the luxury. If anyone has any practical recommendations on where to start, that would be great.

Even down-scaling for the purposes of PlayRaws is a drag. I recall a night scene photo where the signposts were more illegible than others’ takes.

@David_Tschumperle I am confused by this. I would assume the output to be 0,1,2,3.

gmic e {V=[0,1,2,3];V[0,3]}  # 0,1,2
681

@afre - That can be seen as V[start of vector, number of vector]

682

I guess I misunderstood the docs. I thought the q was also an index. Will have to edit my files.

The sub-vector [X[p],X[p+s]…X[p+s*(q-1)]] (of size q) of a vector X is retrieved by X[p,q,s].

I don’t know. Looking at other non-G’MIC code samples, it seems that my interpretation is more common. If we were to go by length, then we would have to make sure we don’t assign one that is out of bounds. It will take more effort to find out when the step size isn’t 1.

683

The fact is that the size of a vector variable in a math expression is a const and must be known at expression parsing/compilation time.

  • If you allow a syntax like V[ind_start,ind_end], you cannot prevent V to have a variable size.
  • You could eventually force ind_start and ind_end to be const values to ensure the size of V is const as well (as size(V) would be ind_end - ind_start + 1).
  • But if you do that, you force ind_start to be a const, which is not cool, because in the current expression V[ind_start,size_V], the argument ind_start can be actually a variable itself (only the size of V must be a const, not the starting indice of the subvector).
684

Thanks for the explanation. It wasn’t a complaint; just confusion. It is easy to overlook even after reading the ref doc. Perhaps you could make the language more explicit. ATM, it doesn’t affect much of my code, at least from the regex search. But I do recall it tripping up others, including in my conversations with @garagecoder.

1 Like
685

I admit the reference doc is a bit tough. It’s open for contributions though :slight_smile:

686

I think I made a mistake - the new command is actually a vector L1 median approximate (i.e. geometric median). It might be possible to improve the accuracy using the algorithm above. If so, that will be a very much faster way to get a patch vector median, so I’m going to have a look at it later.

Edit:
@afre there were indeed a few mistakes. It’s definitely not a true median, but I think a useful approximate. It allows non-integer kernel sizes and operates on vectors rather than per channel (so no “false” colour).I’ve renamed it to “vecsmooth” with command gcd_vecsmooth (for want of a better name). It’s quite fast because only basic algebra and boxfilter are used. Example @ size = 3 for comparison:
elephantmedian

Edit 2: reduced the smoothing to a single iteration (consequently much faster and less memory hungry). That simplifies it to a variance based smooth, hence the new name: gcd_smooth_variance

688

I ended up working on parallel rep_dla and coded it according to the idea that if an area detect a point, then points will start moving. However, it only works for few iterations. :confused:

new_rep_dla 
new_rep_dla:
skip ${1=2},${2=10},${3=1},${4=0},${5=0},${6=1},${7=0}

m "dla_target_2s : n 0,1 s. c *"
m "dla_target_3s : n 0,1 s c add / 3"
m "dla_target_4s_plus : n 0,1 ts={s-1} s c add[^-1] /.. $ts *"
m "dla_target : if s==1 n 0,1 elif s==2 dla_target_2s elif s==3 dla_target_3s elif s>3 dla_target_4s_plus fi"
m "dla_check_variance : tv=0 repeat s sh. $> tv+={iv#-1} rm. if $tv u $tv return fi done u $tv"

if $4
 m "dla_expand : if $5 dilate_circ $5 else rm. fi"
else
 m "dla_expand : if $5 erode_circ $5 else rm. fi"
fi

include_dla_map=0

if $!==1
 dla_check_variance
 if ${}
  dla_target
  if iv#-1 include_dla_map=1 fi
 else
  channels 0 f 0
 fi
else
 include_dla_map=${"is_image_arg $7"}
 world_width=${-max_w}
 world_height=${-max_h}
 if $include_dla_map
  pass$7 0
  dla_target.
  if !iv#-1 include_dla_map=0 f. 0 fi
  r. $world_width,$world_height,100%,100%
 else
  $world_width,$world_height,1,1
 fi
fi

l. . f.. 0

 if $1>=2 noise_poissondisk.. $1,5
 else f.. x%2*y%2
 fi

 if $include_dla_map
  f. i>.5?1
  +dla_expand.
  if $!==3
   f[0] i#-2!=i#-1?i
   rm.
  else
   if $4
    f[0] i*(1-i#-1)
   else
    f[0] i*i#-1
   fi
  fi
  if !$6
   if $4
    *. .001
   else
    replace. 0,.999
   fi
  fi
 else
  f.. "begin(
   const ww=w-1;
   const hh=h-1;
   const cx=(w-1)/2;
   const cy=(h-1)/2;
  );
  xx=(x/ww-.5)*2;
  yy=(y/hh-.5)*2;
  sqrt(xx^2+yy^2)<1?i
  "
  if !$4 f. 1 fi
  eval "
   const ww=w-1;
   const hh=h-1;
   const cx=round((w-1)/2);
   const cy=round((h-1)/2);
   i(#-2,cx,cy)=0;
   i(#-1,cx,cy)=!i(#-1,cx,cy);
  "
 fi

 extract[0] "i==1",1 channels[0] 0,1
 
 if $4
  pt_eval=max(crop(#-1,xp-4,yp-4,0,0,9,9,1,1))
  if $3==0 mode=(i(#-1,xp-1,yp-1,0,0,0,2)||i(#-1,xp-1,yp+1,0,0,0,2))||(i(#-1,xp+1,yp-1,0,0,0,2)||i(#-1,xp+1,yp+1,0,0,0,2))
  elif $3==1 mode=(i(#-1,xp-1,yp,0,0,0,2)||i(#-1,xp+1,yp,0,0,0,2))||(i(#-1,xp,yp-1,0,0,0,2)||i(#-1,xp,yp+1,0,0,0,2))
  elif $3==2 mode=((i(#-1,xp-1,yp-1,0,0,0,2)||i(#-1,xp-1,yp+1,0,0,0,2))||(i(#-1,xp+1,yp-1,0,0,0,2)||i(#-1,xp+1,yp+1,0,0,0,2)))||((i(#-1,xp-1,yp,0,0,0,2)||i(#-1,xp+1,yp,0,0,0,2))||(i(#-1,xp,yp-1,0,0,0,2)||i(#-1,xp,yp+1,0,0,0,2)))
  elif $3==3 mode=altern?(i(#-1,xp-1,yp-1,0,0,0,2)||i(#-1,xp-1,yp+1,0,0,0,2))||(i(#-1,xp+1,yp-1,0,0,0,2)||i(#-1,xp+1,yp+1,0,0,0,2)):(i(#-1,xp-1,yp,0,0,0,2)||i(#-1,xp+1,yp,0,0,0,2))||(i(#-1,xp,yp-1,0,0,0,2)||i(#-1,xp,yp+1,0,0,0,2))
  fi
 else
  pt_eval=min(crop(#-1,xp-4,yp-4,0,0,9,9,1,1))<1
  if $3==0 mode=(i(#-1,xp-1,yp-1,0,0,0,2)<1||i(#-1,xp-1,yp+1,0,0,0,2)<1)||(i(#-1,xp+1,yp-1,0,0,0,2)<1||i(#-1,xp+1,yp+1,0,0,0,2)<1)
  elif $3==1 mode=(i(#-1,xp-1,yp,0,0,0,2)<1||i(#-1,xp+1,yp,0,0,0,2)<1)||(i(#-1,xp,yp-1,0,0,0,2)<1||i(#-1,xp,yp+1,0,0,0,2)<1)
  elif $3==2 mode=((i(#-1,xp-1,yp-1,0,0,0,2)<1||i(#-1,xp-1,yp+1,0,0,0,2)<1)||(i(#-1,xp+1,yp-1,0,0,0,2)<1||i(#-1,xp+1,yp+1,0,0,0,2)<1))||((i(#-2,xp-1,yp,0,0,0,2)<1||i(#-1,xp+1,yp,0,0,0,2)<1)||(i(#-1,xp,yp-1,0,0,0,2)<1||i(#-1,xp,yp+1,0,0,0,2)<1))
  elif $3==3 mode=altern?(i(#-1,xp-1,yp-1,0,0,0,2)<1||i(#-1,xp-1,yp+1,0,0,0,2)<1)||(i(#-1,xp+1,yp-1,0,0,0,2)<1||i(#-1,xp+1,yp+1,0,0,0,2)<1):(i(#-1,xp-1,yp,0,0,0,2)<1||i(#-1,xp+1,yp,0,0,0,2)<1)||(i(#-1,xp,yp-1,0,0,0,2)<1||i(#-2,xp,yp+1,0,0,0,2)<1)
  fi
 fi
 
 f.. ":begin(
   const ww=w#-1;
   const hh=h#-1;
   const lim_atmp=$2;
   const use_altern=$3==3;
   const dlmode=$4>0?1:0;
   newpos_x=[-1,-1,-1,0,0,1,1,1];
   newpos_y=[-1,0,1,-1,1,-1,0,1];
   altern=round(u(0,1));
  );
  attempts=0;
  cond=1;
  coords=I;
  xp=coords[0];
  yp=coords[1];
  "$pt_eval";
  do(
    if(!cond,break(););
    if("$pt_eval",
     do(
      if("$mode",
       i(#-1,xp%ww,yp%hh)=dlmode;
       cond=0;
       break();
      );
      pos_off=round(u(0,7));
      nox=newpos_x[pos_off];
      noy=newpos_y[pos_off];
      xp+=nox;
      yp+=noy;
     );
    );
  ,cond);
  "

endl

This part takes forever:

 f.. ":begin(
   const ww=w#-1;
   const hh=h#-1;
   const lim_atmp=$2;
   const use_altern=$3==3;
   const dlmode=$4>0?1:0;
   newpos_x=[-1,-1,-1,0,0,1,1,1];
   newpos_y=[-1,0,1,-1,1,-1,0,1];
   altern=round(u(0,1));
  );
  attempts=0;
  cond=1;
  coords=I;
  xp=coords[0];
  yp=coords[1];
  "$pt_eval";
  do(
    if(!cond,break(););
    if("$pt_eval",
     do(
      if("$mode",
       i(#-1,xp%ww,yp%hh)=dlmode;
       cond=0;
       break();
      );
      pos_off=round(u(0,7));
      nox=newpos_x[pos_off];
      noy=newpos_y[pos_off];
      xp+=nox;
      yp+=noy;
     );
    );
  ,cond);
  "

xp,yp are x and y pixel coordinate of points.

if pt_eval? is used to detect if a cropped area contains a value, and then it runs another do which is used to add pts to a image, and it’ll break away from the two do loop after. There has to be a better way of doing this as it leads to 100% CPU.

689

Looks complicated. I wonder if you could break it down into smaller pieces?