G'MIC math evaluator now deals with vectors (and matrices, and custom functions!)

So, I’m back here realizing that my Spiral Matrix Transform Filter - Spiral Matrix Transform Filter by Reptorian1125 · Pull Request #173 · dtschump/gmic-community · GitHub - has a serious issue. I realize that big dimensions simply don’t work. Like 4098x4098 according to my test. But different width and height bigger than that does work. I’m not sure if it the mathematics behind the code or there’s a big bug I discovered within vector processing.

EDIT: I find that it’s most likely has to do with spiralbw, some unexpected result shows up there.

EDIT: And using a different method posted at github bug issue, the bug still shows up. I noticed the bug starts exactly at 1 value after 4096x4096 which is 2^12,2^12. So, vector size cannot be greater than 16,777,216‬. The bigger the value after 4096x4096, the buggier it gets.

I’m trying to fix the reverse spiral transform but I don’t know what to do at all. I know that Reptorian’s already sorted out the spiral transform and I can condense it into this:

+spiralbw
f[0] "I(i#1%w,floor(i#1/w))"

However, I really don’t know how this works. As far as I know, this reads like how I imagine the reverse spiral transform to be. Since the index values for where the pixels after the spiral transform should be after the reverse spiral transform are found in the values of the pixels in the spiral, it should be possible to map them back. I don’t know how to do this at all since everything I’ve tried just gives me the spiral transform. How can I swap the pixels in the spiral transform with other pixels so that I can get the original back without having to construct another map or through using dynamic arrays?

More generally, if there is an image specifying the order in which pixels go, how do I transform back to the original order? I can actually get part of the way there in the spiral case if I do this:

rm
50,50,1,1,">begin(index=-1);index+=1;index"
+spiralbw
f[0] "I(i#1%w,floor(i#1/w))"
f[0] ">begin(index=-1);index+=1;
valatspiral=i[#1,index,2];
I[valatspiral]=I[#1,index];"
n 0,255

Edit: I’ve found a way to use the warp command to sort this whole mess out. G’MIC never fails to surprise, and this command is certainly extremely helpful! This script will return the original image (which I’ve assumed is already in the pipeline) through a spiral transform and a reverse spiral transform.

+spiralbw
f[0] "I(i#1%w,floor(i#1/w))"
# Here comes the magic!
100%,100%,1,2
f[2] "[i#1%w,floor(i#1/w)]"
warp[0] [2],2
k[0]

How does the spiral transform works?
Well, to give a clearer explanation, spiralbw gives the value in 0,1,2,3,4,5,6 all around in a spiral fashion. That’s what I call the coordinate map.

Every value of coordinate uses the top left coordinate of pixel which means it starts with 0,0

0 1 2 3
11 12 13 4
10 15 14 5
9 8 7 6

On the left side of i, that translate into

0 1 2 3
3 0 1 0
2 3 2 1
1 0 3 2

On the right side, it translate to

0 0 0 0
2 3 3 1
2 3 3 1
2 2 1 1

The left side is the x coordinate, and the right side is the y coordinate. Combining them gives this.

0,0 1,0 2,0 3,0
3,2 0,3 1,3 0,1
2,2 3,3 2,3 1,1
1,2 0,2 3,1 2,1

That what represents the output. The output is given by these coordinates. Use your mouse to get to where’s the original coordinate. It goes from left to right, top to bottom order. I used c++ to get the reverse transform since I had no idea how to reverse it otherwise.

@David_Tschumperle Warp does seem to return the original image, however there’s some missing information after 2^24 pixel. At least it doesn’t have the issue of the original dynamic array code.

See a close-up here -

image

I know about how the transform itself seems to work, but I kept on getting confused about which pixels it was selecting in trying to set some pixel values to the values of other pixels.

+spiralbw
f[0] "I[i#1]"
..
f. "i(#0,i0#1%w,floor(i0#1/w),0)=i"
1 Like

That’s another way of doing forward spiral transformation.

That being said, for the reverse: The warp approach while easy to do, has issues like a 1 layer dynamic array. I’m not sure how to transform the dynamic array to support double rather than float as reverse transform seem to a bug after 2^24th pixel. That’s the easiest method. The much harder method is to create a dynamic array that stacks layer after layer or overwriting a new layer with image dimension using values. Fill approach with loops, I’d imagine that would take a long time to calculate.

Hmm not sure I understand about the layers, but the example I just gave is forward and reverse (try removing the last instruction). Haven’t tested with large images though…

It does not work with large images. Error shows up after 2^24th pixel But, it’s a easier solution to do forward and reverse transform at once.

Reminds me about using single precision float insteat of double precision float for accumulating single precision float values (just a thought). The accumulation may saturate at a point…

Yes seems so (I think David mentioned it in another thread), but there should still be a way around it. I expect directly calculating 2D coordinates would solve the problem (spiralbw uses 1D indexes). Somebody just needs to work out the maths (tempting…) :slight_smile:

On the math part, the first part of the curves going up to down looks like this - min(max(slope1,slope2),max(slope3,slope4)). slope1, and slope2 are positive slope while the others are negative. So, if one could figure out generate 1d slope generated by two variables being width and height, you solve the problem. But this is just a idea.

What reminded me about the 2^24 stuff is the following:

Let’s make this loop over a 30 MP image where all pixels are equal:

float sum = 0.f;  // single precision float
for (i = 0; i < numberOfPixels; ++i) {
    sum += image[i];
}

If all pixels are equal, the accumulation will saturate at ~16 Mp.

Ingo
1 Like

Run out of time to finish/tidy this tonight, but here’s a proof of concept:

gcd_spiralbw :
sh. 0
eval. "
begin(P=[0,0];M=0;LX=w-1;LY=h-1;D=[1,0];R=[0,-1,1,0];sx=0;sy=0;xl=w-1);
if(
  (P[0]==LX && (P[1]==M || P[1]==LY))
  || (P[0]==M && P[1]==LY)
  , D=R*D);
if((P[0]==M && P[1]==M+1), D=R*D;M++;LX--;LY--);
P+=D;
sx++;
I(#-2,P[0],P[1])=[sx,sy,0];
if(sx==xl,sx=0;sy++);
" rm.

# example use (second warp is reverse):
gmic sp rose +f 0 gcd_spiralbw. warp.. . warp.. .,2
1 Like

image

That looks like a gui message (probably an image with opacity channel?), so I’m not surprised - the code is at “scratchpad” level. Somebody needs to work it to a better state for gui filter :slight_smile: . Sorry, don’t have time for at least a day now!

@garagecoder, you are truly awesome.
I’ve been itching to say it for a while, but I have to admit that your different G’MIC scripts are always very well done. You have no idea how much I love noise_poissondisk :heart_eyes::+1:
hugs!

1 Like

Also, some more info about the 2^{24} = 4096\times4096 limit.

The limit observed by @Reptorian was due to the fact he was using dynamic arrays, where the number of elements is stored as a float value (because it is stored directly as an image value, and images are float-valued in G’MIC).
As a float is 32bits and has “only” 24 bits of signficant digits (see Single-precision floating-point format - Wikipedia), this means it can store accurately integers (so the size of the dynamic array) only up to 2^{24}.

The fact is that such dynamic arrays have not been designed to store this amount of elements.
But there are plenty other solutions to deal with that issue.
If you stay inside the math evaluator (which works only with double, so with 53 bits of significant digits, Double-precision floating-point format - Wikipedia), you can manage integers up to 2^{53} = 67108864\times134217728 which becomes more than enough for 2d image resolutions :slight_smile:

This means, for instance, the code @garagecoder has proposed can be simplified a little bit, using an offset off rather than coordinates sx and sy:

gcd_spiralbw :
  sh. 0
  eval. "
begin(P=[0,0];M=0;LX=w-1;LY=h-1;D=[1,0];R=[0,-1,1,0];off=0);
if(
  (P[0]==LX && (P[1]==M || P[1]==LY)) || (P[0]==M && P[1]==LY)
  , D=R*D);
if((P[0]==M && P[1]==M+1), D=R*D;M++;LX--;LY--);
P+=D;

I(#-2,P[0],P[1])=[off%w,int(off/w),0];
++off;
"
  rm.

while being still valid for images up to resolution 67108864\times134217728.

Another solution would be to create variants of dynamic arrays in the math evaluator, with the array size coded in a variable rather than a pixel value. Each image dimension in G’MIC is coded as an unsigned int, so virtually up to 2^{32}.

1 Like

I can correct the vector size by forcing the input to have two channels and then getting rid of the 0 in the affected vector.

gcd_spiralbw_warp :
to_graya.
sh. 0
eval. "
begin(P=[0,0];M=0;LX=w-1;LY=h-1;D=[1,0];R=[0,-1,1,0];off=0);
if(
(P[0]==LX && (P[1]==M || P[1]==LY)) || (P[0]==M && P[1]==LY), D=R*D);
if((P[0]==M && P[1]==M+1), D=R*D;M++;LX--;LY--);
P+=D;
I(#-2,P[0],P[1])=[off%w,int(off/w)];
++off;
"
rm.

There are still a few issues which can be seen using this:

gcd_spiralbw_warp
f "[i0+i1*w,0]" to_gray
n 0,255

There is always a white dot towards the top left when I run it using the GIMP plugin (I haven’t tested the CLI version). For some image dimensions there are also line artefacts near the centre. It might be easier to figure out what’s happening by commenting out the second line. It looks like some sort of overshoot to me.

This made my day, thankyou! :blush:
G’MIC just seems to fit how I think extremely well. It’s nice the method is still on topic too, using vectors and matrices (but that’s coincidence really).

@Joan_Rake1 thanks for the testing, I’ll have a look at that when I can (I promised to look at other things and got distracted with this already!)

1 Like

Is it possible to convert a vector (e.g. after solving an equation) to a GRAYscale image?
[1,2,1,2,4,2,1,2,1] → (1,2,1,2,4,2,1,2,1)
It will help me visualize the results.

I need something like:
M=[1,2,1,2,4,2,1,2,1] # square matrix as a result of a computation
??? # Conversion to vector/matrix in the RGBA form, Mg=[1,1,1,0;2,2,2,0;1,1,1,0;…]
input[-1] {$Mg} # create an image from the vector Mg

Thanks for the help.