Fast median algorithm

[Siekmanski]

 

You could try a lowpass filter to get rid of some noise.

Maybe we can take advantage of atan(Gx/Gy), with a 3*3 matrix we get 6 directions back.
So it should be posible to find the direction of the pixel and eliminate all the pixels that don't point in the same direction.
This is theoretical of course. 

[guga]

You could try a lowpass filter to get rid of some noise.

Maybe we can take advantage of atan(Gx/Gy), with a 3*3 matrix we get 6 directions back.
So it should be posible to find the direction of the pixel and eliminate all the pixels that don't point in the same direction.
This is theoretical of course. 

Hi Marinus

Lowpass filter is what the crop_watermark routine actually does.  But, this seems to me a bit weird because we can simply revert the formula used in the matrix calculation to best find the numbers we are looking for.

For example, no matter if we are using Sobel, Scharr, Prewitt or whatever, the matrix of a 3x3 algorithm always resumes in:

Gx = A*(M1+M7-M3-M9) + D*(M4-M6) ; Where A and D are the 'magic numbers" used in the operator. For Sobel A = -1, D = -2. For Scharr, A = 47, D = 162 and so on. Those numbers are what forms the matrix and it´s transverse form.
Gy = A*(M1+M3-M7-M9) + D*(M2-M8)

Since the values at Gx and Gy re limited in between -1 to 1, i´m trying to identify the relation (properties) of M1 to M9 according to the values of the magic numbers A and D)

For example. The properties i found so far for Gx and Gy (after solving the math) was:

Gx Equations

M6 > 1/(-D)*((-D)*M4 - 1), M9 < -M1 + M3 - M7
M6 < 1/(-D)*((-D)*M4 - 1), M9 > -M1 + M3 - M7
M6 > 1/(-D)*((-D)*M4 - 1), M9>=1/(-A) (A*M1 + A*M3 + D*M4 + (-D)*M6 + A*M7 + 1)

Gy Equations
M8 > 1/(-D)*((-D)*M2 - 1), M9 > M1 + M3 - M7
M8 < 1/(-D)*((-D)*M2 - 1), M9 <= 1/(-A) (-A*M1 + (-D)*M2 + (-A)*M3 + A*M7 + D*M8 - 1)
M8 > 1/(-D)*((-D)*M2 - 1), M9 >= 1/(-A) (-A*M1 + (-D)*M2 + (-A)*M3 + A*M7 + D*M8 - 1)


If i´m correct, then the M1 value may not be used whatsoever, because for M9 do exist in Gx and Gy, The 1st part of all equations, shows an identify between
"M9 < -M1 + M3 - M7"
"M9 > -M1 + M3 - M7"
"M9 > M1 + M3 - M7"

So, the only logical solution  that M9 exists in all those equations is when M1 = 0. So, M1 may not be used anyway to calculate the final result of Gx and Gy. Which seems to be logical because M1 is the common pixel on all orientations (x, y). So it´s a bit weird we have to include the pixel itself to find how much influence the neighbours pixels are playing with him.

I´m not sure, but, if the logic is correct, then both equations can be resume to:

Gx = A*(0+M7-M3-M9) + D*(M4-M6)
Gy = A*(0+M3-M7-M9) + D*(M2-M8)

Therefore:

Gx = A*(M7-M3-M9) + D*(M4-M6)
Gy = A*(M3-M7-M9) + D*(M2-M8)

This should keep Gx and Gy on their own limits. What i´m trying to find is howM7, M9, M3, M4 etc are related to each other.  See what happens when M9 > or < then a given value or how they behave according to the values of A and D (that are constants) or trying to see if the value at M9 can be equal to M3 or M7 etc etc (so we can fix it directly simply adding or subtracting 1 from a given pixel (or 1/255 normalized), when we find pixels that are equal but they shouldn´t be) etc

This is what i´m currently trying to find if there is some property of the matrix we can use to force both equations to stay within the limits of -1 to 1 without having to be forced to use thresholds or other filters to fix this very same situation.

[daydreamer]

You could try a lowpass filter to get rid of some noise.

Maybe we can take advantage of atan(Gx/Gy), with a 3*3 matrix we get 6 directions back.
So it should be posible to find the direction of the pixel and eliminate all the pixels that don't point in the same direction.
This is theoretical of course. 

Hi Marinus

Lowpass filter is what the crop_watermark routine actually does.  But, this seems to me a bit weird because we can simply revert the formula used in the matrix calculation to best find the numbers we are looking for.

For example, no matter if we are using Sobel, Scharr, Prewitt or whatever, the matrix of a 3x3 algorithm always resumes in:

Gx = A*(M1+M7-M3-M9) + D*(M4-M6) ; Where A and D are the 'magic numbers" used in the operator. For Sobel A = -1, D = -2. For Scharr, A = 47, D = 162 and so on. Those numbers are what forms the matrix and it´s transverse form.
Gy = A*(M1+M3-M7-M9) + D*(M2-M8)

Since the values at Gx and Gy re limited in between -1 to 1, i´m trying to identify the relation (properties) of M1 to M9 according to the values of the magic numbers A and D)

For example. The properties i found so far for Gx and Gy (after solving the math) was:

Gx Equations

M6 > 1/(-D)*((-D)*M4 - 1), M9 < -M1 + M3 - M7
M6 < 1/(-D)*((-D)*M4 - 1), M9 > -M1 + M3 - M7
M6 > 1/(-D)*((-D)*M4 - 1), M9>=1/(-A) (A*M1 + A*M3 + D*M4 + (-D)*M6 + A*M7 + 1)

Gy Equations
M8 > 1/(-D)*((-D)*M2 - 1), M9 > M1 + M3 - M7
M8 < 1/(-D)*((-D)*M2 - 1), M9 <= 1/(-A) (-A*M1 + (-D)*M2 + (-A)*M3 + A*M7 + D*M8 - 1)
M8 > 1/(-D)*((-D)*M2 - 1), M9 >= 1/(-A) (-A*M1 + (-D)*M2 + (-A)*M3 + A*M7 + D*M8 - 1)


If i´m correct, then the M1 value may not be used whatsoever, because for M9 do exist in Gx and Gy, The 1st part of all equations, shows an identify between
"M9 < -M1 + M3 - M7"
"M9 > -M1 + M3 - M7"
"M9 > M1 + M3 - M7"

So, the only logical solution  that M9 exists in all those equations is when M1 = 0. So, M1 may not be used anyway to calculate the final result of Gx and Gy. Which seems to be logical because M1 is the common pixel on all orientations (x, y). So it´s a bit weird we have to include the pixel itself to find how much influence the neighbours pixels are playing with him.

I´m not sure, but, if the logic is correct, then both equations can be resume to:

Gx = A*(0+M7-M3-M9) + D*(M4-M6)
Gy = A*(0+M3-M7-M9) + D*(M2-M8)

Therefore:

Gx = A*(M7-M3-M9) + D*(M4-M6)
Gy = A*(M3-M7-M9) + D*(M2-M8)

This should keep Gx and Gy on their own limits. What i´m trying to find is howM7, M9, M3, M4 etc are related to each other.  See what happens when M9 > or < then a given value or how they behave according to the values of A and D (that are constants) or trying to see if the value at M9 can be equal to M3 or M7 etc etc (so we can fix it directly simply adding or subtracting 1 from a given pixel (or 1/255 normalized), when we find pixels that are equal but they shouldn´t be) etc

This is what i´m currently trying to find if there is some property of the matrix we can use to force both equations to stay within the limits of -1 to 1 without having to be forced to use thresholds or other filters to fix this very same situation.

If you use float pixels0.0,0.0,0.0,0.0black 1.0,1.0,1.0,1.0 =white while doing math instead of with 0.0-255,255,255,255
0.5,0.5,0.5 grey, keep all constants /variable between 0.0-1.0  maybe help
Because 0.5*0.5 =0.25, while 127.0*127.0 becomes lot bigger than 1.0

[guga]

Hi Daydreamer. The equations are already in the normalized version, just to make easier to try finding their properties. If i find a common formula to force the values at Gx and Gy be stuck inside -1 to 1, then we can also use a version using integers to try optimizing the routines.

I´m not so good in math, this is why i´m taking so much time to try finding this stuff. What i´m trying to do is finding a way reversing the routines to see if the math fits onto it. This is hard to do because if i succeed to fix a value for let´s say M9, then the value at M3 comes out of his place and vice-versa. That´s why trying to find the proper arithmetic involving this stuff is necessary, because i can see what causes the extrapolations, and what happens if M9 > M4-M3 etc etc.....

Even if i succeed to find it, i still need to see if it will work as expected. For what i´m seeing so far, the values in the imputed matrix can be changed/]updated to it best fits to the limits.  What i´ll also need to see is that if they needs to be changed only locally (so, per operation involving finding the sobel) or if they needs to be physically changed in the others pixels as well, so the next time the function enter on the loop to find the next position, it is already fixed by the previous loop and so on.

I tried to get help on a math group in facebook, but people simply couldn´t find the answer . So, i´m basically trying to do it by hand and also using www.wolframalpha.com to retrieve all possible solutions that can be used.

[mineiro]

If you use float pixels0.0,0.0,0.0,0.0black 1.0,1.0,1.0,1.0 =white while doing math instead of with 0.0-255,255,255,255
0.5,0.5,0.5 grey, keep all constants /variable between 0.0-1.0  maybe help
Because 0.5*0.5 =0.25, while 127.0*127.0 becomes lot bigger than 1.0

Hi Daydreamer. The equations are already in the normalized version, just to make easier to try finding their properties. If i find a common formula to force the values at Gx and Gy be stuck inside -1 to 1, then we can also use a version using integers to try optimizing the routines.

I'm not following this topic, when I read Daydreamer answer and guga comment Arithmetic Coding comes to my mind. Maybe can be an answer.
The output is between 0 (low) and 1 (high).
So renormalize/rescaling can be done by binary.
0.5 == 1/2, 0.25 == 1/4 and 0.75 == 3/4
Values get converging between high and low; when a situation like
1/4 and 3/4 happens (the MSB aren't equal) is necessary to readjust (rescaling). If they are the same you output then and rescale so value can fit in a range (register size as an example).
The nice of this is that you can separate model (symbols frequency table as an example) from codification.
Zero is seen as an infinite zeros and 1.0 is seen as infinite ones, so you can do this by using binary only.
So, 1.0 can be the same as 255 or FFFFFFFF... . To avoid overflow you can deal with 7FFFFF... .

The most easy to follow text about AC that I have found and have C source code easy to understand was:
http://www.sable.mcgill.ca/~ebodde/pubs/sable-tr-2007-5.pdf
This text starts by using floating point numbers and in next chapter deals with binary numbers only; and explain why is necessary rescaling.
Well, if this works to you, you can check "range encoding" or "srANS"(assymetrical numeral system, Jarek Duda) thats public domain.

[guga]

Hi Mineiro

I´ll read the code, but i still clueless about the proper answer. I´mmpost here in 2 parts, because the text became long. What i found so far is:

I - properties of the matrix involving Sobel and Scharr.


    Scharr Operators works similar to sobel.
   
    In Sobel we get this matrices   
    Sobelx
        -1 0 1
        -2 0 2
        -1 0 1

    Sobely
        -1 -2 -1
        0   0  0
        1   2  1
   
    Scharr Matrices are given by:

    ScharrX
        47  0 -47
        162 0 -162
        47  0 -47

    ScharrY
        47   162   47
        0     0     0
        -47 -162  -47

    Other matrix of scharr is:

    ScharrX
        3  0 -3
        10 0 -10
        3  0 -3

    ScharrY
        3   10   3
        0    0   0
        -3 -10  -3
   

    Laplacian is simple   
        0   -1   0
        -1   4  -1
        0   -1   0


The convolution of the matrices are done like ths:

Matrix Values

A B C
D E F
G H I

Pixels Positions

M1 M2 M3
M4 M5 M6
M7 M8 M9

Convolution is generally done like this:
Convolute = A*M1+D*M4+G*M7+ B*M2+E*M5+H*M8 +C*M3+F*M6+I*M9

In Sobel Operator and Scharr we have one line and one row emptied (with 0) and the opposed side has the exact same values as the other,
but with the sign switched. Thus, to make the math more easy those can be written as:

For Gx
B, E, H = 0
C = -A
F = -D
I = -G

therefore, taking the left side as a reference:

Gx = A*M1+D*M4+G*M7+ B*M2+E*M5+H*M8 +C*M3+F*M6+I*M9

turns onto:

Gx = A*M1+D*M4+G*M7+ 0*M2+0*M5+0*M8 +(-A)*M3+(-D)*M6+(-G)*M9
Gx = A*M1+D*M4+G*M7 -A*M3 -D*M6 -G*M9

We also know that for Gx, A = G (Both are -1 for sobelX and 47 or 3 for Scharr), therefore:

Sobel, A = -1, D = -2
Scharr, A = 47, D = 162

Gx = A*M1+D*M4+A*M7 -A*M3 -D*M6 -A*M9
Gx = A*M1+A*M7-A*M3-A*M9 +D*M4 -D*M6
Gx = A*(M1+M7-M3-M9) + D*(M4-M6)

Since Gx can be a number valid in between -1 to 1 we can define a equation like this:

-1>= A*(M1+M7-M3-M9) + D*(M4-M6) <= 1

Solving it result in:

A*(M1 + M7) + D*M4 + 1 <= A*(M3 + M9) + D*M6

solving again for A, we have:

A*(M1 + M7) - A*(M3 + M9) <=  D*M6 - D*M4 - 1
A*(M1 + M7 - M3 - M9) <=  D*(M6-M4)- 1

A <= (D*(M6-M4) - 1)/ (M1+M7-M3-M9)

If A and D are negative, upon the above equation, we have the following solutions:

M6 = M4 + 1/D, M9 > M1 - M3 + M7
M6 = M4 + 1/D, M9 < M1 - M3 + M7
M6 > M4 + 1/D, M9 > M1 - M3 + M7
M6 > M4 + 1/D, M9 <= (A*(M1+M7-M3) + D*(M4-M6) + 1)/A
M6 < M4 + 1/D, M9 >= (A*(M1+M7-M3) + D*(M4-M6) + 1)/A

Ranaming:

Delta1 = M6-M4
Delta2 = M1+M7-M3

we have:
Delta1 = 1/D, M9 > Delta2
Delta1 = 1/D, M9 < Delta2
Delta1 > 1/D, M9 > Delta2
Delta1 > 1/D, M9 <= (A*Delta2 - D*Delta1 + 1)/A
Delta1 < 1/D, M9 >= (A*Delta2 - D*Delta1 + 1)/A

We now have 2 identities
Delta1 > 1/D, M9 > Delta2
Delta1 > 1/D, M9 <= (A*Delta2 - D*Delta1 + 1)/A

So, for each Delta1 > 1/D, we have the following properties:

M9 > Delta2
A*(M9-Delta2) <= 1 - D*Delta1

Renaming again
K = M9-Delta2, we have:

K > 0
A*K <= 1 - D*Delta1

____________________________________________________________________
____________________________________________________________________

On his turn, Gy, has the following properties:

For Gy
D, E, F = 0
A = -G
B = -H
C = -I

Gy = A*M1+D*M4+G*M7+ B*M2+E*M5+H*M8 +C*M3+F*M6+I*M9
turns onto:


Gy = A*M1+0*M4+(-A)*M7+ B*M2+0*M5+(-B)*M8 +C*M3+0*M6+(-C)*M9
Gy = A*M1-A*M7+ B*M2-B*M8 +C*M3-C*M9
We also know that for Gy, C = A (Both are -1 for sobelY and 47 or 3 for Scharr), therefore:
Gy = A*M1-A*M7+ B*M2-B*M8 +A*M3-A*M9
Gy = A*M1-A*M7 +A*M3-A*M9 +B*M2-B*M8
Gy = A*(M1+M3-M7-M9) + B*(M2-M8)

Since Gy can also be a number valid in between -1 to 1 we can define a equation like this:

-1>= A*(M1+M3-M7-M9) + B*(M2-M8) <= 1

Solving it result in:

A*(M1 + M3) + B*M2 + 1 <= A*(M7 + M9) + B*M8

solving again for A, we have: (Remembering that B = D. I.e: The transversed value of SobelGx)

A <= (D*(M8-M2)-1)/(M1+M3-M7-M9)

Since A and D are negative, upon the above equation, we have these solutions:

1/(-D)*((-D)*M2 - 1)

M8 = M2 + 1/D, M9 > M1 + M3 - M7
M8 = M2 + 1/D, M9 < M1 + M3 - M7
M8 > M2 + 1/D, M9 > M1 + M3 - M7
M8 > M2 + 1/D, M9 <= (A*(M1+M3-M7) + D*(M2-M8) + 1) / A
M8 < M2 + 1/D, M9 >= (A*(M1+M3-M7) + D*(M2-M8) + 1) / A

Ranaming
DeltaY1 = M8-M2
DeltaY2 = M1+M3-M7

we have:
DeltaY1 = 1/D, M9 > DeltaY2
DeltaY1 = 1/D, M9 < DeltaY2
DeltaY1 > 1/D, M9 > DeltaY2
DeltaY1 > 1/D, M9 <= (A*DeltaY2 - D*DeltaY1 + 1) / A
DeltaY1 < 1/D, M9 >= (A*DeltaY2 - D*DeltaY1 + 1) / A

We now have 2 identities
DeltaY1 > 1/D, M9 > DeltaY2
DeltaY1 > 1/D, M9 <= (A*DeltaY2 - D*DeltaY1 + 1) / A

So, for each Delta1 > 1/D, we have the following properties:

M9 > DeltaY2
A*(M9-DeltaY2) <= 1 - D*DeltaY1

Renaming again
Ky = M9-DeltaY2, we have:

Ky > 0
A*Ky <= 1 - D*DeltaY1


Those are the properties i found for the whole matrix operation. It seems correct, but don´t know yet how to apply. For example, M9 or M8 can be negative values ? If they can, then how is it possible since they are, in fact only integers from 0 to 255 (i.e: 0 to 1)?,. Also what happens, for example with the solutions i´ve found like:

M8 = M2 + 1/D, M9 > M1 + M3 - M7

It implies that M9 on this case is always bigger then M1+M3-M7 (also all positive values) only if M8 = M2+1/D. Is that correct ?

Do i need a chunk of Conditional If to evaluate when those conditions are met for each solution ?


Now...Part 2 i got this....(next post)

[guga]

II - Directly reverting the equation at G

I tried to revert the final value of G using the equatiosn i found in the previous post. This is what i got.

____________________________________________________________________
____________________________________________________________________

let´s now calculate the total G from those properties and see what we found

Since G^2 = Gx^+Gy^2, we have:

Gx = A*(M1+M7-M3-M9) + D*(M4-M6)
Gy = A*(M1+M3-M7-M9) + D*(M2-M8)

The maximum value of G^2 = 1, and the minimum is 0 therefore:

0<= (A*(M1+M7-M3-M9) + D*(M4-M6))^2 + (A*(M1+M3-M7-M9) + D*(M2-M8))^2 <= 1

This equation is hard to solve. So, let´s rename the variables for easier understanding and use the proper signs
for A and D (in case both are negative) to see what we found.

renaming:
x = M1
y = M7
z = M3
w = M9
k = M4
t = M6
g = M2
h = M8

we end up with this:

0<= (A*(x+y-z-w) + D*(k-t))^2 + (A*(x+z-y-w) + D*(g-h))^2 <= 1

For Sobel:
A = -1, D = -2
For Scharr:
A = 47, D = 162

Let´ use some random values keeping the sign to see if we find a solution:

Using A = -1, D = -2, we have:

0<= (-1*(x+y-z-w) + (-2)*(k-t))^2 + (-1*(x+z-y-w) + (-2)*(g-h))^2 <= 1

And the possible solutions are:
x = 1/2 (-2 g + 2 h - 2 k + 2 t + 2 w - sqrt(2)), z = 1/2 (-2 g + 2 h + 2 k - 2 t + 2 y)
x = 1/2 (-2 g + 2 h - 2 k + 2 t + 2 w + sqrt(2)), z = 1/2 (-2 g + 2 h + 2 k - 2 t + 2 y)


Using A = -17, D = -49, we have:
0<= (-17*(x+y-z-w) + (-49)*(k-t))^2 + (-17*(x+z-y-w) + (-49)*(g-h))^2 <= 1

And the possible solutions are:
x = 1/34 (-49 g + 49 h - 49 k + 49 t + 34 w - sqrt(2)), z = 1/34 (-49 g + 49 h + 49 k - 49 t + 34 y)
x = 1/34 (-49 g + 49 h - 49 k + 49 t + 34 w + sqrt(2)), z = 1/34 (-49 g + 49 h + 49 k - 49 t + 34 y)


Now with positive values:

Using A = 17, D = 49, we have:
0<= (17*(x+y-z-w) + (49)*(k-t))^2 + (17*(x+z-y-w) + (49)*(g-h))^2 <= 1

And the possible solutions are:
x = 1/34 (-49 g + 49 h - 49 k + 49 t + 34 w - sqrt(2)), z = 1/34 (-49 g + 49 h + 49 k - 49 t + 34 y)
x = 1/34 (-49 g + 49 h - 49 k + 49 t + 34 w + sqrt(2)), z = 1/34 (-49 g + 49 h + 49 k - 49 t + 34 y)

Using A = 1, D = 2, we have:
0<= (1*(x+y-z-w) + (2)*(k-t))^2 + (1*(x+z-y-w) + (2)*(g-h))^2 <= 1

x = 1/2 (-2 g + 2 h - 2 k + 2 t + 2 w - sqrt(2)), z = 1/2 (-2 g + 2 h + 2 k - 2 t + 2 y)
x = 1/2 (-2 g + 2 h - 2 k + 2 t + 2 w + sqrt(2)), z = 1/2 (-2 g + 2 h + 2 k - 2 t + 2 y)

Now mixing the signs:

Using A= 59, D = -30
0<= (59*(x+y-z-w) + (-30)*(k-t))^2 + (59*(x+z-y-w) + (-30)*(g-h))^2 <= 1

And the possible solutions are:
x = 1/118 (30 g - 30 h + 30 k - 30 t + 118 w - sqrt(2)), z = 1/118 (30 g - 30 h - 30 k + 30 t + 118 y)
x = 1/118 (30 g - 30 h + 30 k - 30 t + 118 w + sqrt(2)), z = 1/118 (30 g - 30 h - 30 k + 30 t + 118 y)

And switching the signs:
Using A= -59, D = 30
0<= (-59*(x+y-z-w) + (30)*(k-t))^2 + (-59*(x+z-y-w) + (30)*(g-h))^2 <= 1

And the possible solutions are:
x = 1/118 (30 g - 30 h + 30 k - 30 t + 118 w - sqrt(2)), z = 1/118 (30 g - 30 h - 30 k + 30 t + 118 y)
x = 1/118 (30 g - 30 h + 30 k - 30 t + 118 w + sqrt(2)), z = 1/118 (30 g - 30 h - 30 k + 30 t + 118 y)

Ok, we have a standard solution here. Let´s try identify them:

Sign1 = -1 when A and D have the same sign. So, both are positive or both are negative.
Sign1 = 1 when A and D have different signs. So, one is positive and other negative. (or vice-versa)

Sign2 = 1 when A and D have the same sign. So, both are positive or both are negative.
Sign2 = -1 when A and D have different signs. So, one is positive and the other negative. (or vice-versa)

|xxxx|  = Absolute value of a number (all numbers are turned onto positive)

x = |1/(2*A)| * ( (Sign1)*|D|*g + (Sign2)*|D|*h + (Sign1)*|D|*k + (Sign2)*|D|*t + |(2*A)|*w - sqrt(2) ), z = |1/(2*A)| * ( (Sign1)*|D|*g + (Sign2)*|D|*h + (Sign2)*|D|*k + (Sign1)*|D|*t + |(2*A)|*y )
x = |1/(2*A)| * ( (Sign1)*|D|*g + (Sign2)*|D|*h + (Sign1)*|D|*k + (Sign2)*|D|*t + |(2*A)|*w + sqrt(2) ), z = |1/(2*A)| * ( (Sign1)*|D|*g + (Sign2)*|D|*h + (Sign2)*|D|*k + (Sign1)*|D|*t + |(2*A)|*y )


Putting the proper labels back we have:

M1 = |1/(2*A)| * ( (Sign1)*|D|*M2 + (Sign2)*|D|*M8 + (Sign1)*|D|*M4 + (Sign2)*|D|*M6 + |(2*A)|*M9 - sqrt(2) ), M3 = |1/(2*A)| * ( (Sign1)*|D|*M2 + (Sign2)*|D|*M8 + (Sign2)*|D|*M4 + (Sign1)*|D|*M6 + |(2*A)|*M7 )
M1 = |1/(2*A)| * ( (Sign1)*|D|*M2 + (Sign2)*|D|*M8 + (Sign1)*|D|*M4 + (Sign2)*|D|*M6 + |(2*A)|*M9 + sqrt(2) ), M3 = |1/(2*A)| * ( (Sign1)*|D|*M2 + (Sign2)*|D|*M8 + (Sign2)*|D|*M4 + (Sign1)*|D|*M6 + |(2*A)|*M7 )


In practice, the above equations result in a couple of different values:

A - When both signs of A and D are negative, we have:

    M3 = |1/(2*A)| * ( (Sign1)*|D|*M2 + (Sign2)*|D|*M8 + (Sign2)*|D|*M4 + (Sign1)*|D|*M6 + |(2*A)|*M7 )
    M3 = |1/(2*A)| * ( -|D|*M2 + |D|*M8 + |D|*M4 - |D|*M6 + |(2*A)|*M7 )

    let´s disconsider the Abs sign (we only need the values of A and D, not their signs), so it lead us onto:
    M3 = 1/(2*A) * ( D*(M8-M2) + D*(M4-M6) + 2*A*M7 )

    M1 = |1/(2*A)| * ( (Sign1)*|D|*M2 + (Sign2)*|D|*M8 + (Sign1)*|D|*M4 + (Sign2)*|D|*M6 + |(2*A)|*M9 - sqrt(2) )
    M1 = 1/(2*A) * ( -|D|*M2 + |D|*M8 -|D|*M4 + |D|*M6 + 2*A*M9 - sqrt(2) )

    let´s disconsider the Abs sign (we only need the values of A and D, not their signs), so it lead us onto:
    M1 = 1/(2*A) * ( D*(M8-M2) - D*(M4-M6) + 2*A*M9 - sqrt(2) )

Gx = A*(M1+M7-M3-M9) + D*(M4-M6)



Those are what i´ve found so far, but when i testd, i´ve got weird results.;..Then i stopped this part of the calculations yesterday night to rest a little.

It seems correct, but, the problem is that, when i put them on excel...The result is always Sqrt(2)/2, no matter what are the values of the input. If I accept M3 and M1 to use values calculated from Item A above, even if we find thing like -458 etc.... when applying these values, the resultant G will always be sqrt(2)/2. And we will end up with values of M3 and M1 outside of the limits of 0 to 255 (0 to 1.0)


What i´m missing here ?

[Siekmanski]

Hi guga,

You have made it very complicated.
If you want a certain range back, normalize the edge-detection matrices ( sobel, scharr etc. ),
There are only 6 values for the convolution calculations for Gradient X and 6 for Gradient Y.

For the Sobel Matrix the range = -1020 to 1020
If you want the range -1 to 1, divide the Sobel Matrix members with 1020. ( or multiply with 0.0009803921)
It's handy to have the gray values as real4 in memory.

    Laplacian is simple   
        0   -1   0
        -1   4  -1
        0   -1   0

The Laplacian matrix looks more like a lowpass filter.

[guga]

-765 to 765 ?

How did you found that range ?  Does it have a formula to retrieve it ? Is it the same range as in Scharr or other operators ?

[Siekmanski]

Sorry, just corrected my previous post.
It should be -1020 to 1020.

max. value pixel = 255

-1,-2,-1 = -1*255 + -2*255 + -1*255 = -1020
1,2,1 = 1*255 + 2*255 + 1*255 = 1020

[guga]

Thanks...Now i see what you were talking about normalizing the values.

Indeed, doing this way (rather then trying to find a equation) does the job and forces the values to stay within their own limits. I made a test now, and the resultant image is too dark, but, that´s not a problem, because it seems that the edges are all there properly identified and the limits were all respected . If i want eventually to better see the result, all i have to do is make the result be a bit brighter calculating the difference between the max and minimum sobel found on each image and normalize again (as JJ did to calculate the median). But, normalizing again probably won´t be needed in the scope of the watermark remover, because the shape of the watermark can be identified as well using this method. So, it don´t matter if the resultant image is too dark, because the shape can be created anyway and fixed on the way you showed :) 

The good thing is that now we can calculate the proper normalization ratio using different operators, such as sobel or Scharr etc. To calculate the normalization ratio it is as simple as this:

Ratio = 1/((|ValA|*2+|ValD|)*255)

For Sobel:
ValA = -1 (in abs, only the number and not the sign)
ValD = -2 (in abs, only the number and not the sign)

For Scharr:
ValA = 47 (in abs, only the number and not the sign)
ValD = 162 (in abs, only the number and not the sign)


And the apply those values to Gx or Gy and multiply by 255 to we get the proper normalization back in respect to -255 to 255, right ?


Gx = (A*(M1+M7-M3-M9) + D*(M4-M6)) * Ratio * 255
Gy = (A*(M1+M3-M7-M9) + D*(M2-M8)) * Ratio * 255

On this way, Gx and Gy will be normalized to -1 to 1 and also stay within -255 to 255 if we multiply the result by 255, right ?


Or, simply doing this, right (When M1, M2...xxx are the integers)?

Ratio = 255/((|ValA|*2+|ValD|)*255) = 1/((|ValA|*2+|ValD|)

[mineiro]

Matrix Values

A B C
D E F
G H I

Pixels Positions

M1 M2 M3
M4 M5 M6
M7 M8 M9

Convolution is generally done like this:
Convolute = A*M1+D*M4+G*M7+ B*M2+E*M5+H*M8 +C*M3+F*M6+I*M9

I don't undestand, I need read all this post.
My thinking of whats happening was divide to conquer(1/2, 1/4 and 3/4, ...):
123456789
1*X^8 + 2*X^7 + 3*X^6 + ... + 9*X^0 . Changing base X can be an option.
So, forgot my comment, I was thinking other things. Sorry.

[guga]

Hi Marinus

Tested it and the result is amazing. Even though the resultant image is too dar, after i equalized it, it enhances all the edges nearly to perfection.

I´ll make some adjusts to better see the edges (lighter them up a bit) to see what is the result without equalizing. (I equalized only to see the results, it´s not part ot the watermark remover algo)

So, adjusting the sobel/Scharr etc to it fits to their own limits are the way to go. No matter if we us it for tyhe watermark remover or sobel itself. The edges are thin, but they are all there. Also, there´s no need to smooth the image as it is commonly done with other apps, because it can detect everything regardless the quality of the image.

[guga]

Sorry, just corrected my previous post.
It should be -1020 to 1020.

max. value pixel = 255

-1,-2,-1 = -1*255 + -2*255 + -1*255 = -1020
1,2,1 = 1*255 + 2*255 + 1*255 = 1020

Hi Marinus...Another question.

This is to compute the ratio of a 3x3 matrix. How to do it for 5x5 matrix or 9x9 etc ? You also need to sum up the values of the cols ? But, what if we have 2 or more cols after the midle one ? we sum them all too ?

Like this:

   -5  -4  0   4   5
   -8 -10  0  10   8
  -10 -20  0  20  10
   -8 -10  0  10   8
   -5  -4  0   4   5

Is correct to do this ?
2*(5*255 + 8*255) + 10*255 + 2*(4*255+10*255) + 20*255
which resumes to:
((5+8+4+10)*2  + (10+20) ) * 255 = 21420

[Siekmanski]

It's as simple as this, no matter the size of the matrices:

e.g. a maximum edge:

Code Select Expand


Do the Matrix convolution with the maximum pixel value (255).

         [ 255 0 0 ]
Pixels = [ 255 0 0 ]
         [ 255 0 0 ]


     [+1  0 -1 ]            [ 255 0 0 ]
Gx = [+2  0 -2 ] * Pixels = [ 510 0 0 ] = 1*255 + 0*0 + -1*0 + 2*255 + 0*0 + -2*0 + 1*255 + 0*0 + -1*0 = 1020 ( we found a maximum edge on the X-axis )
     [+1  0 -1 ]            [ 255 0 0 ]


     [+1 +2 +1 ]            [ 255 0 0 ]
Gy = [ 0  0  0 ] * Pixels = [   0 0 0 ] = 1*255 + 2*0 + 1*0 + 0*255 + 0*0 + 0*0 + -1*255 + -2*0 + -1*0 = 0 ( we found no edge on the Y-axis )
     [-1 -2 -1 ]            [-255 0 0 ]

As you can see here: 1020 is a maximum edge but also -1020 is a maximum edge.

         [ 0 0 255 ]
Pixels = [ 0 0 255 ] * Gx = -1020
         [ 0 0 255 ]


G = sqrt( Gx * Gx + Gy * Gy ) = 1020

The result of an Edge is now between 0 - 1020

If your result is nearer to 0, the edge will be weaker.

The difference between the minimum and maximum of GradientX and GradientY is the range.
The range for Sobel is -1020 to 1020.

Normalize the Sobel Matrix to the maximum range you need.
If you want the range -255 to 255,
Divide all Sobel Matrix members by 1020 * 255, write it back to a New Matrix ( normalized to your chosen range )

If you want to calculate the gradient's direction, atan( Gy / Gx )
You could also consider to calculate a range between -PI/2 and PI/2 ( -90 to 90 degrees )

Note: as you can see, you only need to convolute 6 values per Sobel matrix.