// return the median value in a vector of 27 floats pointed to by a
float heapMedian3( float *a )
{
float left[14], right[14], median, *p;
unsigned char nLeft, nRight;
// pick first value as median candidate
p = a;
median = *p++;
nLeft = nRight = 1;
for(;;)
{
// get next value
float val = *p++;
// if value is smaller than median, append to left heap
if( val < median )
{
// move biggest value to the heap top
unsigned char child = nLeft++, parent = (child - 1) / 2;
while( parent && val > left[parent] )
{
left[child] = left[parent];
child = parent;
parent = (parent - 1) / 2;
}
left[child] = val;
// if left heap is full
if( nLeft == 14 )
{
// for each remaining value
for( unsigned char nVal = 27 - (p - a); nVal; --nVal )
{
// get next value
val = *p++;
// if value is to be inserted in the left heap
if( val < median )
{
child = left[2] > left[1] ? 2 : 1;
if( val >= left[child] )
median = val;
else
{
median = left[child];
parent = child;
child = parent*2 + 1;
while( child < 14 )
{
if( child < 13 && left[child+1] > left[child] )
++child;
if( val >= left[child] )
break;
left[parent] = left[child];
parent = child;
child = parent*2 + 1;
}
left[parent] = val;
}
}
}
return median;
}
}
// else append to right heap
else
{
// move smallest value to the heap top
unsigned char child = nRight++, parent = (child - 1) / 2;
while( parent && val < right[parent] )
{
right[child] = right[parent];
child = parent;
parent = (parent - 1) / 2;
}
right[child] = val;
// if right heap is full
if( nRight == 14 )
{
// for each remaining value
for( unsigned char nVal = 27 - (p - a); nVal; --nVal )
{
// get next value
val = *p++;
// if value is to be inserted in the right heap
if( val > median )
{
child = right[2] < right[1] ? 2 : 1;
if( val <= right[child] )
median = val;
else
{
median = right[child];
parent = child;
child = parent*2 + 1;
while( child < 14 )
{
if( child < 13 && right[child+1] < right[child] )
++child;
if( val <= right[child] )
break;
right[parent] = right[child];
parent = child;
child = parent*2 + 1;
}
right[parent] = val;
}
}
}
return median;
}
}
}
}
/* This source code is in the public domain */
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }
# Note: The code needs more than 2 elements to work
float lefselect(float a[], const int n, const int k) {
int l=0, m = n-1, i=l, j=m;
float x;
while (l<m) {
if( a[k] < a ) F_SWAP(a,a[k]);
if( a[j] < a ) F_SWAP(a,a[j]);
if( a[j] < a[k] ) F_SWAP(a[k],a[j]);
x=a[k];
while (j>k & i<k) {
do i++; while (a<x);
do j--; while (a[j]>x);
F_SWAP(a,a[j]);
}
i++; j--;
if (j<k) {
while (a<x) i++;
l=i; j=m;
}
if (k<i) {
while (x<a[j]) j--;
m=j; i=l;
}
}
return a[k];
}
;
; build: asmc64 -q -pe median.asm
;
include stdio.inc
include stdlib.inc
include tchar.inc
type typedef long_t ; item size
compare proto :ptr, :ptr {
mov eax,[rcx]
sub eax,[rdx]
sbb ecx,ecx
sbb ecx,-1
and eax,eax
cmovnz eax,ecx
}
memxchg proto :ptr, :ptr {
mov eax,[rcx]
mov r8d,[rdx]
mov [rdx],eax
mov [rcx],r8d
}
.data
array dd 1, 18, 0, 0, 0, 0, 2, 2, 5, 99, 600, 785, 758 , 17, 11, 21, 23, 24, 66, 21, 14, 57, 17
.code
median proc uses rsi rdi rbx p:ptr, n:size_t
local result:ptr
mov rax,rdx
.if rax > 1
mov rsi,rcx
dec rax
mov ecx,type
mul rcx
lea rdi,[rsi+rax]
xor r9d,r9d ; stack level
lea rax,[rdi+rcx] ; middle from (hi - lo) / 2
sub rax,rsi
xor edx,edx
div rcx
shr rax,1
mul rcx
add rax,rsi
mov result,rax
.while 1
mov rcx,type
lea rax,[rdi+rcx]
sub rax,rsi
.ifnz
xor edx,edx
div rcx
shr rax,1
mul rcx
.endif
lea rbx,[rsi+rax]
.ifsd compare(rsi, rbx) > 0
memxchg(rsi, rbx)
.endif
.ifsd compare(rsi, rdi) > 0
memxchg(rsi, rdi)
.endif
.ifsd compare(rbx, rdi) > 0
memxchg(rbx, rdi)
.endif
mov r10,rsi
mov r11,rdi
.while 1
add r10,type
.if r10 < rdi
.continue .ifsd compare(r10, rbx) <= 0
.endif
.while 1
sub r11,type
.break .if r11 <= rbx
.break .ifsd compare(r11, rbx) <= 0
.endw
mov rcx,r11
mov rax,r10
.break .if rcx < rax
memxchg(rcx, rax)
.if rbx == r11
mov rbx,r10
.endif
.endw
add r11,type
.while 1
sub r11,type
.break .if r11 <= rsi
.break .ifd compare(r11, rbx)
.endw
mov rdx,r10
mov rax,r11
sub rax,rsi
mov rcx,rdi
sub rcx,rdx
.ifs rax < rcx
mov rcx,r11
.if rdx < rdi
push rdx
push rdi
inc r9d
.endif
.if rsi < rcx
mov rdi,rcx
.continue
.endif
.else
mov rcx,r11
.if rsi < rcx
push rsi
push rcx
inc r9d
.endif
.if rdx < rdi
mov rsi,rdx
.continue
.endif
.endif
.break .if !r9d
dec r9d
pop rdi
pop rsi
.endw
mov rcx,result
mov eax,[rcx]
.endif
ret
median endp
main proc
printf( "median: %d\n", median( &array, lengthof(array) ) )
ret
main endp
end _tstart
Here's one, (q)sort and grab the middle
without using sorting algorithms
float Dansby_kth(float array[], const unsigned int length, const unsigned int kTHvalue)
{
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }
unsigned int left = 0;
unsigned int left2 = 0;
unsigned int right = length - 1;
unsigned int right2 = length - 1;
unsigned int kthPlus = kTHvalue;
unsigned int kthMinus = kTHvalue;
while (right > left)
{
//Median of 3 optimization
if( array[right] < array[left] ) F_SWAP(array[left],array[right]);
if( array[right] < array[kTHvalue] ) F_SWAP(array[kTHvalue],array[right]);
if( array[kTHvalue] < array[left] ) F_SWAP(array[left],array[kTHvalue]);
//Median of 3 optimization
const float temp = array[kTHvalue];//this is pivot
kthPlus = kTHvalue + 1;
kthMinus = kTHvalue - 1;
while ((right2 > kTHvalue) && (left2 < kTHvalue))
{
do
{
left2++;
}while (array[left2] < temp);
do
{
right2--;
}while (array[right2] > temp);
F_SWAP(array[left2],array[right2]);
}
left2++;
right2--;
if (right2 < kthMinus)
{
while (array[left2] < temp)
{
left2++;
}
left = left2;
right2 = right;
kthMinus --;
}
else
if (kthPlus < left2)
{
while (array[right2] > temp)
{
right2--;
}
right = right2;
left2 = left;
kthPlus ++;
}
else
if( array[left] < array[right] )
{
F_SWAP(array[right],array[left]);
}
}
#undef F_SWAP
return array[kTHvalue];
}
//points to the center value (median)
int middleElement = elementsinarray / 2;
Dansby_kth(redArray, elementsinarray,middleElement);
Dansby_kth(greenArray, elementsinarray,middleElement);
Dansby_kth(blueArray, elementsinarray,middleElement);
But, it is in 64 bits and also uses qsort. I´m looking for some that does not need to use sorting algorithms inside.
but do you have it for 32 Bits so i can test it ?
float Dansby_kth(float array[], const unsigned int length, const unsigned int kTHvalue)
{
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }
unsigned int left = 0;
unsigned int left2 = 0;
unsigned int right = length - 1;
unsigned int right2 = length - 1;
unsigned int kthPlus = kTHvalue;
unsigned int kthMinus = kTHvalue;
while (right > left)
{
//Median of 3 optimization
if (array[right] < array[left]) F_SWAP(array[left], array[right]);
if (array[right] < array[kTHvalue]) F_SWAP(array[kTHvalue], array[right]);
if (array[kTHvalue] < array[left]) F_SWAP(array[left], array[kTHvalue]);
//Median of 3 optimization
const float temp = array[kTHvalue];//this is pivot
kthPlus = kTHvalue + 1;
kthMinus = kTHvalue - 1;
while ((right2 > kTHvalue) && (left2 < kTHvalue))
{
do
{
left2++;
} while (array[left2] < temp);
do
{
right2--;
} while (array[right2] > temp);
F_SWAP(array[left2], array[right2]);
}
left2++;
right2--;
if (right2 < kthMinus)
{
while (array[left2] < temp)
{
left2++;
}
left = left2;
right2 = right;
kthMinus--;
}
else
if (kthPlus < left2)
{
while (array[right2] > temp)
{
right2--;
}
right = right2;
left2 = left;
kthPlus++;
}
else
if (array[left] < array[right])
{
F_SWAP(array[right], array[left]);
}
}
#undef F_SWAP
return array[kTHvalue];
}
Ok, i´ll give a try and test the Heapmedian3 to see if it works as expected. However, as you said, it needs an extra buffer in memory which could be overkilling when we need to calculate the median of a larger array then only 27 bytes longs (say a array of floats, ints etc) with 100 Mb for example.
I[´ll try porting it to see if it works and also check for speed.
I was trying to port Dansby_kth that is very fast and don´t uses extra buffer, but it is producing a wrong result and i can´t fix that :(
Do you have any idea of what could be possibly wrong ?
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }
float Dansby_kth(float array[], const unsigned int length)
{
unsigned int kTHvalue = length / 2 - 1;
unsigned int left = 0;
unsigned int left2 = 0;
unsigned int right = length - 1;
unsigned int right2 = length - 1;
unsigned int kthPlus = kTHvalue;
unsigned int kthMinus = kTHvalue;
while (right > left)
{
if( array[right] < array[left] ) F_SWAP(array[left],array[right]);
if( array[right] < array[kTHvalue] ) F_SWAP(array[kTHvalue],array[right]);
if( array[kTHvalue] < array[left] ) F_SWAP(array[left],array[kTHvalue]);
const float temp = array[kTHvalue];//this is pivot
kthPlus = kTHvalue + 1;
kthMinus = kTHvalue - 1;
while ((right2 > kTHvalue) && (left2 < kTHvalue))
{
do
{
left2++;
}while (array[left2] < temp);
do
{
right2--;
}while (array[right2] > temp);
F_SWAP(array[left2],array[right2]);
}
left2++;
right2--;
if (right2 < kthMinus)
{
while (array[left2] < temp)
{
left2++;
}
left = left2;
right2 = right;
kthMinus --;
}
else
if (kthPlus < left2)
{
while (array[right2] > temp)
{
right2--;
}
right = right2;
left2 = left;
kthPlus ++;
}
else
if( array[left] < array[right] )
{
F_SWAP(array[right],array[left]);
}
}
return array[kTHvalue];
}
Ok, i´ll give a try and test the Heapmedian3 to see if it works as expected. However, as you said, it needs an extra buffer in memory which could be overkilling when we need to calculate the median of a larger array then only 27 bytes longs (say a array of floats, ints etc) with 100 Mb for example.
I[´ll try porting it to see if it works and also check for speed.
You may need to save the array for later use so it depends on the situation. In a benchmark test where the function is called repeatably that's at least the case.QuoteI was trying to port Dansby_kth that is very fast and don´t uses extra buffer, but it is producing a wrong result and i can´t fix that :(
Do you have any idea of what could be possibly wrong ?
It's difficult to say what the right answer is here. If the array is even (1,2,3,4) the answer is 2 or 3 so you need to make a choice there I guess. The returned result seems to be 18 and I assume the correct answer 17. In that case it overshoots the index by one. Reducing the index by one seems at least on this array to give the right answer.Code: [Select]#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }
float Dansby_kth(float array[], const unsigned int length)
{
unsigned int kTHvalue = length / 2 - 1;
unsigned int left = 0;
unsigned int left2 = 0;
unsigned int right = length - 1;
unsigned int right2 = length - 1;
unsigned int kthPlus = kTHvalue;
unsigned int kthMinus = kTHvalue;
while (right > left)
{
if( array[right] < array[left] ) F_SWAP(array[left],array[right]);
if( array[right] < array[kTHvalue] ) F_SWAP(array[kTHvalue],array[right]);
if( array[kTHvalue] < array[left] ) F_SWAP(array[left],array[kTHvalue]);
const float temp = array[kTHvalue];//this is pivot
kthPlus = kTHvalue + 1;
kthMinus = kTHvalue - 1;
while ((right2 > kTHvalue) && (left2 < kTHvalue))
{
do
{
left2++;
}while (array[left2] < temp);
do
{
right2--;
}while (array[right2] > temp);
F_SWAP(array[left2],array[right2]);
}
left2++;
right2--;
if (right2 < kthMinus)
{
while (array[left2] < temp)
{
left2++;
}
left = left2;
right2 = right;
kthMinus --;
}
else
if (kthPlus < left2)
{
while (array[right2] > temp)
{
right2--;
}
right = right2;
left2 = left;
kthPlus ++;
}
else
if( array[left] < array[right] )
{
F_SWAP(array[right],array[left]);
}
}
return array[kTHvalue];
}
EDIT: Spotted an error -> jb changed to jbe
New attachment: see reply#13
( and hope it's error free? ) :bgrin::biggrin: No, I'm afraid
.data?
MedianBuffer db 1024 dup (?) ; minimum buffer size must be as large as the biggest number in the data set
Must be:.data
MedianBuffer db 1024 dup (0) ; minimum buffer size must be as large as the biggest number in the data set
CalculateMedian proc uses ebx esi edi
mov esi,offset DataSet
mov edi,offset MedianBuffer
mov ecx,PopulationSize
PopulateLP:
mov eax, [esi]
inc byte ptr [edi+eax]
add esi,4
dec ecx
jnz PopulateLP
mov ecx,PopulationSize
mov ebx,ecx
shr ebx,1
and ecx,1 ; Check even or uneven
jnz UnEven
dec ebx ; Make even
UnEven:
xor edx,edx
mov eax,-1
FindMedianLP:
inc eax
movzx ecx,byte ptr[edi+eax]
add edx,ecx
cmp edx,ebx
jbe FindMedianLP
;----------------------
.if edx > ebx
mov ecx, edx
sub ecx, ebx
.if ecx > 0
jnz Done
.endif
.endif
;-----------------------
mov edx,PopulationSize
and edx,1 ; Check even or uneven
jnz Done
mov edx,eax
FindSecondMiddleNumberLP:
inc eax
movzx ecx,byte ptr[edi+eax]
test ecx,ecx
jz FindSecondMiddleNumberLP
add eax,edx
shr eax,1
Done: ; Result in eax
ret
CalculateMedian endp
cmp ecx,PopulationSize/2
jb average
align 4
CalculateMedian proc uses ebx esi edi
mov esi,offset DataSet
mov edi,offset MedianBuffer
mov ecx,PopulationSize
PopulateLP:
mov eax,[esi]
inc byte ptr [edi+eax]
add esi,4
dec ecx
jnz PopulateLP
mov ecx,PopulationSize
mov ebx,ecx
shr ebx,1
and ecx,1 ; Check even or uneven
jnz UnEven
dec ebx ; Make even
UnEven:
xor edx,edx
mov eax,-1
FindMedianLP:
inc eax
movzx ecx,byte ptr[edi+eax]
add edx,ecx
cmp edx,ebx
jbe FindMedianLP
mov edx,PopulationSize
and edx,1 ; Check even or uneven
jnz Done
mov edx,eax
FindSecondMiddleNumberLP:
inc eax
movzx ecx,byte ptr[edi+eax]
cmp ecx,PopulationSize/2
jb Average
test ecx,ecx
jz FindSecondMiddleNumberLP
Average:
add eax,edx
shr eax,1
Done: ; Result in eax
ret
(...)
FindMedianLP:
inc eax
movzx ecx,byte ptr[edi+eax]
add edx,ecx
cmp edx,ebx
jbe FindMedianLP
(....)
MedianBuffer dd 256 dup (0)
(...)
FindMedianLP:
inc eax
add ecx [edi+eax] ; or add ecx [MedianBuffer+eax] to avoid using more registers.
cmp edx,ebx
jbe FindMedianLP
[MedianPopulateTable: D$ 0 #256] ; minimum buffer size must be as large as the biggest number in the data set. In my case it is a range of pixels so, 256 values only. In this case, i´m using the table as a 256 dword (D$ in RosAsm syntax) array
Proc CalculateMedian3:
Arguments @Input, @ArraySize
Local @IsArraySizeOdd, @MiddleTablePos, @MedianLeftVal, @iCounter
Uses ecx
mov ecx D@ArraySize | mov D@iCounter ecx ; "ArraySize" is the same as PopulationSize. I just renamed it as array so i could better follow what it is doing in the other routines i´m testing it with.
mov ecx D@Input
mov D@IsArraySizeOdd &TRUE
; Populate Buffer Table
@PopulateLP:
mov eax D$ecx
inc D$MedianPopulateTable+eax*4
add ecx 4
dec D@iCounter | jne @PopulateLP
mov ecx D@ArraySize | mov D@MiddleTablePos ecx | shr D@MiddleTablePos 1
Test_If cl cl ; Is number even ? Set IsArraySizeOdd to FALSE and decrease the position of the array table to be compared with
dec D@MiddleTablePos
mov D@IsArraySizeOdd &FALSE ; set Flag to not odd number (So, it's a even)
Test_End
; FindMedianLP
xor ecx ecx
mov eax 0-1
Do
inc eax
add ecx D$MedianPopulateTable+eax*4
Loop_Until ecx > D@MiddleTablePos
If D@IsArraySizeOdd = &FALSE ; is Array Size even ?
; FindSecondMiddleNumberLP
mov D@MedianLeftVal eax
Do
inc eax
Loop_Until D$MedianPopulateTable+eax*4 <> 0
add eax D@MedianLeftVal
shr eax 1
End_If
; Result in eax
EndP
Now it is need to find the median of all pixels on all images and create another image containing only the medians values (for x and y pos)
I´m not sure i understood what is wrong with the code as suggested by HSE. I´m testing it here and didn´t found incorrect values so far. What i´m missing ?maybe you should have a big dynamic alloc memory buffer already cleared,so that simple move on to next 256 bytes,next time its called?,your going to use dynamic alloc anyway to process big image data
For a tiny table of pixels (256 integer values only from 0 to 255) the algorithm seems to work fine, so far. I ported it to RosAsm trying to make it a bit faster (not measured to see what happened), using local variables so i could try to understand what the algo is doing. . Also it only needs to preserve ecx, since the data is stored in local vars, so no extra push/pop chains at the start and end of the function. And i set the data in the table to dwords, to replace routines like this:Code: [Select](...)
FindMedianLP:
inc eax
movzx ecx,byte ptr[edi+eax]
add edx,ecx
cmp edx,ebx
jbe FindMedianLP
To something like:Code: [Select](....)
MedianBuffer dd 256 dup (0)
(...)
FindMedianLP:
inc eax
add ecx [edi+eax] ; or add ecx [MedianBuffer+eax] to avoid using more registers.
cmp edx,ebx
jbe FindMedianLP
Can you explain to me how does it is working ? I´m not fully understood why.
And also, does it needs to use a table ? This is because even if the algo can be improved for speed, the table should byte zeroed whenever the function is called and if we are dealing with huge amount of data like pixels in videos or in my case, hundreds of images to be analysed), then, cleaning the table could slow down the algo a little bit.
And another question,. how to make it work with SSE on Float or Real Numbers ? I mean, the algo is working on a array of integers (gray maps), but it won´t work on another table that is a array of floats to used to compute the median of Sobel data.
Maybe other assemblers behave different?:thumbsup: Could be that!
movzx ecx,byte ptr[edi+eax]
cmp ecx,PopulationSize/2 <--- ????
jb Average
.if edx > ebx
mov ecx, edx
sub ecx, ebx
cmp ecx , 0
jg Done
.endif
Now it is need to find the median of all pixels on all images and create another image containing only the medians values (for x and y pos)
- If we are talking about 32-bit colour images, does it mean you need the routine at byte size level?
- Have you thought of using the arithmetic medium instead? It's definitely much faster, and in practice there might be little difference.
Below an illustration of the difference between median and average. The only little problem is that your data are not sorted :cool:
[ImgDataHeader:
ImgDataHeader.Size: D$ 0 ; Total size of our structure with the arrays included
ImgDataHeader.Count: D$ 0 ; total amount of valid data structures. So, the total of images to be processed that fits to a inputed size. Ex: Count only the images on a directory that have resolution of 960*720
ImgDataHeader.DataSize: D$ 0 ; Size of each data chunk. Width*Height*4
ImgDataHeader.MaskWidth: D$ 0 ; Width of the image that will be later converted to a mask
ImgDataHeader.MaskHeight: D$ 0] ; Width of the image that will be later converted to a mask
Immediately followed by a array of N "FullImageData" structures that are only pointers to the processed data. N is given by "ImgDataHeader.Count"
[FullImageData:
FullImageData.pGrayMap: D$ 0 ; Pointer to graymap
FullImageData.pSobelX: D$ 0 ; Pointer to SobelXMap
FullImageData.pSobelY: D$ 0 ; Pointer to SobelYMap
FullImageData.pSobel: D$ 0] ; Pointer to SobelMap
(....)
FullImageData.pGrayMap100: D$ 0 ; Pointer to graymap of the 100th image. Ex: Points to GreymapImage1 (The 1st greymap)
FullImageData.pSobelX100: D$ 0 ; Pointer to SobelXMap of the 100th image Ex: Points to SobelXMap1 (The 1st SobelXMap1)
FullImageData.pSobelY100: D$ 0 ; Pointer to SobelYMap of the 100th image Ex: Points to SobelYMap1 (The 1st SobelYMap1)
FullImageData.pSobel100: D$ 0] ; Pointer to SobelMap of the 100th image Ex: Points to SobelMap1 (The 1st SobelMap1)
]
Immediatelly followed by the maps, Starting with the greymap
[GreymapImage1: D$ 0#(960*720*4)
(...)
GreymapImage100: D$ 0#(960*720*4)
(...)
[SobelXMap1: D$ 0#(960*720*4)
SobelYMap1: D$ 0#(960*720*4)
SobelMap1: D$ 0#(960*720*4)
(...)
SobelXMap100: D$ 0#(960*720*4)
SobelYMap100: D$ 0#(960*720*4)
SobelMap100: D$ 0#(960*720*4)
I´m not sure i understood what is wrong with the code as suggested by HSE. I´m testing it here and didn´t found incorrect values so far. What i´m missing ?maybe you should have a big dynamic alloc memory buffer already cleared,so that simple move on to next 256 bytes,next time its called?,your going to use dynamic alloc anyway to process big image data
For a tiny table of pixels (256 integer values only from 0 to 255) the algorithm seems to work fine, so far. I ported it to RosAsm trying to make it a bit faster (not measured to see what happened), using local variables so i could try to understand what the algo is doing. . Also it only needs to preserve ecx, since the data is stored in local vars, so no extra push/pop chains at the start and end of the function. And i set the data in the table to dwords, to replace routines like this:Code: [Select](...)
FindMedianLP:
inc eax
movzx ecx,byte ptr[edi+eax]
add edx,ecx
cmp edx,ebx
jbe FindMedianLP
To something like:Code: [Select](....)
MedianBuffer dd 256 dup (0)
(...)
FindMedianLP:
inc eax
add ecx [edi+eax] ; or add ecx [MedianBuffer+eax] to avoid using more registers.
cmp edx,ebx
jbe FindMedianLP
Can you explain to me how does it is working ? I´m not fully understood why.
And also, does it needs to use a table ? This is because even if the algo can be improved for speed, the table should byte zeroed whenever the function is called and if we are dealing with huge amount of data like pixels in videos or in my case, hundreds of images to be analysed), then, cleaning the table could slow down the algo a little bit.
And another question,. how to make it work with SSE on Float or Real Numbers ? I mean, the algo is working on a array of integers (gray maps), but it won´t work on another table that is a array of floats to used to compute the median of Sobel data.
thats what I mean with SSE MAXPS and MINPS is:
if arraya>arrayb swap(arraya,arrayb) :
arrayb=MAXSS(arraya,arrayb)
arraya=MINSS(arraya,arrayb)
so unroll it MAXPS and MINPS does 4 of these conditional swaps in parallel, which could help with big image data,and MAXSS and MINSS to sort the final few numbers
# W_m is the cropped watermark
W_m = poisson_reconstruct(cropped_gx, cropped_gy)
# Get the number of images in the folder
num_images = len(gxlist)
J, img_paths = get_cropped_images(
'images/fotolia_processed', num_images, start, end, cropped_gx.shape)
Hi NiDud. When the number is even you simply compute the average of the 2 located in the middle. Like this:
Array = (1,2,3,4)
You will calculate the median using 2 and 3
Sum them up and divide by 2
So. The media for this is: (2+3)/2 = 2.5
So, you are computing the width / 2 and whatever values are in position (width/2) and (Width/2)+1 you simply calculate the average of those 2.
instead using a table of bytes, i simply convert them to Dwords to gain speed in other parts of the routines, to avoid having to use movzx eax B$ecx+....... So a simple mov eax D$xxx is enough to eliminate the needs of using extra code.
I have forgotten to say: very smart idea!
Maybe other assemblers behave different?:thumbsup: Could be that!
CalculateMedian proc p:ptr, n:dword
local count[256]:dword
xor eax,eax
mov ecx,256
mov edx,edi
lea edi,count
rep stosd
mov edi,edx
mov edx,p
mov ecx,n
.repeat
movzx eax,byte ptr [edx]
inc count[eax*4]
add edx,1
.untilcxz
mov ecx,n
shr ecx,1
xor edx,edx
mov eax,-1
.repeat
inc eax
add edx,count[eax*4]
.until edx > ecx
ret
CalculateMedian endp
Let's make a test case: test1000.zip are 1,000 random DWORDs between 0 and 125,248. The graph shows their distribution when sorted.Nice.
CalculateMedian proc uses ebx esi edi
mov esi,offset DataSet
mov edi,offset MedianBuffer ; Zero the MedianBuffer before using this routine.
mov ecx,PopulationSize
PopulateLP:
mov eax,[esi]
inc byte ptr [edi+eax] ; Populate the Median buffer with the values as index numbers
add esi,4 ; and keep track for duplicates in one go.
dec ecx
jnz PopulateLP
xor esi,esi ; Uneven marker.
mov ecx,PopulationSize
mov ebx,ecx
shr ebx,1 ; The middle of the PopulationSize.
and ecx,1 ; Check even or uneven.
jnz UnEven
dec ebx ; Make even.
inc esi ; Set marker true = even.
UnEven:
xor edx,edx ; Median counter.
mov eax,-1 ; Number counter.
FindMedianLP:
inc eax
movzx ecx,byte ptr[edi+eax] ; Count of single and/or duplicate numbers. ( zero is not a number )
add edx,ecx ; Add it to the Median counter.
cmp edx,ebx ; Compare if we reached the middle of the PopulationSize.
jbe FindMedianLP ; Repeat until we reached the middle.
dec esi ; Check even or uneven.
js Ready ; Ready if uneven.
cmp ecx,1 ; Check if next number is a duplicate. ( no averaging needed )
ja Ready ; Ready if next number is the same value as the previous number.
mov edx,eax ; Save the first middle number to calculate the average of 2 middle numbers.
FindSecondMiddleNumberLP:
inc eax
movzx ecx,byte ptr[edi+eax]
test ecx,ecx ; If zero, its not a number.
jz FindSecondMiddleNumberLP
add eax,edx ; Second middle number found, add it to the first middle number.
shr eax,1 ; And get the average.
Ready:
ret
CalculateMedian endp
align 4
CalculateMedian proc uses ebx esi edi
mov esi,offset DataSet
mov edi,offset MedianBuffer ; Zero the MedianBuffer before using this routine.
mov ecx,PopulationSize
PopulateLP:
mov eax,[esi]
inc dword ptr[edi+eax*4] ; Populate the Median buffer with the values as index numbers
add esi,4 ; and keep track for duplicates in one go.
dec ecx
jnz PopulateLP
xor esi,esi ; Uneven marker.
mov ecx,PopulationSize
mov ebx,ecx
shr ebx,1 ; The middle of the PopulationSize.
and ecx,1 ; Check even or uneven.
jnz UnEven
dec ebx ; Make even.
inc esi ; Set marker true = even.
UnEven:
xor edx,edx ; Median counter.
mov eax,-1 ; Number counter.
FindMedianLP:
inc eax
mov ecx,[edi+eax*4] ; Count of single and/or duplicate numbers. ( zero is not a number )
add edx,ecx ; Add it to the Median counter.
cmp edx,ebx ; Compare if we reached the middle of the PopulationSize.
jbe FindMedianLP ; Repeat until we reached the middle.
dec esi ; Check even or uneven.
js Ready ; Ready if uneven.
cmp ecx,1 ; Check if next number is a duplicate. ( no averaging needed )
ja Ready ; Ready if next number is the same value as the previous number.
mov edx,eax ; Save the first middle number to calculate the average of 2 middle numbers.
FindSecondMiddleNumberLP:
inc eax
mov ecx,[edi+eax*4]
test ecx,ecx ; If zero, its not a number.
jz FindSecondMiddleNumberLP
add eax,edx ; Second middle number found, add it to the first middle number.
shr eax,1 ; And get the average.
Ready:
ret
CalculateMedian endp
You can not divide a number by 1000, because it is the index number.Yes, I can. Exactly because the index number is the number. I'm scaling down data, and loosing precision, nothing else. :biggrin:
It fails, because there are numbers higher than 1023 in the file.:arrow_right: It fails because fails. Maximum number is 125.
Then he needs only an index buffer of 256 * 1 dword = 1024 bytes and can handle 4294967295 duplicate colors.:thumbsup:
My Median routine is based on integer numbers ( used as indexed numbers ) and will not work with Real4 numbers. ( @HSE this is why a data set of Real4 values will fail. :thumbsup: )
Pixel data are integer data, that's why I came up with this Algorithm using indexed numbers.
mov ecx, edx
sub ecx, ebx
cmp ecx , 1
jg Done
Hi guga,
So, if I did understand it correctly:
1. Convert 100 images to gray values.
2. Get the median of 100 pixels. ( from the same x/y position of all those 100 images )
3. If the image size = 960*720 pixels you'll have at the end 691200 median values for your Sobel function.
Is the above correct?
I asked you this because you use pixel data which is in 4 bytes format per ARGB value.
Every gray pixel value has a range of 0-255 because each color component in the gray RGB value have the same value.
So you can store a gray value as a byte.
My Median routine is based on integer numbers ( used as indexed numbers ) and will not work with Real4 numbers. ( @HSE this is why a data set of Real4 values will fail. :thumbsup: )
Pixel data are integer data, that's why I came up with this Algorithm using indexed numbers.
But you want the Median values to be in Real4 format?Yes, please. If you can, could you try a version to work with Real4 ? On this way i could properly use it in step "5" (and others) that are using Real4 data rather then integers to detect the watermark. This is because i´m not sure if calculate the median previously biased only on the integer values of the grey will also works for the rest of the detect watermark algorithm. So i could try later both methods and see which one is getting closer to a result according to the google article.
You can convert and save the calculated Median integer value as a real4 value at the end of this Algorithm and use it for the Sobel function. ( and gain 0.5 value precision for double middle numbers )
And from that point you can do your Quartile Sobel Real4 calculations in fast SIMD instructions SSE SSE2 etc.
def estimate_watermark(foldername):
"""
Given a folder, estimate the watermark (grad(W) = median(grad(J))) follow Eq.4
Also, give the list of gradients, so that further processing can be done on it
"""
if not os.path.exists(foldername):
warnings.warn("Folder does not exist.", UserWarning)
return None
images = []
for r, dirs, files in os.walk(foldername):
# Get all the images
for file in files:
img = cv2.imread(os.sep.join([r, file]))
if img is not None:
images.append(img)
else:
print("%s not found." % (file))
# Compute gradients
print("Computing gradients.")
gradx = list(map(lambda x: cv2.Sobel(
x, cv2.CV_64F, 1, 0, ksize=KERNEL_SIZE), images))
grady = list(map(lambda x: cv2.Sobel(
x, cv2.CV_64F, 0, 1, ksize=KERNEL_SIZE), images))
# Compute median of grads
print("Computing median gradients.")
Wm_x = np.median(np.array(gradx), axis=0)
Wm_y = np.median(np.array(grady), axis=0)
return (Wm_x, Wm_y, gradx, grady)
gx, gy, gxlist, gylist = estimate_watermark('images/fotolia_processed')
def poisson_reconstruct(gradx, grady, kernel_size=KERNEL_SIZE, num_iters=11, h=0.1,
boundary_image=None, boundary_zero=True):
"""
Iterative algorithm for Poisson reconstruction.
Given the gradx and grady values, find laplacian, and solve for image
Also return the squared difference of every step.
h = convergence rate
"""
fxx = cv2.Sobel(gradx, cv2.CV_64F, 1, 0, ksize=kernel_size)
fyy = cv2.Sobel(grady, cv2.CV_64F, 0, 1, ksize=kernel_size)
laplacian = fxx + fyy
m, n, p = laplacian.shape
if boundary_zero == True:
est = np.zeros(laplacian.shape)
else:
assert(boundary_image is not None)
assert(boundary_image.shape == laplacian.shape)
est = boundary_image.copy()
est[1:-1, 1:-1, :] = np.random.random((m-2, n-2, p))
loss = []
for i in range(num_iters):
old_est = est.copy()
est[1:-1, 1:-1, :] = 0.25*(est[0:-2, 1:-1, :] + est[1:-1, 0:-2, :] +
est[2:, 1:-1, :] + est[1:-1, 2:, :] - h*h*laplacian[1:-1, 1:-1, :])
error = np.sum(np.square(est-old_est))
loss.append(error)
return (est)
gx, gy, gxlist, gylist = estimate_watermark('images/fotolia_processed') ; <---- Already did previously (or trying to finish it to see if it ok)
# est = poisson_reconstruct(gx, gy, np.zeros(gx.shape)[:,:,0])
cropped_gx, cropped_gy = crop_watermark(gx, gy) ; <---------- Is really necessary crop the watermark here ??? It can b do0ne previously since the very 1st graymap to we work only on the already cropped tables. (Assuming it is necessary to crop btw)
# random photo
img = cv2.imread('images/fotolia_processed/fotolia_137840668.jpg') ; <--- read the image from where the watermark should be removed
im, start, end = watermark_detector(img, cropped_gx, cropped_gy) ; < ---- try to detect the watermark on this image as well. (This routine is where it will use canny)
# # Save result of watermark detection
# plt.subplot(1, 2, 1)
# plt.imshow(img)
# plt.subplot(1, 2, 2)
# plt.imshow(im)
# plt.savefig('images/results/watermark_detect_result.png')
# We are done with watermark estimation
# W_m is the cropped watermark
W_m = poisson_reconstruct(cropped_gx, cropped_gy) ; <--- Use poisson only on the target image i suppose
(...)
def crop_watermark(gradx, grady, threshold=0.4, boundary_size=2):
"""
Crops the watermark by taking the edge map of magnitude of grad(W)
Assumes the gradx and grady to be in 3 channels
@param: threshold - gives the threshold param
@param: boundary_size - boundary around cropped image
"""
W_mod = np.sqrt(np.square(gradx) + np.square(grady))
# Map image values to [0, 1]
W_mod = PlotImage(W_mod)
# Threshold the image with threshold=0.4
W_gray = image_threshold(np.average(W_mod, axis=2), threshold=threshold)
x, y = np.where(W_gray == 1)
# Boundary of cropped image (contain watermark)
xm, xM = np.min(x) - boundary_size - 1, np.max(x) + boundary_size + 1
ym, yM = np.min(y) - boundary_size - 1, np.max(y) + boundary_size + 1
return gradx[xm:xM, ym:yM, :], grady[xm:xM, ym:yM, :]
# return gradx[xm:xM, ym:yM, 0], grady[xm:xM, ym:yM, 0]
I thought this is what you wanted to do:
Instead of X and Y convolution matrices you like to use X and Y median matrices to detect the edges in the image?
I do not understand the part, you talk about 100 different pictures.
I'm not really sure what to make of the watermark removal part, or is it only about edge detection?
Ah, there are more than 256 different gray colors..... need to convert them all. :biggrin:
Float_Two dq 2.0
Float_One_255 dq 0.00392156862745098 ; (1/255)
Float_SobleVarG dq 180.3122292025696 ; 255/sqrt(2). This is to speed up the calculation rather then using divide.
SobelGetGPLusEx proc near
FractionY = dword ptr -14h
FractionX = dword ptr -10h
Divisor = dword ptr -0Ch
DataCheck = dword ptr -8
pReturn = dword ptr -4
pMatrix = dword ptr 8
pOutSobelX = dword ptr 0Ch
pOutSobelY = dword ptr 10h
pOutSobel = dword ptr 14h
push ebp
mov ebp, esp
sub esp, 14h
push edi
push edx
push ecx
finit
mov edi, [ebp+pMatrix]
; To calculate Gx^2 later. Therefore Gx = M3+M9 + 2*(M6-M4) - (M7+M1)
fld dword ptr [edi+20] ; M6
fsub dword ptr [edi+12] ; M4
fmul Float_Two
fadd dword ptr [edi+8] ; M3
fadd dword ptr [edi+32] ; M9
fsub dword ptr [edi+24] ; M7
fsub dword ptr [edi+0] ; M1
lea ecx, [ebp+DataCheck]
fistp dword ptr [ecx]
cmp [ebp+DataCheck], 0
jnz short loc_4C7649
fldz
fstp [ebp+FractionX]
jmp short loc_4C76AB
; ---------------------------------------------------------------------------
loc_4C7649:
cmp [ebp+DataCheck], -255 ; Blacks. Edge too dark. Limit was extrapolated
jge short loc_4C7675
xor edx, edx
mov ecx, 255
mov eax, [ebp+DataCheck]
neg eax
div ecx
inc eax
imul eax, 255
mov [ebp+Divisor], eax
fild [ebp+DataCheck]
fidiv [ebp+Divisor]
fstp [ebp+FractionX]
jmp short loc_4C76AB
; ---------------------------------------------------------------------------
loc_4C7675:
cmp [ebp+DataCheck], 255 ; Whites. Edge too brigth. Limit was extrapolated
jle short loc_4C769F
xor edx, edx
mov ecx, 255
mov eax, [ebp+DataCheck]
div ecx
inc eax
imul eax, 255
mov [ebp+Divisor], eax
fild [ebp+DataCheck]
fidiv [ebp+Divisor]
fstp [ebp+FractionX]
jmp short loc_4C76AB
; ---------------------------------------------------------------------------
loc_4C769F:
fild [ebp+DataCheck]
fmul Float_One_255
fstp [ebp+FractionX]
loc_4C76AB:
mov eax, [ebp+pOutSobelX]
mov ecx, [ebp+FractionX]
mov [eax], ecx
; To calculate Gy^2 later. Therefore Gy = M7+M9 + 2*(M8-M2) - (M3+M1)
fld dword ptr [edi+28] ; M8
fsub dword ptr [edi+4] ; M2
fmul Float_Two
fadd dword ptr [edi+24] ; M7
fadd dword ptr [edi+32] ; M9
fsub dword ptr [edi+8] ; M3
fsub dword ptr [edi+0] ; M1
lea ecx, [ebp+DataCheck]
fistp dword ptr [ecx]
cmp [ebp+DataCheck], 0
jnz short loc_4C76DD
fldz
fstp [ebp+FractionY]
jmp short loc_4C773F
; ---------------------------------------------------------------------------
loc_4C76DD:
cmp [ebp+DataCheck], -255 ; Blacks. Edge too dark. Limit was extrapolated
jge short loc_4C7709
xor edx, edx
mov ecx, 255
mov eax, [ebp+DataCheck]
neg eax
div ecx
inc eax
imul eax, 255
mov [ebp+Divisor], eax
fild [ebp+DataCheck]
fidiv [ebp+Divisor]
fstp [ebp+FractionY]
jmp short loc_4C773F
; ---------------------------------------------------------------------------
loc_4C7709:
cmp [ebp+DataCheck], 255 ; Whites. Edge too brigth. Limit was extrapolated
jle short loc_4C7733
xor edx, edx
mov ecx, 255
mov eax, [ebp+DataCheck]
div ecx
inc eax
imul eax, 255
mov [ebp+Divisor], eax
fild [ebp+DataCheck]
fidiv [ebp+Divisor]
fstp [ebp+FractionY]
jmp short loc_4C773F
; ---------------------------------------------------------------------------
loc_4C7733:
; Soble = sqrt((255*FractionX)^2+(255*FractionY)^2)) = G = sqrt(Gx^2+Gy^2)
; since FractionX and FractionY can have a maximum of 1 and -1, therefore sobleMax = (255/sqrt(2)) * sqrt(FractionX^2+FractionY^2)
fild [ebp+DataCheck]
fmul Float_One_255
fstp [ebp+FractionY]
loc_4C773F:
mov eax, [ebp+pOutSobelY]
mov ecx, [ebp+FractionY]
mov [eax], ecx
fld [ebp+FractionX]
fmul st, st
fld [ebp+FractionY]
fmul st, st
faddp st(1), st
fsqrt
fmul Float_SobleVarG
lea edx, [ebp+pReturn]
fistp dword ptr [edx]
mov eax, [ebp+pReturn]
cmp eax, 255
jbe short loc_4C776F
mov eax, 255
loc_4C776F:
mov ecx, [ebp+pOutSobel]
mov [ecx], eax
pop ecx
pop edx
pop edi
mov esp, ebp
pop ebp
retn 10h
SobelGetGPLusEx endp
pMatrix(In):
Pointer to a 3*3 Matrix of Dwords containing the gray values
pOutSobelX(Out):
Pointer to a variable that will store the SobelX value. The variable must be is Dword/Float Sized (It exports a float, but, they are the same size ;) ). The saved value is a Real4 from -1 to 1. So they are, actually, the normalized grey. i.e: Color/255
If the value is negative, it means the pixel is dark. If the value is positive, it means, the pixel is brighter. So, at the end, basically, it goes from 0 to 127 (means dark) and 128 to 255 meaning bright.
The sign here (negative or positive) was kept because we can use the values of Gx and Gy to compute the direction of the pixel with atan(Gy/Gx)
pOutSobelY(Out):
Pointer to a variable that will store the SobelY value. The variable must be is Dword/Float Sized (It exports a float, but, they are the same size ;) ). The saved value is a Real4 from -1 to 1. So they are, actually, the normalized grey. i.e: Color/255
If the value is negative, it means the pixel is dark. If the value is positive, it means, the pixel is brighter. So, at the end, basically, it goes from 0 to 127 (means dark) and 128 to 255 meaning bright.
The sign here (negative or positive) was kept because we can use the values of Gx and Gy to compute the direction of the pixel with atan(Gy/Gx)
pOutSobel(Out):
Pointer to a variable that will store the Sobel value. The variable must be a Dword. The saved value is a integer from 0 to 255.
TmpMatrix dd 212, 10, 15
dd 105, 5, 18
dd 215, 178, 0
lea eax D$ebp+SobelX
lea ebx D$ebp+SobelY
lea ecx D$ebp+Sobel
invoke SobelGetGPlusEx, offset TmpMatrix, ecx, ebx, eax
Proc SobelGetGPLusEx:
Arguments @pMatrix, @pOutSobelX, @pOutSobelY, @pOutSobel
Local @pReturn, @DataCheck, @Divisor, @FractionX, @FractionY
Uses edi, edx, ecx
finit
mov edi D@pMatrix
; To calculate Gx^2 later. Therfore Gx = M3+M9 + 2*(M6-M4) - (M7+M1)
fld F$edi+FloatMatricesInt.M6Dis | fsub F$edi+FloatMatricesInt.M4Dis | fmul R$Float_Two
fadd F$edi+FloatMatricesInt.M3Dis | fadd F$edi+FloatMatricesInt.M9Dis
fsub F$edi+FloatMatricesInt.M7Dis | fsub F$edi+FloatMatricesInt.M1Dis
lea ecx D@DataCheck | fistp F$ecx
If D@DataCheck = 0
fldz | fstp F@FractionX
Else_If D@DataCheck <s 0-255 ; Blacks. Edge too dark
xor edx edx | mov ecx 255 | mov eax D@DataCheck | neg eax | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionX
Else_If D@DataCheck >s 255 ; Whites. Edge too brigth
xor edx edx | mov ecx 255 | mov eax D@DataCheck | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionX
Else
fild D@DataCheck | fmul R$Float_One_255 | fstp F@FractionX
End_If
mov eax D@pOutSobelX | mov ecx D@FractionX | mov D$eax ecx
; To calculate Gy^2 later. Therefore Gy = M7+M9 + 2*(M8-M2) - (M3+M1)
fld F$edi+FloatMatricesInt.M8Dis | fsub F$edi+FloatMatricesInt.M2Dis | fmul R$Float_Two
fadd F$edi+FloatMatricesInt.M7Dis | fadd F$edi+FloatMatricesInt.M9Dis
fsub F$edi+FloatMatricesInt.M3Dis | fsub F$edi+FloatMatricesInt.M1Dis
lea ecx D@DataCheck | fistp F$ecx
If D@DataCheck = 0
fldz | fstp F@FractionY
Else_If D@DataCheck <s 0-255 ; Blacks
xor edx edx | mov ecx 255 | mov eax D@DataCheck | neg eax | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionY; | fmul R$Float255
Else_If D@DataCheck >s 255 ; Whites
xor edx edx | mov ecx 255 | mov eax D@DataCheck | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionY; | fmul R$Float255
Else
fild D@DataCheck | fmul R$Float_One_255 | fstp F@FractionY
End_If
mov eax D@pOutSobelY | mov ecx D@FractionY | mov D$eax ecx
; Soble = sqrt((255*FractionX)^2+(255*FractionY)^2)) = G = sqrt(Gx^2+Gy^2)
; since FractionX and FractionY can have a maximum of 1 and -1, therefore sobleMax = (255/sqrt(2)) * sqrt(FractionX^2+FractionY^2)
fld F@FractionX | fmul ST0 ST0 | fld F@FractionY | fmul ST0 ST0 | faddp ST1 ST0 | fsqrt | fmul R$Float_SobleVarG
lea edx D@pReturn
fistp F$edx
mov eax D@pReturn
If eax > 255 ; Can be replaced with movzx eax ax. But i´m keping the comparition to check for errors while debugging.
mov eax 255
End_If
mov ecx D@pOutSobel | mov D$ecx eax
EndP
; variables
[Float_Two: R$ 2]
[FloatOne_255: R$ (1/255)]
[Float_SobleVarG: R$ 180.3122292025696187222153123367365]; 255/sqrt(2))
; equates
[FloatMatricesInt.M1Dis 0
FloatMatricesInt.M2Dis 4
FloatMatricesInt.M3Dis 8
FloatMatricesInt.M4Dis 12
FloatMatricesInt.M5Dis 16
FloatMatricesInt.M6Dis 20
FloatMatricesInt.M7Dis 24
FloatMatricesInt.M8Dis 28
FloatMatricesInt.M9Dis 32]
[Size_Of_FloatMatricesInt 36]
998 µs for MbMedian M=124998
233 µs for CalcMedian M=124998
1034 µs for ArraySort M=124998
#elements=10000, maxY=501499
1662 µs for MbMedian M=125749
92 µs for CalcMedian M=110213 <<<<<<<<<<<<<
1567 µs for ArraySort M=125749
#elements=20000, maxY=501499
3334 µs for MbMedian M=126000
138 µs for CalcMedian M=126000
2549 µs for ArraySort M=126000
#elements=40000, maxY=501499
Hi Jochen,
Don't know... could be a bug...
It is not so fast when there are great gaps between the numbers.
But this routine is special designed for gray color bitmap pixel data.
- they are normally close to each other and there are many duplicate values and have a small range 0-255
- Big chance the whole buffer is populated at all offsets.
- it can process many thousands of pixel values very fast.
Hi Jochen,
Don't know... could be a bug...
It's a very nice algo, Marinus, similar to a bucket sort. And I found the culprit: I had zeroed the MedianBuffer using the number of elements instead of the sizeof the buffer. Now it works!
For Guga's application, with 0...255 as value range, it's fine. For other data sets with different ranges, especially with negative numbers, it will be difficult.
Just some ideas:
For example, to get rid of the errors at the image borders:
- image = 100 * 100 pixels
- X and Y matrices = 3 * 3
Allocate a work buffer of 102 * 102 pixels with zeros.
Copy the bitmap data exactly in the middle so we have created a 1 pixel space around the image.
Now we can use the 3 * 3 matrix convolution maps without restoring the image border pixels.
The 1 pixel extra boundary space is only read, not written.
To prevent overshoots in Gx and Gy, we could normalize the matrix members.
If I'm correct, we will have clamped G values between 0-255 after sqrt(Gx^2+Gy^2).
This way we can also experiment with other matrix sizes and/or for example Gaussian ( evenly distributed ) convolution maps.
This is all theoretical, of course, and it has to be proven to work in the real world. :biggrin:
It is not so fast when there are great gaps between the numbers.
But this routine is special designed for gray color bitmap pixel data.
- they are normally close to each other and there are many duplicate values and have a small range 0-255
- Big chance the whole buffer is populated at all offsets.
- it can process many thousands of pixel values very fast.
Think it's the other way around. The sort algo only have to loop half the array if all are equal. Yours have to loop the whole array for the count and (depending on the value) another half to calculate so that's 1:3 in favor of the sort algo.
500 dup(110)
total [0 .. 3], 1++
282577 cycles 0.asm: median
619579 cycles 1.asm: CalculateMedian
The sort algo gets slower as more swaps are needed and yours a fixed overhang by clearing the count buffer.
The overhang is only clearing a buffer of 256 items put that against thousands of pixels.
The median could be in the last part of the image.
For example:
9,10,11,12,13,1,2,3,4,5,6,8 -> median is 7 ( the average of 6 and 8 )
Hi guga,
COOL! :cool: :cool: :cool:
We don't need other graphics libraries.
GdiPlus is our friend.
We can load images, manipulate the bitmap data and save the new images with it.
The overhang is only clearing a buffer of 256 items put that against thousands of pixels.
I'm not sure how it's going to be used but I assumed the result was the median value of one specific pixel from each bitmap. If that's the case the input is 100 bytes.
Also, the Real4 version would also be needed, for what i saw in the others routines
Hi Siekmanski
I found an issue with the median calculation. When analyzing the routine in the python version, i saw that it also takes the median for negative values as well.
The values on SobelX and SobelY can result in either positive, negative or zeroes. How to make your algo also works for signed integers ? Ex: say we have a sequence of: -100, -20, -5, 0, 36, 158, -12, 14, 66 etc, etc
Or this: -12, 0, 145 . Median = 0. (limited, therefore from -255 to 255)
-100, -20, -5, 0, 36, 158, -12, 14, 66 => Median = 0
-20, -5, 0, 36, 158, -12, 14, 66 => Median = 7
Also, the Real4 version would also be needed, for what i saw in the others routines
def crop_watermark(gradx, grady, threshold=0.4, boundary_size=2):
"""
Crops the watermark by taking the edge map of magnitude of grad(W)
Assumes the gradx and grady to be in 3 channels
@param: threshold - gives the threshold param
@param: boundary_size - boundary around cropped image
"""
W_mod = np.sqrt(np.square(gradx) + np.square(grady)) <---- This is, actually the G. Always positive.
# Map image values to [0, 1]
W_mod = PlotImage(W_mod) ; What is this ?
# Threshold the image with threshold=0.4
W_gray = image_threshold(np.average(W_mod, axis=2), threshold=threshold)
x, y = np.where(W_gray == 1)
# Boundary of cropped image (contain watermark)
xm, xM = np.min(x) - boundary_size - 1, np.max(x) + boundary_size + 1
ym, yM = np.min(y) - boundary_size - 1, np.max(y) + boundary_size + 1
return gradx[xm:xM, ym:yM, :], grady[xm:xM, ym:yM, :] ; < ????????
# return gradx[xm:xM, ym:yM, 0], grady[xm:xM, ym:yM, 0]
def image_threshold(image, threshold=0.5):
'''
Threshold the image to make all its elements greater than threshold*MAX = 1
'''
m, M = np.min(image), np.max(image)
im = PlotImage(image)
im[im >= threshold] = 1
im[im < 1] = 0
return im
Also, the Real4 version would also be needed, for what i saw in the others routines
You can forget the Real4 version, Marinus' algo is based on the number itself being an index. That works fine for a byte (0...255), but a Real4 has a range of 0...FFFFFFFFh, so you would need an index buffer of 4GB. Use a sort algo instead.
3571 µs for converting from REAL4 to int32
648 ms for MbMedian M=249925309
720 ms for ArraySort M=249925309
#elements=5120000, maxY=767 Mio
xor ecx, ecx
.Repeat
fld REAL4 ptr [esi]
fmul FP4(10.0)
fistp DWORD ptr [edi]
inc ecx
.Until ecx>=eax
As far as I understand so far, correct me if I'm wrong....
- Perform edge detection on all 100 images with the Sobel formula.
- Find the medians, collecting 100 values at all the same x,y locations from 100 edge detected images.
- Use those median values to construct a new image.
- Detect the location of the watermark.
- Perform some magic, to remove the watermark?
:biggrin:
Do you have the links to the Google articles?
def estimate_watermark(foldername):
(...)
# Compute median of grads
print("Computing median gradients.")
Wm_x = np.median(np.array(gradx), axis=0)
Wm_y = np.median(np.array(grady), axis=0)
return (Wm_x, Wm_y, gradx, grady)
def crop_watermark(gradx, grady, threshold=0.4, boundary_size=2):
"""
Crops the watermark by taking the edge map of magnitude of grad(W)
Assumes the gradx and grady to be in 3 channels
@param: threshold - gives the threshold param
@param: boundary_size - boundary around cropped image
"""
W_mod = np.sqrt(np.square(gradx) + np.square(grady)) <-------
I found it where VisualStudio generated negative values. It was the way the python version was created. It does not calculate the media of G directly, it calculate the median of Gx and later the median of Gy and only after it, it compute the value of G = sqrt (MedianGx^2 + MedianGy^2)
To work with range -256 to 256 you can add 256 to values then you range became 0 to 512. Later substract 256 to result.
CalculateMedian proc uses edi p:ptr, n:dword
local count[256]:byte
xor eax,eax
mov ecx,256
lea edi,count
rep stosb
mov edx,p
mov ecx,n
lea edi,count[128]
.repeat
movsx eax,byte ptr [edx]
inc byte ptr [edi+eax]
add edx,1
.untilcxz
mov ecx,n
shr ecx,1
xor edx,edx
mov eax,-129
.repeat
inc eax
movzx edi,count[eax+128]
add edx,edi
.until edx > ecx
ret
CalculateMedian endp
If you want to make integers from real4 in the range 0-1, you have to multiply previously by 100, 1000 or precision you requiere. Later divide result for that number.
.686
.xmm
.model flat, stdcall
compare macro a, b
movss xmm0,[a]
comiss xmm0,[b]
exitm<>
endm
swapif macro a, b
ifdef _MAXSS
movss xmm0,[a]
movss xmm1,[b]
minss xmm1,[a]
maxss xmm0,[b]
movss [b],xmm0
movss [a],xmm1
else
compare(a, b)
.ifa
movss xmm1,[b]
movss [b],xmm0
movss [a],xmm1
.endif
endif
exitm<>
endm
swap macro a, b
movss xmm0,[a]
movss xmm1,[b]
movss [b],xmm0
movss [a],xmm1
exitm<>
endm
.code
option loopalign:4
median proc uses esi edi ebx p:ptr, n:dword
local level:dword
mov eax,n
.if eax > 1
mov esi,p
lea edi,[esi+eax*4-4]
mov level,0
.while 1
lea eax,[edi+4]
sub eax,esi
shr eax,3
lea ebx,[esi+eax*4]
swapif(esi, ebx)
swapif(esi, edi)
swapif(ebx, edi)
mov ecx,esi
mov edx,edi
.while 1
add ecx,4
.if ecx < edi
compare(ecx, ebx)
.continue .ifna
.endif
.while 1
sub edx,4
.break .if edx <= ebx
compare(edx, ebx)
.break .ifna
.endw
.break .if edx < ecx
swap(edx, ecx)
.if ebx == edx
mov ebx,ecx
.endif
.endw
add edx,4
.while 1
sub edx,4
.break .if edx <= esi
compare(edx, ebx)
.break .ifnz
.endw
mov eax,edx
mov edx,ecx
mov ecx,edi
sub ecx,edx
sub eax,esi
.if eax < ecx
.if edx < edi
push edx
push edi
inc level
.endif
.if esi < ecx
mov edi,ecx
.continue
.endif
.else
add eax,esi
.if esi < eax
push esi
push eax
inc level
.endif
.if edx < edi
mov esi,edx
.continue
.endif
.endif
.break .if !level
dec level
pop edi
pop esi
.endw
mov ecx,p
mov edx,n
shr edx,1
movss xmm0,[ecx+edx*4]
.ifnc
addss xmm0,[ecx+edx*4-4]
mov eax,2.0
movd xmm1,eax
divss xmm0,xmm1
.endif
.endif
ret
median endp
end
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : error A2070:invalid instruction operands
compare(1): Macro Called From
swapif(9): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : error A2070:invalid instruction operands
compare(2): Macro Called From
swapif(9): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : error A2008:syntax error : .
swapif(10): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : error A2070:invalid instruction operands
swapif(11): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : error A2070:invalid instruction operands
swapif(12): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : error A2070:invalid instruction operands
swapif(13): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
C:\masm32\examples\NidudMedian\3dframes2.asm(135) : fatal error A1011:directive must be in control block
swapif(14): Macro Called From
C:\masm32\examples\NidudMedian\3dframes2.asm(135): Main Line Code
align 4
CalculateMedianFP6 proc
; Fill the first 6 xmm registers with the unsorted 6 values
movaps xmm0,oword ptr [GrayPixels]
movlps xmm4,qword ptr [GrayPixels+16]
movaps xmm1,xmm0
movaps xmm2,xmm0
shufps xmm1,xmm1,00111001b
shufps xmm2,xmm2,01001110b
movaps xmm3,xmm0
movaps xmm5,xmm4
shufps xmm3,xmm3,10010011b
shufps xmm5,xmm5,00111001b
; Branchless floating point sort algorithm, only using registers.
; 6 numbers = max 12 swaps
movss xmm6,xmm1
movss xmm7,xmm4
minps xmm1,xmm2
minps xmm4,xmm5
maxps xmm2,xmm6
maxps xmm5,xmm7
movss xmm6,xmm0
movss xmm7,xmm3
minps xmm0,xmm2
minps xmm3,xmm5
maxps xmm2,xmm6
maxps xmm5,xmm7
movss xmm6,xmm0
movss xmm7,xmm3
minps xmm0,xmm1
minps xmm3,xmm4
maxps xmm1,xmm6
maxps xmm4,xmm7
movss xmm6,xmm1
movss xmm7,xmm0
minps xmm1,xmm4
minps xmm0,xmm3
maxps xmm4,xmm6
maxps xmm3,xmm7
movss xmm6,xmm2
movss xmm7,xmm1
minps xmm2,xmm5
minps xmm1,xmm3
maxps xmm5,xmm6
maxps xmm3,xmm7
movss xmm6,xmm2
minps xmm2,xmm4
maxps xmm4,xmm6
movss xmm6,xmm2
minps xmm2,xmm3
maxps xmm3,xmm6
addss xmm2,xmm3
mulss xmm2,real4 ptr [Half] ; Median in xmm2
ret
CalculateMedianFP6 endp
Hi guga,
A branchless floating point sort algorithm, only using registers.
To calculate the Median of the 6 numbers from a Sobel Matrix [ 3 * 3 ]Code: [Select]align 4
CalculateMedianFP6 proc
; Fill the first 6 xmm registers with the unsorted 6 values
movaps xmm0,oword ptr [GrayPixels]
movlps xmm4,qword ptr [GrayPixels+16]
movaps xmm1,xmm0
movaps xmm2,xmm0
shufps xmm1,xmm1,00111001b
shufps xmm2,xmm2,01001110b
movaps xmm3,xmm0
movaps xmm5,xmm4
shufps xmm3,xmm3,10010011b
shufps xmm5,xmm5,00111001b
; Branchless floating point sort algorithm, only using registers.
; 6 numbers = max 12 swaps
movss xmm6,xmm1
movss xmm7,xmm4
minps xmm1,xmm2
minps xmm4,xmm5
maxps xmm2,xmm6
maxps xmm5,xmm7
movss xmm6,xmm0
movss xmm7,xmm3
minps xmm0,xmm2
minps xmm3,xmm5
maxps xmm2,xmm6
maxps xmm5,xmm7
movss xmm6,xmm0
movss xmm7,xmm3
minps xmm0,xmm1
minps xmm3,xmm4
maxps xmm1,xmm6
maxps xmm4,xmm7
movss xmm6,xmm1
movss xmm7,xmm0
minps xmm1,xmm4
minps xmm0,xmm3
maxps xmm4,xmm6
maxps xmm3,xmm7
movss xmm6,xmm2
movss xmm7,xmm1
minps xmm2,xmm5
minps xmm1,xmm3
maxps xmm5,xmm6
maxps xmm3,xmm7
movss xmm6,xmm2
minps xmm2,xmm4
maxps xmm4,xmm6
movss xmm6,xmm2
minps xmm2,xmm3
maxps xmm3,xmm6
addss xmm2,xmm3
mulss xmm2,real4 ptr [Half] ; Median in xmm2
ret
CalculateMedianFP6 endp
I suceeded to dl your AsmC to make it run, but don´t know how to generate the executable file so i can see what are those macros used in 'median' function all about. Can you make a executable file for the SSE/Float version one (32 bits, pls), so i can test ?
;
; build: asmc -pe medianf.asm
;
.686
.xmm
.model flat, stdcall
compare macro a, b
movss xmm0,[a]
comiss xmm0,[b]
exitm<>
endm
swapif macro a, b
ifdef _MAXSS
movss xmm0,[a]
movss xmm1,[b]
minss xmm1,[a]
maxss xmm0,[b]
movss [b],xmm0
movss [a],xmm1
else
compare(a, b)
.ifa
movss xmm1,[b]
movss [b],xmm0
movss [a],xmm1
.endif
endif
exitm<>
endm
swap macro a, b
movss xmm0,[a]
movss xmm1,[b]
movss [b],xmm0
movss [a],xmm1
exitm<>
endm
.data
array dd -1.0, 18.0, 0.0, 0.0, 0.0, 0.0, 2.0, 2.0,
5.0, 99.0, 600.0, 785.0, 758.0, 17.0, 11.0, 21.0,
23.0, 24.0, 66.0, 21.0, 14.0, 57.0, 17.0
.code
;option loopalign:4
median proc uses esi edi ebx p:ptr, n:dword
local level:dword
mov eax,n
.if eax > 1
mov esi,p
lea edi,[esi+eax*4-4]
mov level,0
.while 1
lea eax,[edi+4]
sub eax,esi
shr eax,3
lea ebx,[esi+eax*4]
swapif(esi, ebx)
swapif(esi, edi)
swapif(ebx, edi)
mov ecx,esi
mov edx,edi
.while 1
add ecx,4
.if ecx < edi
compare(ecx, ebx)
.continue .ifna
.endif
.while 1
sub edx,4
.break .if edx <= ebx
compare(edx, ebx)
.break .ifna
.endw
.break .if edx < ecx
swap(edx, ecx)
.if ebx == edx
mov ebx,ecx
.endif
.endw
add edx,4
.while 1
sub edx,4
.break .if edx <= esi
compare(edx, ebx)
.break .ifnz
.endw
mov eax,edx
mov edx,ecx
mov ecx,edi
sub ecx,edx
sub eax,esi
.if eax < ecx
.if edx < edi
push edx
push edi
inc level
.endif
.if esi < ecx
mov edi,ecx
.continue
.endif
.else
add eax,esi
.if esi < eax
push esi
push eax
inc level
.endif
.if edx < edi
mov esi,edx
.continue
.endif
.endif
.break .if !level
dec level
pop edi
pop esi
.endw
mov ecx,p
mov edx,n
shr edx,1
movss xmm0,[ecx+edx*4]
.ifnc
addss xmm0,[ecx+edx*4-4]
mov eax,2.0
movd xmm1,eax
divss xmm0,xmm1
.endif
.endif
ret
median endp
main proc
median( &array, 23 )
ret
main endp
end main
Hi Nidud
Can you please make a small test ?
Check the median (The integer version: CalculateMedian) for these, pls ? -216,-151,136
The correct value is -151.
I´m not sure if i ported correctly your version, because the result i´ve got is wrong
To calculate the edges on the X axis
[ X X X ]
Gx = [ 0 0 0 ]
[ X X X ]
To calculate the edges on the Y axis
[ X 0 X ]
Gy = [ X 0 X ]
[ X 0 X ]
Hi guga,
-> I found it where VisualStudio generated negative values. It was the way the python version was created. It does not calculate the media of G directly, it calculate the median of Gx and later the median of Gy and only after it, it compute the value of G = sqrt (MedianGx^2 + MedianGy^2)
This comment from you is the reason I wrote a 6 numbers sorting routine.
The Median Matrix holds only positive numbers ( pixel values are never negative )
Therefore the calculation for the average middle numbers is correct.
And don't forget you only have to sort 6 values per matrix!!!Code: [Select]To calculate the edges on the X axes
[ X X X ]
Gx = [ 0 0 0 ]
[ X X X ]
To calculate the edges on the Y axes
[ X 0 X ]
Gy = [ X 0 X ]
[ X 0 X ]
You only need these 6 numbers. ( X's in the matrices )
For the Sobel Matrix convolution as for the Median Matrix calculations.
( You can save the Gray Color pixel conversion calculations as floating point )
( So you need no conversions from integer and floating point between all the calculation )
For the Median Matrix calculation you only read the 6 [X] pixel values ( always positive ) sort them, get the Median and that's the Gx and Gy.
No convolution needed as with the Sobel coefficients.
The G values should always be positive. ( they represent the new pixel color )
gradx = list(map(lambda x: cv2.Sobel(
x, cv2.CV_64F, 1, 0, ksize=KERNEL_SIZE), images))
grady = list(map(lambda x: cv2.Sobel(
x, cv2.CV_64F, 0, 1, ksize=KERNEL_SIZE), images))
# Compute median of grads
print("Computing median gradients.")
Wm_x = np.median(np.array(gradx), axis=0)
Wm_y = np.median(np.array(grady), axis=0)
W_mod = np.sqrt(np.square(gradx) + np.square(grady))
Just step for step I think, before we lose track.That image processing is that similar to search and replace black pixels surrounding lighter pixels with grey ?
This weekend I have some more time.
If you like it ( otherwise you might think I will hijack your project, of course I don't want that ), I can code the image loader, save some bytes in memory etc.
Show the original image on screen and next to it show the different edge detection images routines.
Later we can figure out ( if we succeed... ) how to do the watermark stuff.
Haven't looked into that part yet.
Just step for step I think, before we lose track.
This weekend I have some more time.
If you like it ( otherwise you might think I will hijack your project, of course I don't want that ), I can code the image loader, save some bytes in memory etc.
Show the original image on screen and next to it show the different edge detection images routines.
Later we can figure out ( if we succeed... ) how to do the watermark stuff.
Haven't looked into that part yet.
That image processing is that similar to search and replace black pixels surrounding lighter pixels with grey ?
postprocessed antialias, compared to line/oval /text outputted with antialias for example gdi+
Cool :cool:, we are a team now. :thumbsup:Great ! :thumbsup: :thumbsup: :thumbsup: :thumbsup: :thumbsup:
The goal is to get the location of the watermark with the most suitable edge detection technique.
In a movie the watermark is static and always the same color while the background changes between dark to light.
So when taking the median of multiple frames, the watermark will show up as a brighter color as the surrounding background, I think?
My idea is:
- load an image.
- convert the colors to gray colors to memory in real4 format.
- use real4 format for the sobel or median sobel.
and save the real4 G results ( are always positive ) as bytes in memory ( for 100 images )
- then calculate the 100 pixels median of each x/y position in the 100 images.
Guga,
I doubt there is an exhaustively good way to remove a watermark, just for practice I have done it with a few photos and the best after tweaking what you can is done by cloning areas close to where the water mark is. While the human eye is tolerant on many things, field depth is not one of them and cloning from about the same field depth keeps at least that aspect looking OK. On an image it is basically texture that your eye will pick if you get it wrong.
If you can sample the colour(s) used in the water mark you can try reducing that specific set of colours from the masked areas of the water mark. I have recently got a piece of Russian software that appears to be designed for altering images of people (mainly pretty girls) but it can do some useful things, it seems to do some form of AI cloning the close area and filling in the gaps. Its not perfect but it does look OK if you have a bit of practice.
Congrats guga,enough Iterations you can sell it as new Godzilla :greenclp::greensml: :greensml: :greensml: :greensml:
can't it be fixed with some filter out single alone pixels?
; Equates used as flags to the type of fix we want onto the sobel operator.
[SOBEL_FIX_TRUNCATE 0] ; Normal result. Only truncated all above the limits of -255 and 255
[SOBEL_FIX_SIMPLE 1] ; common adjust of the edges. Much better.
; Equates used on the matrix
[FloatMatricesInt.M1Dis 0
FloatMatricesInt.M2Dis 4
FloatMatricesInt.M3Dis 8
FloatMatricesInt.M4Dis 12
FloatMatricesInt.M5Dis 16
FloatMatricesInt.M6Dis 20
FloatMatricesInt.M7Dis 24
FloatMatricesInt.M8Dis 28
FloatMatricesInt.M9Dis 32]
[Size_Of_FloatMatricesInt 36]
[Float_SobleVarG: R$ 180.3122292025696187222153123367365]; 255/sqrt(2)) . This is a variable used to speed up the math when computing the G
;;
SobelGetGPLusEx
Calculates the Sobel operator for a given amount of pixels stored opn a 3x3 matrix
Parameters:
pMatrix(In): The Inputed pixels stored on a 3x3 matrix where the Sobel edges are calculated. The input values must be in Real4 format from 0 to 255. (No integer yet. I´ll adapt to handle integers later)
To make easier understand, the matrix is written using equates and the values must be previously calculated and filled.
The format of the matrix is:
M1 M2 M3
M4 M5 M6
M7 M8 M9
So, at Point M1, Pixel at pos X0, Y0
Point M2, Pixel at pos X1, Y0
Point M3, Pixel at pos X2, Y0
Point M4, Pixel at pos X0, Y1
Point M5, Pixel at pos X1, Y1
Point M6, Pixel at pos X2, Y1
Point M7, Pixel at pos X0, Y2
Point M8, Pixel at pos X1, Y2
Point M9, Pixel at pos X2, Y2
pOutSobelX(out): Pointer to a variable that will Store the SobelX value (Magnitude in X direction. I.e.: Gx).
The stored value can be negative, positive or zero. All in Real4 format from -1 to 1
pOutSobelY(out): Pointer to a variable that will Store the SobelY value (Magnitude in Y direction. I.e.: Gy).
The stored value can be negative, positive or zero. All in Real4 format from -1 to 1
pOutSobel(Out): Pointer to a variable that will Store the Sobel value (Total Magnitude of the pixel. I.e.: G).
The output is a positive Integer value from 0 to 255. The magnitude of Sobel operator is given by:
G = sqrt(Gx^2+Gy^2)
FixType(In): A flag able to fix the Sobel operator when it exceeds the limits of -1 to 1 (Normalized)
It´s common to the Sobel Operator extrapolates the resultant limits (-255 to 255. Or -1 to 1, normalized)
We can see values (unormalized) in the order of -600, 482, 265, -258, -780 etc etc
Currently, in order to fix that 2 Flags can be used:
SOBEL_FIX_TRUNCATE (Value = 0). Simply truncates the values of Gx and Gy to -1 (if the result is too dark) or 1 (If the result is too bright)
SOBEL_FIX_SIMPLE (Value = 1). Fix the values adjusting the extrapolation to it stays within the limits.
The math envolving is this fix is: FixedValue = OldValue/(255*(floor(|OldValue|/255)+1))
pDegree(Out/Optional): A Pointer to a variable to store the Angle formed by the pixel. Sobel Operator can be used to calculate the angle (direction) of the Pixel.
It outputs the angle in Degrees. The format of the output is always in real4 format from 0º to 360º.
This parameter can be &NULL if you don´t want to exprt the angle.
Return Values:
The function does not return any value.
Example of usage:
[MyMatrix: F$ 25, 200, 30
100, 45, 75
0, 81, 255]
call SobelGetGPlusEx MyMatrix, MySobelX, MySobelY, MySobel, SOBEL_FIX_SIMPLE, &NULL
;;
Proc SobelGetGPLusEx:
Arguments @pMatrix, @pOutSobelX, @pOutSobelY, @pOutSobel, @FixType, @pDegree
Local @pReturn, @DataCheck, @Divisor, @FractionX, @FractionY, @Floor, @YAxis, @XAxis
Structure @TempAxis 16, @SobelYAxisDis 0, @SobelXAxisDis 8
Uses edi, edx, ecx, eax
finit
mov edi D@pMatrix
; To calculate Gx^2 later. Therefore Gx = M3+M9 + 2*(M6-M4) - (M7+M1)
fld F$edi+FloatMatricesInt.M6Dis | fsub F$edi+FloatMatricesInt.M4Dis | fmul R$Float_Two
fadd F$edi+FloatMatricesInt.M3Dis | fadd F$edi+FloatMatricesInt.M9Dis
fsub F$edi+FloatMatricesInt.M7Dis | fsub F$edi+FloatMatricesInt.M1Dis
lea ecx D@DataCheck | fistp F$ecx
.If D@DataCheck = 0
fldz | fstp F@FractionX
.Else_If D@DataCheck <s 0-255 ; Blacks. Edge too dark.
If D@FixType = SOBEL_FIX_TRUNCATE ; Truncate the value to -1. The default Sobel behaviour (It truncates to -255 or 255. But, since we are using normalized version we set the limits to -1 and 1)
fld R$Float_MinusOne | fstp F@FractionX
Else ; SOBEL_FIX_SIMPLE. FixedValue = OldValue/(255*(floor(|OldValue|/255)+1))
xor edx edx | mov ecx 255 | mov eax D@DataCheck | neg eax | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionX
End_If
.Else_If D@DataCheck >s 255 ; Whites. Edge too brigth.
If D@FixType = SOBEL_FIX_TRUNCATE ; Truncate the value to -1. The default Sobel behaviour (It truncates to -255 or 255. But, since we are using normalized version we set the limits to -1 and 1)
fld1 | fstp F@FractionX
Else ; SOBEL_FIX_SIMPLE. FixedValue = OldValue/(255*(floor(OldValue/255)+1))
xor edx edx | mov ecx 255 | mov eax D@DataCheck | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionX
End_If
.Else
fild D@DataCheck | fmul R$Float_One_255 | fstp F@FractionX
.End_If
mov eax D@pOutSobelX | mov ecx D@FractionX | mov D$eax ecx
; To calculate Gy^2 later. Therefore Gy = M7+M9 + 2*(M8-M2) - (M3+M1)
fld F$edi+FloatMatricesInt.M8Dis | fsub F$edi+FloatMatricesInt.M2Dis | fmul R$Float_Two
fadd F$edi+FloatMatricesInt.M7Dis | fadd F$edi+FloatMatricesInt.M9Dis
fsub F$edi+FloatMatricesInt.M3Dis | fsub F$edi+FloatMatricesInt.M1Dis
lea ecx D@DataCheck | fistp F$ecx
.If D@DataCheck = 0
fldz | fstp F@FractionY
.Else_If D@DataCheck <s 0-255 ; Blacks. Edge too dark. FixedValue = OldValue/(255*(floor(|OldValue|/255)+1))
If D@FixType = SOBEL_FIX_TRUNCATE ; Truncate the value to -1. The default Sobel behaviour (It truncates to -255 or 255. But, since we are using normalized version we set the limits to -1 and 1)
fld R$Float_MinusOne | fstp F@FractionY
Else ; SOBEL_FIX_SIMPLE ; FixedValue = OldValue/(255*(floor(|OldValue|/255)+1))
xor edx edx | mov ecx 255 | mov eax D@DataCheck | neg eax | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionY
End_If
.Else_If D@DataCheck >s 255 ; Whites. Edge too brigth. FixedValue = OldValue/(255*(floor(OldValue/255)+1))
If D@FixType = SOBEL_FIX_TRUNCATE ; Truncate the value to -1. The default Sobel behaviour (It truncates to -255 or 255. But, since we are using normalized version we set the limits to -1 and 1)
fld1 | fstp F@FractionY
Else ; SOBEL_FIX_SIMPLE; FixedValue = OldValue/(255*(floor(OldValue/255)+1))
xor edx edx | mov ecx 255 | mov eax D@DataCheck | div ecx | inc eax | imul eax 255 | mov D@Divisor eax
fild D@DataCheck | fidiv F@Divisor | fstp F@FractionY
End_If
.Else
fild D@DataCheck | fmul R$Float_One_255 | fstp F@FractionY
.End_If
mov eax D@pOutSobelY | mov ecx D@FractionY | mov D$eax ecx
; Output the Angle (in degrees) if needed
If D@pDegree <> 0
lea eax D@SobelYAxisDis | fld F@FractionY | fstp R$eax
lea edx D@SobelXAxisDis | fld F@FractionX | fstp R$edx
call atan2 eax, edx, &True
mov eax D@pDegree | fstp F$eax
End_If
; And finally create the Sobel G after Gx and Gy are already fixed
; Soble = sqrt((255*FractionX)^2+(255*FractionY)^2)) = G = sqrt(Gx^2+Gy^2)
; since FractionX and FractionY can have a maximum of 1 and -1, therefore sobleMax = (255/sqrt(2)) * sqrt(FractionX^2+FractionY^2)
fld F@FractionX | fmul ST0 ST0 | fld F@FractionY | fmul ST0 ST0 | faddp ST1 ST0 | fsqrt | fmul R$Float_SobleVarG
lea edx D@pReturn
fistp F$edx
mov eax D@pReturn
If eax > 255
mov eax 255
End_If
mov ecx D@pOutSobel | mov D$ecx eax
EndP
[Float_AtanPiFactor: R$ (180/3.1415926535897932384626433832795)]
[Float360: R$ 360]
; Macros used to simulate If/Else/EndIf using Fpu
[Fpu_If | fld #3 | fld #1 | fcompp | fstsw ax | fwait | sahf | jn#2 R0>>]
[Fpu_Else_If | jmp R5>> | R0: | fld #3 | fld #1 | fcompp | fstsw ax | fwait | sahf | jn#2 R0>>]
[Fpu_Else | jmp R5>> | R0:]
[Fpu_End_If | R0: | R5:]
Proc atan2:
Arguments @pY, @pX, @ConvDegree
Structure @TempStorage 16, @HueDis 0
Uses eax, ebx, ecx
mov ebx D@pY
mov ecx D@pX
fld R$ebx
fld R$ecx
fpatan
fstsw ax
wait
shr ax 1
jnb L2>
fclex | stc | xor eax eax | ExitP
L2:
.If D@ConvDegree = &TRUE
fmul R$Float_AtanPiFactor | fst R@HueDis
Fpu_If R@HueDis < R$FloatZero
fadd R$Float360
Fpu_Else_If R@HueDis >= R$Float360
fsub R$Float360
Fpu_End_If
.End_If
clc
mov eax &TRUE
EndP
:thumbsup:
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. :badgrin:
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:thumbsup:
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. :badgrin:
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
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.
Laplacian is simple
0 -1 0
-1 4 -1
0 -1 0
Matrix ValuesI don't undestand, I need read all this post.
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
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
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 ]
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 = [ 255 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 = [ 255 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 ]
[ 0 0 255 ]
Pixels = [ 0 0 255 ] * Gx = -1020
[ 0 0 255 ]
[+1 0 -1 ] [ 0 0 255 ]
G = [+2 0 -2 ] * Pixels = [ 0 0 255 ] = 1*0 + 0*0 + -1*255 + 2*0 + 0*0 + -2*255 + 1*0 + 0*0 + -1*255 = -1020 ( we found a maximum edge on the Gx axis but in the opposed sign, because its´in the opposite direction)
[+1 0 -1 ] [ 0 0 255 ]
Proc SetSobelData:
Arguments @pColorData, @ImgWidth, @ImgHeight, @pOutSobelX, @pOutSobelY
Uses esi, edx, ecx, ebx, edi
mov esi D@pColorData
; initialize min and max values
mov D$EdgeSobelPlus+EdgeSobelPlus.SobelGMinDis 255
mov D$EdgeSobelPlus+EdgeSobelPlus.SobelGMaxDis 0
; initialize SobelPlus structure
call PrepareSobelEdges EdgeSobelPlus, 0, 0, EDGE_FILTER_SOBEL
(...)
mov ebx D@Imgwidth | shl ebx 2 | mov D@NextLine ebx
mov D@Y 0 | mov eax D@Y
...While eax < D@Imgheight
mov D@X 0
mov edx D@X
..While edx < D@Imgwidth
; fill our matrix with the inputed pixels
call FillMatrix3_3 esi, edx, D@Y, TmpMatrix, D@ImgWidth, D@ImgHeight
call FindSobelEdgesEx TmpMatrix, EdgeSobelPlus, &FALSE
; I´ll Optimize this later removing the labels and using registers, instead
mov eax D$EdgeSobelPlus+EdgeSobelPlus.SobelGxDis | mov edi D@CurSobelX | mov D$edi eax | add D@CurSobelX 4
mov eax D$EdgeSobelPlus+EdgeSobelPlus.SobelGyDis | mov edi D@CurSobelY | mov D$edi eax | add D@CurSobelY 4
mov eax D$EdgeSobelPlus+EdgeSobelPlus.SobelGDis | mov edi D@CurSobel | add D@CurSobel 4
On eax < D$EdgeSobelPlus+EdgeSobelPlus.SobelGMinDis, mov D$EdgeSobelPlus+EdgeSobelPlus.SobelGMinDis eax
On eax > D$EdgeSobelPlus+EdgeSobelPlus.SobelGMaxDis, mov D$EdgeSobelPlus+EdgeSobelPlus.SobelGMaxDis eax
; Not necessary to set the output here for the scope of the watermar detector. I putted here only because i wanted to see the results 1st
mov B$edi+BGRA.RedDis al
mov B$edi+BGRA.GreenDis al
mov B$edi+BGRA.BlueDis al
mov B$edi+BGRA.AlphaDis 255
inc edx
..End_While
inc D@Y
mov eax D@Y
add esi D@NextLine
...End_While
(...)
EndP
;;
PrepareSobelEdges
This function precalculates the values of the normalization to be used on a edge detection function according to the operators values you feed in
Parameters:
pSobelStruct - Pointer to a EdgeSobelPlus Structure that will hold and deal with the several values retrieved during a sobel, scharr or any other edge operators we want
SobelValA - Is the 'magic' numbers used in an edge operator. ex: In sobel, this value is 0-1, for scharr, this value is 47 and so on. It is the number related to the M1 position on the matrix index.
You can set whatever value you want to create your own Edge operator indexes. To do that, you must define the TypeFlag member as EDGE_FILTER_USER.
Otherwise, if you want a predefined operator (sobel), you can set this member to 0 and use the EDGE_FILTER_SOBEL equate in the TypeFlag parameter.
SobelValB - Is the 'magic' numbers used in an edge operator. ex: In sobel, this value is 0-2, for scharr, this value is 162 and so on. It is the number related to the M4 position on the matrix index.
You can set whatever value you want to create your own Edge operator indexes. To do that, you must define the TypeFlag member as EDGE_FILTER_USER.
Otherwise, if you want a predefined operator (sobel), you can set this member to 0 and use the EDGE_FILTER_SOBEL equate in the TypeFlag parameter
TypeFlag - A set of predefined Flags to use for specific operators. Currently, you can use any of these 3 equates:
Equate Name Equate Value Meaning
EDGE_FILTER_USER 0 Use this flag, if you want to create your own pre-defined operator. You must fill the values on SobelValA and SobelValB
EDGE_FILTER_SOBEL 1 Use this flag to enable only the Sobel Operator. The values on SobelValA and SobelValB will be ignored (you can, as well, simply set them to zero. But for convention, better is use 0 insetad)
EDGE_FILTER_SCHARR 2 Use this flag to enable only the Sobel Operator. The values on SobelValA and SobelValB will be ignored (you can, as well, simply set them to zero. But for convention, better is use 0 insetad)
;;
Proc PrepareSobelEdges:
Arguments @pSobelStruct, @SobelValA, @SobelValD, @TypeFlag
Local @MatrixValA, @MatrixValD
Uses edi
mov eax D@SobelValA | mov D@MatrixValA eax
mov eax D@SobelValD | mov D@MatrixValD eax
.If D@TypeFlag = EDGE_FILTER_SOBEL
mov D@MatrixValA 0-1
mov D@MatrixValD 0-2
.Else_If D@TypeFlag = EDGE_FILTER_SCHARR
mov D@MatrixValA 47
mov D@MatrixValD 161; 162 (official value). 161 is the best value i found because the sum of MatrixA+matrixB applying the formula below gives 255
.Else
; if both are zero, avoid division by zero on the normalization ratio
mov edi D@MatrixValA
.Test_If_Not_And eax eax, edi edi
mov edi D@pSobelStruct
fldz | fst F$edi+EdgeSobelPlus.MatrixValADis | fst F$edi+EdgeSobelPlus.MatrixValBDis | fstp F$edi+EdgeSobelPlus.NormalizeRatioDis
ExitP
.Test_End
.End_If
mov edi D@pSobelStruct ; All Ok, let´s calculate and put their values
fild D@MatrixValA | fst F$edi+EdgeSobelPlus.MatrixValADis | fabs | fmul R$Float_Two
fild D@MatrixValD | fst F$edi+EdgeSobelPlus.MatrixValBDis | fabs | faddp ST1 ST0
fmul R$Float_255 | fld1 | fdivrp ST0 ST1 | fstp F$edi+EdgeSobelPlus.NormalizeRatioDis
EndP
[EdgeSobelPlus:
EdgeSobelPlus.ValA: F$ 0 ; Is the magic number used in Sobel in Position M1. Ex: 0-1
EdgeSobelPlus.ValB: F$ 0 ; Is the magic number used in Sobel in Position M4. Ex: 0-2
EdgeSobelPlus.NormalizeRatio: F$ 0 ; Is the ratio used to normalize the edge detection function.
EdgeSobelPlus.SobelGx: F$ 0 ; The value of Sobel Gx. Real4 format
EdgeSobelPlus.SobelGy: F$ 0 ; The value of Sobel Gy. Real4 format
EdgeSobelPlus.SobelG: F$ 0 ; The value of Sobel G. Real4 format
EdgeSobelPlus.Angle: F$ 0 ; The Angle (in degrees)of Sobel G. Real4 format
EdgeSobelPlus.SobelGxMin: F$ 0 ; The Minimum Value of Sobel Gx inside a image. Real4 format
EdgeSobelPlus.SobelGyMin: F$ 0 ; The Minimum Value of Sobel Gy inside a image. Real4 format
EdgeSobelPlus.SobelGMin: F$ 0 ; The Minimum Value of Sobel G inside a image. Real4 format
EdgeSobelPlus.SobelGxMax: F$ 0 ; The Maximum Value of Sobel Gx inside a image. Real4 format
EdgeSobelPlus.SobelGyMax: F$ 0 ; The Maximum Value of Sobel Gy inside a image. Real4 format
EdgeSobelPlus.SobelGMax: F$ 0 ; The Maximum Value of Sobel G inside a image. Real4 format
EdgeSobelPlus.AngleMin: F$ 0 ; The Maximum Value of Sobel G inside a image. Real4 format
EdgeSobelPlus.AngleMax: F$ 0 ; The Maximum Value of Sobel G inside a image. Real4 format
]
; Equates used
[EdgeSobelPlus.MatrixValADis 0
EdgeSobelPlus.MatrixValBDis 4
EdgeSobelPlus.NormalizeRatioDis 8
EdgeSobelPlus.SobelGxDis 12
EdgeSobelPlus.SobelGyDis 16
EdgeSobelPlus.SobelGDis 20
EdgeSobelPlus.AngleDis 24
EdgeSobelPlus.SobelGxMinDis 28
EdgeSobelPlus.SobelGyMinDis 32
EdgeSobelPlus.SobelGMinDis 36
EdgeSobelPlus.SobelGxMaxDis 40
EdgeSobelPlus.SobelGyMaxDis 44
EdgeSobelPlus.SobelGMaxDis 48
EdgeSobelPlus.AngleMinDis 52
EdgeSobelPlus.AngleMaxDis 56]
[Size_of_EdgeSobelPlus 60]
;;
pMatrix (in)- Pointer to the inputed pixels stored on a 3x3 matrix
pSobelStruct (in/out) - Pointer to a EdgeSobelPlus structure that will grab the values predefined to calcaulate the sobel, and fill thenecessay data accordlying.
The values that will be filled here are: SobelGx, SobelGy, Sobel G and Angle (in degrees). All value are stored on Real4 format, except Sobel G that is exported as integer
Sobel Gx and Sobel Gy (or shcarr gx/gy or whatever other operator you used), will be stored on the normalized values between -1 to 1 and always biased on 255.
So the range is, in fact from -255 to 255. To retrieve this values (the same as the pixel) all you have to do later is multiply the values found to 255.
UseDegree - If yopu want the function to export the angle of the pixel set this flag to TRUE. Otherwise set it to FALSE.
;;
Proc FindSobelEdgesEx:
Arguments @pMatrix, @pSobelStruct, @UseDegree
Local @DataCheck
Structure @TempAxis 16, @SobelYAxisDis 0, @SobelXAxisDis 8
Uses esi, edi, edx, ecx, eax
finit
mov esi D@pMatrix
; calculate normalization is already done in PrepareSobelEdges function
;Gx = A*(M1+M7-M3-M9) + D*(M4-M6); where A = 47, D = 162
mov edi D@pSobelStruct
fld F$esi+FloatMatricesInt.M1Dis | fadd F$esi+FloatMatricesInt.M7Dis | fsub F$esi+FloatMatricesInt.M3Dis | fsub F$esi+FloatMatricesInt.M9Dis | fmul F$edi+EdgeSobelPlus.MatrixValADis;D@MatrixValA
fld F$esi+FloatMatricesInt.M4Dis | fsub F$esi+FloatMatricesInt.M6Dis | fmul F$edi+EdgeSobelPlus.MatrixValBDis | faddp ST1 ST0
; normalize to 1020, 765 and multiply the result to 255 to we get the proper integer are already predefined in PrepareSobelEdges function
fmul F$edi+EdgeSobelPlus.NormalizeRatioDis | fstp F$edi+EdgeSobelPlus.SobelGxDis
; Gy = A*(M1+M3-M7-M9) + B*(M2-M8)
fld F$esi+FloatMatricesInt.M1Dis | fadd F$esi+FloatMatricesInt.M3Dis | fsub F$esi+FloatMatricesInt.M7Dis | fsub F$esi+FloatMatricesInt.M9Dis | fmul F$edi+EdgeSobelPlus.MatrixValADis; | fmul D@MatrixValA
fld F$esi+FloatMatricesInt.M2Dis | fsub F$esi+FloatMatricesInt.M8Dis | fmul F$edi+EdgeSobelPlus.MatrixValBDis | faddp ST1 ST0
; normalize to 1020, 765 etc and multiply the result to 255 to we get the proper integer are already predefined in PrepareSobelEdges function
fmul F$edi+EdgeSobelPlus.NormalizeRatioDis | fstp F$edi+EdgeSobelPlus.SobelGyDis
; Calculate the angle if needed
If D@UseDegree <> 0
lea eax D@SobelYAxisDis | fld F$edi+EdgeSobelPlus.SobelGyDis | fstp R$eax
lea edx D@SobelXAxisDis | fld F$edi+EdgeSobelPlus.SobelGxDis | fstp R$edx
call atan2 eax, edx, &True
fstp F$edi+EdgeSobelPlus.AngleDis
End_If
; And finally create the Sobel G after Gx and Gy are already fixed
; Soble = sqrt((255*FractionX)^2+(255*FractionY)^2)) = G = sqrt(Gx^2+Gy^2)
; since FractionX and FractionY can have a maximum of 1 and -1, therefore sobleMax = (255/sqrt(2)) * sqrt(FractionX^2+FractionY^2)
fld F$edi+EdgeSobelPlus.SobelGxDis | fmul ST0 ST0 | fld F$edi+EdgeSobelPlus.SobelGyDis | fmul ST0 ST0 | faddp ST1 ST0 | fsqrt | fmul R$Float_SobleVarG
fistp F$edi+EdgeSobelPlus.SobelGDis
mov eax D$edi+EdgeSobelPlus.SobelGDis
If eax > 255
mov eax 255
End_If
mov D$edi+EdgeSobelPlus.SobelGDis eax
EndP
Putting the 255 on the right side and do the same thing as above, produces always the same result, right ?
.QuotePutting the 255 on the right side and do the same thing as above, produces always the same result, right ?
The sign tells you if the edge is left or right, for Gx
Up or down for Gy.
I'm not that far yet with my programming.
Still at the stage of generating different gray color conversions. :tongue: ( 4 pixels at once, and handle image widths not divisible by 4 )
I have written the code this way, that every stage can access all next stages etc. to see all the possibilities on screen.
stage 1: different gray color conversions
stage 2: different lowpass filters
stage 3: different types of edge detection operators ( Sobel, Scharr etc )
stage 4: creating a Gx,Gy table from one of the previous stages
stage 5: finding the directions per edge pixel and give it a color for that direction. ( stage 5: just for the fun, it should look awesome )
But first step is, getting the math correct ( think I got it covered, else, have to study it more ) then after that, optimize it.
Still thinking what to do with G, maybe clamp it to get a higher intensity color, I will experiment with that when I reach that stage.
[NativeMatrix:
WSData_NTSC_RGB_C.Name: D$ Sz_NTSC_RGB
WsData_WhiteRef_NTSC_RGB: D$ WS_C_OBS2
WSData_NTSC_RGB_C_Red_M1: R$ 0.6068909 WSData_NTSC_RGB_C_Green_M2: R$ 0.1735011 WSData_NTSC_RGB_C_Blue_M3: R$ 0.200348
WSData_NTSC_RGB_C_Red_M4: R$ 0.2989164 WSData_NTSC_RGB_C_Green_M5: R$ 0.586599 WSData_NTSC_RGB_C_Blue_M6: R$ 0.1144846
WSData_NTSC_RGB_C_Red_M7: R$ 0 WSData_NTSC_RGB_C_Green_M8: R$ 0.0660957 WSData_NTSC_RGB_C_Blue_M9: R$ 1.1162243
WSData_GammaDecEnc_NTSC_RGB_C.Gamma: R$ 2.2
WSData_GammaDecEnc_NTSC_RGB_C.Offset: R$ 0.055
WSData_Best_RGB_D50.Name: D$ Sz_Best_RGB
WsData_WhiteRef_Best_RGB: D$ WS_D50_OBS2
WSData_Best_RGB_D50_Red_M1: R$ 0.6326696 WSData_Best_RGB_D50_Green_M2: R$ 0.2045558 WSData_Best_RGB_D50_Blue_M3: R$ 0.1269946
WSData_Best_RGB_D50_Red_M4: R$ 0.2284569 WSData_Best_RGB_D50_Green_M5: R$ 0.7373523 WSData_Best_RGB_D50_Blue_M6: R$ 0.0341908
WSData_Best_RGB_D50_Red_M7: R$ 0 WSData_Best_RGB_D50_Green_M8: R$ 0.0095142 WSData_Best_RGB_D50_Blue_M9: R$ 0.8156958
WSData_GammaDecEnc_Best_RGB_D50.Gamma: R$ 2.2
WSData_GammaDecEnc_Best_RGB_D50.Offset: R$ 0.055
WSData_Beta_RGB_D50.Name: D$ Sz_Beta_RGB
WsData_WhiteRef_Beta_RGB: D$ WS_D50_OBS2
WSData_Beta_RGB_D50_Red_M1: R$ 0.6712537 WSData_Beta_RGB_D50_Green_M2: R$ 0.1745834 WSData_Beta_RGB_D50_Blue_M3: R$ 0.1183829
WSData_Beta_RGB_D50_Red_M4: R$ 0.3032726 WSData_Beta_RGB_D50_Green_M5: R$ 0.6637861 WSData_Beta_RGB_D50_Blue_M6: R$ 0.0329413
WSData_Beta_RGB_D50_Red_M7: R$ 0 WSData_Beta_RGB_D50_Green_M8: R$ 0.040701 WSData_Beta_RGB_D50_Blue_M9: R$ 0.784509
WSData_GammaDecEnc_Beta_RGB_D50.Gamma: R$ 2.2
WSData_GammaDecEnc_Beta_RGB_D50.Offset: R$ 0.055
WSData_ColorMatch_RGB_D50.Name: D$ Sz_ColorMatch
WsData_WhiteRef_ColorMatch_RGB: D$ WS_D50_OBS2
WSData_ColorMatch_RGB_D50_Red_M1: R$ 0.5093439 WSData_ColorMatch_RGB_D50_Green_M2: R$ 0.3209071 WSData_ColorMatch_RGB_D50_Blue_M3: R$ 0.133969
WSData_ColorMatch_RGB_D50_Red_M4: R$ 0.274884 WSData_ColorMatch_RGB_D50_Green_M5: R$ 0.6581315 WSData_ColorMatch_RGB_D50_Blue_M6: R$ 0.0669845
WSData_ColorMatch_RGB_D50_Red_M7: R$ 0.0242545 WSData_ColorMatch_RGB_D50_Green_M8: R$ 0.108782 WSData_ColorMatch_RGB_D50_Blue_M9: R$ 0.6921735
WSData_GammaDecEnc_ColorMatch_RGB_D50.Gamma: R$ 1.8
WSData_GammaDecEnc_ColorMatch_RGB_D50.Offset: R$ 0.055
WSData_Don_RGB_4_D50.Name: D$ Sz_Don_RGB4
WsData_WhiteRef_Don_RGB_4: D$ WS_D50_OBS2
WSData_Don_RGB_4_D50_Red_M1: R$ 0.6457711 WSData_Don_RGB_4_D50_Green_M2: R$ 0.1933511 WSData_Don_RGB_4_D50_Blue_M3: R$ 0.1250978
WSData_Don_RGB_4_D50_Red_M4: R$ 0.2783496 WSData_Don_RGB_4_D50_Green_M5: R$ 0.6879702 WSData_Don_RGB_4_D50_Blue_M6: R$ 0.0336802
WSData_Don_RGB_4_D50_Red_M7: R$ 0.0037113 WSData_Don_RGB_4_D50_Green_M8: R$ 0.0179861 WSData_Don_RGB_4_D50_Blue_M9: R$ 0.8035126
WSData_GammaDecEnc_Don_RGB_4_D50.Gamma: R$ 2.2
WSData_GammaDecEnc_Don_RGB_4_D50.Offset: R$ 0.055
WSData_ECI_RGB_D50.Name: D$ Sz_ECI_RGB
WsData_WhiteRef_ECI_RGB: D$ WS_D50_OBS2
WSData_ECI_RGB_D50_Red_M1: R$ 0.6502043 WSData_ECI_RGB_D50_Green_M2: R$ 0.1780774 WSData_ECI_RGB_D50_Blue_M3: R$ 0.1359383
WSData_ECI_RGB_D50_Red_M4: R$ 0.3202499 WSData_ECI_RGB_D50_Green_M5: R$ 0.6020711 WSData_ECI_RGB_D50_Blue_M6: R$ 0.077679
WSData_ECI_RGB_D50_Red_M7: R$ 0 WSData_ECI_RGB_D50_Green_M8: R$ 0.067839 WSData_ECI_RGB_D50_Blue_M9: R$ 0.757371
WSData_GammaDecEnc_ECI_RGB_D50.Gamma: R$ 1.8
WSData_GammaDecEnc_ECI_RGB_D50.Offset: R$ 0.055
WSData_Ekta_Space_PS5_D50.Name: D$ Sz_Ekta_Space_PS5
WsData_WhiteRef_Ekta_Space_PS5: D$ WS_D50_OBS2
WSData_Ekta_Space_PS5_D50_Red_M1: R$ 0.5938914 WSData_Ekta_Space_PS5_D50_Green_M2: R$ 0.2729801 WSData_Ekta_Space_PS5_D50_Blue_M3: R$ 0.0973485
WSData_Ekta_Space_PS5_D50_Red_M4: R$ 0.2606286 WSData_Ekta_Space_PS5_D50_Green_M5: R$ 0.7349465 WSData_Ekta_Space_PS5_D50_Blue_M6: R$ 0.0044249
WSData_Ekta_Space_PS5_D50_Red_M7: R$ 0 WSData_Ekta_Space_PS5_D50_Green_M8: R$ 0.0419969 WSData_Ekta_Space_PS5_D50_Blue_M9: R$ 0.7832131
WSData_GammaDecEnc_Ekta_Space_PS5_D50.Gamma: R$ 2.2
WSData_GammaDecEnc_Ekta_Space_PS5_D50.Offset: R$ 0.055
WsData_ProPhoto_RGB.Name: D$ Sz_ProPhoto
WsData_WhiteRef_ProPhoto_RGB: D$ WS_D50_OBS2
WSData_ProPhoto_RGB_D50_Red_M1: R$ 0.7976749 WSData_ProPhoto_RGB_D50_Green_M2: R$ 0.1351917 WSData_ProPhoto_RGB_D50_Blue_M3: R$ 0.0313534
WSData_ProPhoto_RGB_D50_Red_M4: R$ 0.2880402 WSData_ProPhoto_RGB_D50_Green_M5: R$ 0.7118741 WSData_ProPhoto_RGB_D50_Blue_M6: R$ 0.0000857
WSData_ProPhoto_RGB_D50_Red_M7: R$ 0 WSData_ProPhoto_RGB_D50_Green_M8: R$ 0 WSData_ProPhoto_RGB_D50_Blue_M9: R$ 0.82521
WSData_GammaDecEnc_ProPhoto_RGB_D50.Gamma: R$ 1.8
WSData_GammaDecEnc_ProPhoto_RGB_D50.Offset: R$ 0.055
WSData_Wide_Gamut_RGB_D50.Name: D$ Sz_Wide_Gamut
WsData_WhiteRef_Wide_Gamut_RGB: D$ WS_D50_OBS2
WSData_Wide_Gamut_RGB_D50_Red_M1: R$ 0.7161046 WSData_Wide_Gamut_RGB_D50_Green_M2: R$ 0.1009296 WSData_Wide_Gamut_RGB_D50_Blue_M3: R$ 0.1471858
WSData_Wide_Gamut_RGB_D50_Red_M4: R$ 0.2581874 WSData_Wide_Gamut_RGB_D50_Green_M5: R$ 0.7249378 WSData_Wide_Gamut_RGB_D50_Blue_M6: R$ 0.0168748
WSData_Wide_Gamut_RGB_D50_Red_M7: R$ 0 WSData_Wide_Gamut_RGB_D50_Green_M8: R$ 0.0517813 WSData_Wide_Gamut_RGB_D50_Blue_M9: R$ 0.7734287
WSData_GammaDecEnc_Wide_Gamut_RGB_D50.Gamma: R$ 2.2
WSData_GammaDecEnc_Wide_Gamut_RGB_D50.Offset: R$ 0.055
WSData_Adobe_RGB_1998.Name: D$ Sz_Adobe_RGB_1998
WsData_WhiteRef_Adobe_RGB_1998: D$ WS_D65_OBS2
WSData_Adobe_RGB_1998_D65_Red_M1: R$ 0.5767309 WSData_Adobe_RGB_1998_D65_Green_M2: R$ 0.185554 WSData_Adobe_RGB_1998_D65_Blue_M3: R$ 0.1881851
WSData_Adobe_RGB_1998_D65_Red_M4: R$ 0.2973769 WSData_Adobe_RGB_1998_D65_Green_M5: R$ 0.6273491 WSData_Adobe_RGB_1998_D65_Blue_M6: R$ 0.075274
WSData_Adobe_RGB_1998_D65_Red_M7: R$ 0.0270343 WSData_Adobe_RGB_1998_D65_Green_M8: R$ 0.0706872 WSData_Adobe_RGB_1998_D65_Blue_M9: R$ 0.9911085
WSData_GammaDecEnc_Adobe_RGB_1998_D65.Gamma: R$ 2.19921875
WSData_GammaDecEnc_Adobe_RGB_1998_D65.Offset: R$ 0.055
WSData_AppleRGB_D65.Name: D$ Sz_AppleRGB
WsData_WhiteRef_AppleRGB: D$ WS_D65_OBS2
WSData_AppleRGB_D65_Red_M1: R$ 0.4497288 WSData_AppleRGB_D65_Green_M2: R$ 0.3162486 WSData_AppleRGB_D65_Blue_M3: R$ 0.1844926
WSData_AppleRGB_D65_Red_M4: R$ 0.2446525 WSData_AppleRGB_D65_Green_M5: R$ 0.6720283 WSData_AppleRGB_D65_Blue_M6: R$ 0.0833192
WSData_AppleRGB_D65_Red_M7: R$ 0.0251848 WSData_AppleRGB_D65_Green_M8: R$ 0.1411824 WSData_AppleRGB_D65_Blue_M9: R$ 0.9224628
WSData_GammaDecEnc_AppleRGB_D65.Gamma: R$ 1.8
WSData_GammaDecEnc_AppleRGB_D65.Offset: R$ 0.055
WSData_Bruce_RGB_D65.Name: D$ Sz_Bruce_RGB
WsData_WhiteRef_Bruce_RGB: D$ WS_D65_OBS2
WSData_Bruce_RGB_D65_Red_M1: R$ 0.4674162 WSData_Bruce_RGB_D65_Green_M2: R$ 0.2944512 WSData_Bruce_RGB_D65_Blue_M3: R$ 0.1886026
WSData_Bruce_RGB_D65_Red_M4: R$ 0.2410115 WSData_Bruce_RGB_D65_Green_M5: R$ 0.6835475 WSData_Bruce_RGB_D65_Blue_M6: R$ 0.075441
WSData_Bruce_RGB_D65_Red_M7: R$ 0.0219101 WSData_Bruce_RGB_D65_Green_M8: R$ 0.0736128 WSData_Bruce_RGB_D65_Blue_M9: R$ 0.9933071
WSData_GammaDecEnc_Bruce_RGB_D65.Gamma: R$ 2.2
WSData_GammaDecEnc_Bruce_RGB_D65.Offset: R$ 0.055
WSData_PAL_SECAM_RGB_D65.Name: D$ Sz_PAL_SECAM
WsData_WhiteRef_PAL_SECAM_RGB: D$ WS_D65_OBS2
WSData_PAL_SECAM_RGB_D65_Red_M1: R$ 0.430619 WSData_PAL_SECAM_RGB_D65_Green_M2: R$ 0.3415419 WSData_PAL_SECAM_RGB_D65_Blue_M3: R$ 0.1783091
WSData_PAL_SECAM_RGB_D65_Red_M4: R$ 0.2220379 WSData_PAL_SECAM_RGB_D65_Green_M5: R$ 0.7066384 WSData_PAL_SECAM_RGB_D65_Blue_M6: R$ 0.0713237
WSData_PAL_SECAM_RGB_D65_Red_M7: R$ 0.0201853 WSData_PAL_SECAM_RGB_D65_Green_M8: R$ 0.1295504 WSData_PAL_SECAM_RGB_D65_Blue_M9: R$ 0.9390943
WSData_GammaDecEnc_PAL_SECAM_RGB_D65.Gamma: R$ 2.2
WSData_GammaDecEnc_PAL_SECAM_RGB_D65.Offset: R$ 0.055
WSData_SMPTE_C_RGB_D65.Name: D$ Sz_SMPTE_C
WsData_WhiteRef_SMPTE_C_RGB: D$ WS_D65_OBS2
WSData_SMPTE_C_RGB_D65_Red_M1: R$ 0.3935891 WSData_SMPTE_C_RGB_D65_Green_M2: R$ 0.3652497 WSData_SMPTE_C_RGB_D65_Blue_M3: R$ 0.1916312
WSData_SMPTE_C_RGB_D65_Red_M4: R$ 0.2124132 WSData_SMPTE_C_RGB_D65_Green_M5: R$ 0.7010437 WSData_SMPTE_C_RGB_D65_Blue_M6: R$ 0.0865431
WSData_SMPTE_C_RGB_D65_Red_M7: R$ 0.0187423 WSData_SMPTE_C_RGB_D65_Green_M8: R$ 0.1119313 WSData_SMPTE_C_RGB_D65_Blue_M9: R$ 0.9581564
WSData_GammaDecEnc_SMPTE_C_RGB_D65.Gamma: R$ 2.2
WSData_GammaDecEnc_SMPTE_C_RGB_D65.Offset: R$ 0.055
WSData_sRGB_D65.Name: D$ Sz_sRGB
WsData_WhiteRef_sRGB: D$ WS_D65_OBS2
WSData_sRGB_D65_Red_M1: R$ 0.4124564 WSData_sRGB_D65_Green_M2: R$ 0.3575761 WSData_sRGB_D65_Blue_M3: R$ 0.1804375
WSData_sRGB_D65_Red_M4: R$ 0.2126729 WSData_sRGB_D65_Green_M5: R$ 0.7151522 WSData_sRGB_D65_Blue_M6: R$ 0.0721749
WSData_sRGB_D65_Red_M7: R$ 0.0193339 WSData_sRGB_D65_Green_M8: R$ 0.1191920 WSData_sRGB_D65_Blue_M9: R$ 0.9503041
WSData_GammaDecEnc_sRGB_D65.Gamma: R$ 2.4
WSData_GammaDecEnc_sRGB_D65.Offset: R$ 0.055
WSData_sRGB_D65_Red_HDTV.Name: D$ Sz_sRGB_HDTV
WsData_WhiteRef_sRGB_HDTV: D$ WS_D65_OBS2
WSData_sRGB_D65_Red_HDTV_M1: R$ 0.4124564 WSData_sRGB_D65_Green_HDTV_M2: R$ 0.3575761 WSData_sRGB_D65_Blue_HDTV_M3: R$ 0.1804375
WSData_sRGB_D65_Red_HDTV_M4: R$ 0.2126729 WSData_sRGB_D65_Green_HDTV_M5: R$ 0.7151522 WSData_sRGB_D65_Blue_HDTV_M6: R$ 0.0721749
WSData_sRGB_D65_Red_HDTV_M7: R$ 0.0193339 WSData_sRGB_D65_Green_HDTV_M8: R$ 0.1191920 WSData_sRGB_D65_Blue_HDTV_M9: R$ 0.9503041
WSData_GammaDecEnc_sRGB_D65_HDTV.Gamma: R$ 2.2
WSData_GammaDecEnc_sRGB_D65_HDTV.Offset: R$ 0.099
WSData_CIE_RGB_E.Name: D$ Sz_CIE_RGB
WsData_WhiteRef_CIE_RGB: D$ WS_E_OBS2
WSData_CIE_RGB_E_Red_M1: R$ 0.488718 WSData_CIE_RGB_E_Green_M2: R$ 0.3106803 WSData_CIE_RGB_E_Blue_M3: R$ 0.2006017
WSData_CIE_RGB_E_Red_M4: R$ 0.1762044 WSData_CIE_RGB_E_Green_M5: R$ 0.8129847 WSData_CIE_RGB_E_Blue_M6: R$ 0.0108109
WSData_CIE_RGB_E_Red_M7: R$ 0 WSData_CIE_RGB_E_Green_M8: R$ 0.0102048 WSData_CIE_RGB_E_Blue_M9: R$ 0.9897952
WSData_GammaDecEnc_CIE_RGB_E.Gamma: R$ 2.2
WSData_GammaDecEnc_CIE_RGB_E.Offset: R$ 0.055]
[Sz_NTSC_RGB: B$ "NTSC RGB", 0
Sz_Best_RGB: B$ "Best RGB", 0
Sz_Beta_RGB: B$ "Beta RGB", 0
Sz_ColorMatch: B$ "ColorMatch", 0
Sz_Don_RGB4: B$ "Don RGB4", 0
Sz_ECI_RGB: B$ "ECI RGB", 0
Sz_Ekta_Space_PS5: B$ "Ekta Space PS5", 0
Sz_ProPhoto: B$ "ProPhoto", 0
Sz_Wide_Gamut: B$ "Wide Gamut", 0
Sz_Adobe_RGB_1998: B$ "Adobe RGB (1998)", 0
Sz_AppleRGB: B$ "Apple RGB", 0
Sz_Bruce_RGB: B$ "Bruce RGB", 0
Sz_PAL_SECAM: B$ "PAL SECAM", 0
Sz_SMPTE_C: B$ "SMPTE-C", 0
Sz_sRGB: B$ "sRGB", 0
Sz_sRGB_HDTV: B$ "sRGB HDTV", 0
Sz_CIE_RGB: B$ "CIE RGB", 0]
[WHITE_SPACE_A_X_OBS2: R$ 109.850 WHITE_SPACE_A_Y_OBS2: R$ 100 WHITE_SPACE_A_Z_OBS2: R$ 35.585 ; Incandescent/tungsten
WHITE_SPACE_A_X_OBS10: R$ 111.144 WHITE_SPACE_A_Y_OBS10: R$ 100 WHITE_SPACE_A_Z_OBS10: R$ 35.200 ; Incandescent/tungsten
WHITE_SPACE_B_X_OBS2: R$ 99.0927 WHITE_SPACE_B_Y_OBS2: R$ 100 WHITE_SPACE_B_Z_OBS2: R$ 85.313 ; Old direct sunlight at noon
WHITE_SPACE_B_X_OBS10: R$ 99.178 WHITE_SPACE_B_Y_OBS10: R$ 100 WHITE_SPACE_B_Z_OBS10: R$ 84.3493 ; Old direct sunlight at noon
WHITE_SPACE_C_X_OBS2: R$ 98.074 WHITE_SPACE_C_Y_OBS2: R$ 100 WHITE_SPACE_C_Z_OBS2: R$ 118.232 ; Old daylight
WHITE_SPACE_C_X_OBS10: R$ 97.285 WHITE_SPACE_C_Y_OBS10: R$ 100 WHITE_SPACE_C_Z_OBS10: R$ 116.145 ; Old daylight
WHITE_SPACE_D50_X_OBS2: R$ 96.422 WHITE_SPACE_D50_Y_OBS2: R$ 100 WHITE_SPACE_D50_Z_OBS2: R$ 82.521 ; ICC profile PCS
WHITE_SPACE_D50_X_OBS10: R$ 96.720 WHITE_SPACE_D50_Y_OBS10: R$ 100 WHITE_SPACE_D50_Z_OBS10: R$ 81.427 ; ICC profile PCS
WHITE_SPACE_D55_X_OBS2: R$ 95.682 WHITE_SPACE_D55_Y_OBS2: R$ 100 WHITE_SPACE_D55_Z_OBS2: R$ 92.149 ; Mid-morning daylight
WHITE_SPACE_D55_X_OBS10: R$ 95.799 WHITE_SPACE_D55_Y_OBS10: R$ 100 WHITE_SPACE_D55_Z_OBS10: R$ 90.926 ; Mid-morning daylight
WHITE_SPACE_D65_X_OBS2: R$ 95.047 WHITE_SPACE_D65_Y_OBS2: R$ 100 WHITE_SPACE_D65_Z_OBS2: R$ 108.883 ; Daylight, sRGB, Adobe-RGB
WHITE_SPACE_D65_X_OBS10: R$ 94.811 WHITE_SPACE_D65_Y_OBS10: R$ 100 WHITE_SPACE_D65_Z_OBS10: R$ 107.304 ; Daylight, sRGB, Adobe-RGB
WHITE_SPACE_D75_X_OBS2: R$ 94.972 WHITE_SPACE_D75_Y_OBS2: R$ 100 WHITE_SPACE_D75_Z_OBS2: R$ 122.638 ; North sky daylight
WHITE_SPACE_D75_X_OBS10: R$ 94.416 WHITE_SPACE_D75_Y_OBS10: R$ 100 WHITE_SPACE_D75_Z_OBS10: R$ 120.641 ; North sky daylight
WHITE_SPACE_E_X_OBS2: R$ 100 WHITE_SPACE_E_Y_OBS2: R$ 100 WHITE_SPACE_E_Z_OBS2: R$ 100 ; Equal energy
WHITE_SPACE_E_X_OBS10: R$ 100 WHITE_SPACE_E_Y_OBS10: R$ 100 WHITE_SPACE_E_Z_OBS10: R$ 100 ; Equal energy
WHITE_SPACE_F1_X_OBS2: R$ 92.834 WHITE_SPACE_F1_Y_OBS2: R$ 100 WHITE_SPACE_F1_Z_OBS2: R$ 103.665 ; Daylight Fluorescent
WHITE_SPACE_F1_X_OBS10: R$ 94.791 WHITE_SPACE_F1_Y_OBS10: R$ 100 WHITE_SPACE_F1_Z_OBS10: R$ 103.191 ; Daylight Fluorescent
WHITE_SPACE_F2_X_OBS2: R$ 99.187 WHITE_SPACE_F2_Y_OBS2: R$ 100 WHITE_SPACE_F2_Z_OBS2: R$ 67.395 ; Cool fluorescent
WHITE_SPACE_F2_X_OBS10: R$ 103.280 WHITE_SPACE_F2_Y_OBS10: R$ 100 WHITE_SPACE_F2_Z_OBS10: R$ 69.026 ; Cool fluorescent
WHITE_SPACE_F3_X_OBS2: R$ 103.754 WHITE_SPACE_F3_Y_OBS2: R$ 100 WHITE_SPACE_F3_Z_OBS2: R$ 49.861 ; White Fluorescent
WHITE_SPACE_F3_X_OBS10: R$ 108.968 WHITE_SPACE_F3_Y_OBS10: R$ 100 WHITE_SPACE_F3_Z_OBS10: R$ 51.965 ; White Fluorescent
WHITE_SPACE_F4_X_OBS2: R$ 109.147 WHITE_SPACE_F4_Y_OBS2: R$ 100 WHITE_SPACE_F4_Z_OBS2: R$ 38.813 ; Warm White Fluorescent
WHITE_SPACE_F4_X_OBS10: R$ 114.961 WHITE_SPACE_F4_Y_OBS10: R$ 100 WHITE_SPACE_F4_Z_OBS10: R$ 40.963 ; Warm White Fluorescent
WHITE_SPACE_F5_X_OBS2: R$ 90.872 WHITE_SPACE_F5_Y_OBS2: R$ 100 WHITE_SPACE_F5_Z_OBS2: R$ 98.723 ; Daylight Fluorescent
WHITE_SPACE_F5_X_OBS10: R$ 93.369 WHITE_SPACE_F5_Y_OBS10: R$ 100 WHITE_SPACE_F5_Z_OBS10: R$ 98.636 ; Daylight Fluorescent
WHITE_SPACE_F6_X_OBS2: R$ 97.309 WHITE_SPACE_F6_Y_OBS2: R$ 100 WHITE_SPACE_F6_Z_OBS2: R$ 60.191 ; Lite White Fluorescent
WHITE_SPACE_F6_X_OBS10: R$ 102.148 WHITE_SPACE_F6_Y_OBS10: R$ 100 WHITE_SPACE_F6_Z_OBS10: R$ 62.074 ; Lite White Fluorescent
WHITE_SPACE_F7_X_OBS2: R$ 95.044 WHITE_SPACE_F7_Y_OBS2: R$ 100 WHITE_SPACE_F7_Z_OBS2: R$ 108.755 ; Daylight fluorescent, D65 simulator
WHITE_SPACE_F7_X_OBS10: R$ 95.792 WHITE_SPACE_F7_Y_OBS10: R$ 100 WHITE_SPACE_F7_Z_OBS10: R$ 107.687 ; Daylight fluorescent, D65 simulator
WHITE_SPACE_F8_X_OBS2: R$ 96.413 WHITE_SPACE_F8_Y_OBS2: R$ 100 WHITE_SPACE_F8_Z_OBS2: R$ 82.333 ; Sylvania F40, D50 simulator
WHITE_SPACE_F8_X_OBS10: R$ 97.115 WHITE_SPACE_F8_Y_OBS10: R$ 100 WHITE_SPACE_F8_Z_OBS10: R$ 81.135 ; Sylvania F40, D50 simulator
WHITE_SPACE_F9_X_OBS2: R$ 100.365 WHITE_SPACE_F9_Y_OBS2: R$ 100 WHITE_SPACE_F9_Z_OBS2: R$ 67.868 ; Cool White Fluorescent
WHITE_SPACE_F9_X_OBS10: R$ 102.116 WHITE_SPACE_F9_Y_OBS10: R$ 100 WHITE_SPACE_F9_Z_OBS10: R$ 67.826 ; Cool White Fluorescent
WHITE_SPACE_F10_X_OBS2: R$ 96.174 WHITE_SPACE_F10_Y_OBS2: R$ 100 WHITE_SPACE_F10_Z_OBS2: R$ 81.712 ; Ultralume 50, Philips TL85
WHITE_SPACE_F10_X_OBS10: R$ 99.001 WHITE_SPACE_F10_Y_OBS10: R$ 100 WHITE_SPACE_F10_Z_OBS10: R$ 83.134 ; Ultralume 50, Philips TL85
WHITE_SPACE_F11_X_OBS2: R$ 100.966 WHITE_SPACE_F11_Y_OBS2: R$ 100 WHITE_SPACE_F11_Z_OBS2: R$ 64.370 ; Ultralume 40, Philips TL84
WHITE_SPACE_F11_X_OBS10: R$ 103.866 WHITE_SPACE_F11_Y_OBS10: R$ 100 WHITE_SPACE_F11_Z_OBS10: R$ 65.627 ; Ultralume 40, Philips TL84
WHITE_SPACE_F12_X_OBS2: R$ 108.046 WHITE_SPACE_F12_Y_OBS2: R$ 100 WHITE_SPACE_F12_Z_OBS2: R$ 39.228 ; Ultralume 30, Philips TL83
WHITE_SPACE_F12_X_OBS10: R$ 111.428 WHITE_SPACE_F12_Y_OBS10: R$ 100 WHITE_SPACE_F12_Z_OBS10: R$ 40.353] ; Ultralume 30, Philips TL83
[WhiteSpaceXDis 0
WhiteSpaceYDis 8
WhiteSpaceZDis 16]
[Size_Of_WhiteSpace 24]
[WS_A_OBS2 0
WS_A_OBS10 1
WS_B_OBS2 2
WS_B_OBS10 3
WS_C_OBS2 4
WS_C_OBS10 5
WS_D50_OBS2 6
WS_D50_OBS10 7
WS_D55_OBS2 8
WS_D55_OBS10 9
WS_D65_OBS2 10
WS_D65_OBS10 11
WS_D75_OBS2 12
WS_D75_OBS10 13
WS_E_OBS2 14
WS_E_OBS10 15
WS_F1_OBS2 16
WS_F1_OBS10 17
WS_F2_OBS2 18
WS_F2_OBS10 19
WS_F3_OBS2 20
WS_F3_OBS10 21
WS_F4_OBS2 22
WS_F4_OBS10 23
WS_F5_OBS2 24
WS_F5_OBS10 25
WS_F6_OBS2 26
WS_F6_OBS10 27
WS_F7_OBS2 28
WS_F7_OBS10 29
WS_F8_OBS2 30
WS_F8_OBS10 31
WS_F9_OBS2 32
WS_F9_OBS10 33
WS_F10_OBS2 34
WS_F10_OBS10 35
WS_F11_OBS2 36
WS_F11_OBS10 37
WS_F12_OBS2 38
WS_F12_OBS10 39]
; Native Working Spaces and their Matrices
[CS_MATRIX_NTSC_RGB_C 0]
[CS_MATRIX_BEST_RGB_D50 1]
[CS_MATRIX_BETA_RGB_D50 2]
[CS_MATRIX_COLORMATCH_RGB_D50 3]
[CS_MATRIX_DON_RGB_4_D50 4]
[CS_MATRIX_ECI_RGB_D50 5]
[CS_MATRIX_EKTA_SPACE_PS5_D50 6]
[CS_MATRIX_PROPHOTO_RGB_D50 7]
[CS_MATRIX_WIDE_GAMUT_RGB_D50 8]
[CS_MATRIX_ADOBE_RGB_1998_D65 9]
[CS_MATRIX_APPLERGB_D65 10]
[CS_MATRIX_BRUCE_RGB_D65 11]
[CS_MATRIX_PAL_SECAM_RGB_D65 12]
[CS_MATRIX_SMPTE_C_RGB_D65 13]
[CS_MATRIX_SRGB_D65 14]
[CS_MATRIX_SRGB_D65_HDTV 15]
[CS_MATRIX_CIE_RGB_E 16]
[NativeMatrix.NameDis 0
NativeMatrix.WhiteRefDis 4
NativeMatrix.RedM1Dis 8
NativeMatrix.GreenM2Dis 16
NativeMatrix.BlueM3Dis 24
NativeMatrix.RedM4Dis 32
NativeMatrix.GreenM5Dis 40
NativeMatrix.BlueM6Dis 48
NativeMatrix.RedM7Dis 56
NativeMatrix.GreenM8Dis 64
NativeMatrix.BlueM9Dis 72
NativeMatrix.GammaDecEnc.GammaDis 80
NativeMatrix.GammaDecEnc.OffsetDis 88
Size_Of_NativeMatrix 96]
;;
PrepareSobelEdges
This function precalculates the values of the normalization to be used on a edge detection function according to the operators values you feed in
Parameters:
pSobelStruct - Pointer to a EdgeSobelPlus Structure that will hold and deal with the several values retrieved during a sobel, scharr or any other edge operators we want
SobelValA - Is the 'magic' numbers used in an edge operator. ex: In sobel, this value is 0-1, for scharr, this value is 47 and so on. It is the number related to the M1 position on the matrix index.
You can set whatever value you want to create your own Edge operator indexes. To do that, you must define the TypeFlag member as EDGE_FILTER_USER.
Otherwise, if you want a predefined operator (sobel), you can set this member to 0 and use the EDGE_FILTER_SOBEL equate in the TypeFlag parameter.
SobelValB - Is the 'magic' numbers used in an edge operator. ex: In sobel, this value is 0-2, for scharr, this value is 162 and so on. It is the number related to the M4 position on the matrix index.
You can set whatever value you want to create your own Edge operator indexes. To do that, you must define the TypeFlag member as EDGE_FILTER_USER.
Otherwise, if you want a predefined operator (sobel), you can set this member to 0 and use the EDGE_FILTER_SOBEL equate in the TypeFlag parameter
TypeFlag - A set of predefined Flags to use for specific operators. Currently, you can use any of these 3 equates:
Equate Name Equate Value Meaning
EDGE_FILTER_USER 0 Use this flag, if you want to create your own pre-defined operator. You must fill the values on SobelValA and SobelValB
EDGE_FILTER_SOBEL 1 Use this flag to enable only the Sobel Operator. The values on SobelValA and SobelValB will be ignored (you can, as well, simply set them to zero. But for convention, better is use 0 insetad)
EDGE_FILTER_SCHARR 2 Use this flag to enable only the Sobel Operator. The values on SobelValA and SobelValB will be ignored (you can, as well, simply set them to zero. But for convention, better is use 0 insetad)
;;
; Equates used for the EdgeSobelPlus structure
[EdgeSobelPlus.MatrixValADis 0
EdgeSobelPlus.MatrixValBDis 4
EdgeSobelPlus.NormalizeRatioDis 8
EdgeSobelPlus.SobelGxDis 12
EdgeSobelPlus.SobelGyDis 16
EdgeSobelPlus.SobelGDis 20
EdgeSobelPlus.AngleDis 24
EdgeSobelPlus.SobelGxMinDis 28
EdgeSobelPlus.SobelGyMinDis 32
EdgeSobelPlus.SobelGMinDis 36
EdgeSobelPlus.SobelGxMaxDis 40
EdgeSobelPlus.SobelGyMaxDis 44
EdgeSobelPlus.SobelGMaxDis 48
EdgeSobelPlus.AngleMinDis 52
EdgeSobelPlus.AngleMaxDis 56]
[Size_of_EdgeSobelPlus 60]
; Equates Flags used to make easier the computation
[EDGE_FILTER_USER 0
EDGE_FILTER_SOBEL 1
EDGE_FILTER_SCHARR 2]
Proc PrepareSobelEdges:
Arguments @pSobelStruct, @SobelValA, @SobelValD, @TypeFlag
Local @MatrixValA, @MatrixValD
Uses edi
mov eax D@SobelValA | mov D@MatrixValA eax
mov eax D@SobelValD | mov D@MatrixValD eax
.If D@TypeFlag = EDGE_FILTER_SOBEL
mov D@MatrixValA 0-1
mov D@MatrixValD 0-2
.Else_If D@TypeFlag = EDGE_FILTER_SCHARR
mov D@MatrixValA 47
mov D@MatrixValD 161; 162 (official value). 161 is the best value i found because the sum of MatrixA+matrixB applying the formula below gives 255
.Else
; if both are zero, avoid division by zero on the normalization ratio. The ratio is, therefore: zero.
mov edi D@MatrixValA
.Test_If_Not_And eax eax, edi edi
mov edi D@pSobelStruct
fldz | fst F$edi+EdgeSobelPlus.MatrixValADis | fst F$edi+EdgeSobelPlus.MatrixValBDis | fstp F$edi+EdgeSobelPlus.NormalizeRatioDis
ExitP
.Test_End
.End_If
mov edi D@pSobelStruct ; All Ok, let´s calculate and put their values
fild D@MatrixValA | fst F$edi+EdgeSobelPlus.MatrixValADis | fabs | fmul R$Float_Two
fild D@MatrixValD | fst F$edi+EdgeSobelPlus.MatrixValBDis | fabs | faddp ST1 ST0
fmul R$Float_255 | fld1 | fdivrp ST0 ST1 | fstp F$edi+EdgeSobelPlus.NormalizeRatioDis
EndP
Example of the HDTV luminosity conversion with the 4 pixel boundary.
I see some black borders around it.
; Blue Green Red ( Alpha, unused = 0.0 )
Luma_PAL_SECAM_NTSC real4 0.1140, 0.5870, 0.2990, 0.0 ; rec601 luma (Y') component
Luma_HDTV real4 0.0722, 0.7152, 0.2126, 0.0 ; ITU-R BT.709 standard luma (Y') component
Luma_HDR real4 0.0593, 0.6780, 0.2627, 0.0 ; ITU-R BT.2100 standard luma (Y') component
Hypnotoad :biggrin: :biggrin: :biggrin: :biggrin: :greensml: :greensml: :greensml: :greensml:Retinex!?
(https://i.ibb.co/x8pNLMv/Fixing-Light-Default-User-Defined4.jpg) (https://ibb.co/x8pNLMv)
Hypnotoad :biggrin: :biggrin: :biggrin: :biggrin: :greensml: :greensml: :greensml: :greensml:Oops ended up in game forum,preview of the soon to be released "froggy 2020" :greenclp: :greenclp: :greenclp:
(https://i.ibb.co/x8pNLMv/Fixing-Light-Default-User-Defined4.jpg) (https://ibb.co/x8pNLMv)
Hypnotoad :biggrin: :biggrin: :biggrin: :biggrin: :greensml: :greensml: :greensml: :greensml:Oops ended up in game forum,preview of the soon to be released "froggy 2020" :greenclp: :greenclp: :greenclp:
(https://i.ibb.co/x8pNLMv/Fixing-Light-Default-User-Defined4.jpg) (https://ibb.co/x8pNLMv)
[-1 -2 -1 ]
SobelX = [ 0 0 0 ] = 0 degree
[+1 +2 +1 ]
[-1 0 +1 ]
SobelY = [-2 0 +2 ] = 90 degree
[-1 0 +1 ]
[ P1 P2 P3 ]
Pixels = [ P4 P5 P6 ]
[ P7 P8 P9 ]
Convolution with simplified math:
GradientX = (P7 + 2*P8 + P9) - (P1 + 2*P2 + P3) ; - if edge is on the left side, + if on the right side
GradientY = (P3 + 2*P6 + P9) - (P1 + 2*P4 + P7) ; - if edge is on the top side, + if on the bottom side
Convolution without multiplications:
GradientX = (P7 + P8 + P8 + P9) - (P1 + P2 + P2 + P3) ; - if edge is on the left side, + if on the right side
GradientY = (P3 + P6 + P6 + P9) - (P1 + P4 + P4 + P7) ; - if edge is on the top side, + if on the bottom side
Great :thumbsup: :thumbsup: :thumbsup: :thumbsup:
Also we need to normalize the result, right ?
I have not really kept up with this topic but it sounds like an algorithm that does sharpening and smoothing, depending on the settings. I use an ffmpeg option for doing this type of work on 1080 and 4k video.
Great :thumbsup: :thumbsup: :thumbsup: :thumbsup:
Also we need to normalize the result, right ?
Maybe, maybe not, that's the question. :biggrin:
Have to implement it first, and see if threshold clamping is an option.
Just finished a routine to save the gray color conversions to memory.
It handles 4 pixels at once and takes care of any possible "rest pixels".
It has a Stride for the width, wich is always a multiple of 16 bytes, no matter if the bitmap width is not divisible by 4.
Works like charm and is very fast.
I wrote it for SSE2 but it can be made faster for SSE4.1 if needed.
And this routine is a good candidate for multithreading if it is needed in realtime video processing.
Hi Hutch,
It's for edge detection in bitmaps.
guga want's to use it for removing watermarks.
I collaborate with him to learn a few new things.
15 25 99 100 in xmm0
movss xmm1,FLT4(1.5) ; move a single float in xmm1
shufps xmm1,xmm1,0 ; xmm1 = 1.5, 1.5, 1.5, 1.5
mulps xmm0,xmm1 ; xmm0 = 22.5, 37.5, 148.5, 150.0
subss xmm0,xmm1
addss xmm0,xmm1
subps xmm0,xmm1
addps xmm0,xmm1
Code: [Select]15 25 99 100 in xmm0
movss xmm1,FLT4(1.5) ; move a single float in xmm1
shufps xmm1,xmm1,0 ; xmm1 = 1.5, 1.5, 1.5, 1.5
mulps xmm0,xmm1 ; xmm0 = 22.5, 37.5, 148.5, 150.0
subtract and add single floats:Code: [Select]subss xmm0,xmm1
addss xmm0,xmm1
subtract and add 4 floats at once:Code: [Select]subps xmm0,xmm1
addps xmm0,xmm1
Hi Guga,
I think I understand what you are doing and an edge detection algo is a lot more useful than just removing water marks, it can also be used for both smoothing rough images or sharpening slightly blurry ones. What I suggested originally was a technique that I have been using for years, close proximity cloning to get the field depth correct and match the texture where possible.
(http://masm32.com/private/both.jpg)
The original on the left, the modified on the right. I used some Russian software to remove the buttons and manually cloned the background to get rid of the insignia. My comment is you can do both as both techniques have their advantages and between them the watermarks you want to remove can be filled in with a texture that maintains the correct field depth.
Hi Guga,
The app is called "Movavi Photo Editor 6". It does exactly what I did not need, tweaking pretty girls to look better but its fun to play with. I did the cloning with my very old and now buggy MicroGrafx Picture Publisher as it has exactly the right tools for doing this.
My suggestion was once you get the edge detection algo up and going, try sampling close to the watermark area to get the texture and field depth then selectively clone small parts to get rid of any blurring. I have seen robot watermark removal and some of it is real junk, a blurred mess that stands out very badly due to no compatibility of field depth.
:biggrin:Wow. Tks you, steve
> irfanview when i need to decrease the size without quality loss
I could lead you astray here, my JpgTool complete with sauce (whoops I mean source) does that just fine. With input from Marinus and Vortex the GDI+ code works really well.
BitmapData struct
dwWidth dd ?
dwHeight dd ?
Stride dd ?
PixelFormat dd ?
Scan0 dd ?
Reserved dd ?
BitmapData ends
.data?
align 4
GDIplusBitmapData BitmapData <?>
.data
align 4
SourceBitmapData dd NULL
.code
align 4
ReleaseBitmapMemory proc mem_ptr:DWORD
.if mem_ptr != NULL
invoke VirtualFree,mem_ptr,0,MEM_RELEASE
mov mem_ptr,NULL
.else
mov eax,1
.endif
ret
ReleaseBitmapMemory endp
align 4
LoadBitmapData proc uses ebx esi edi ImageName:DWORD,pImageData:DWORD
mov ReturnMessage,E_FAIL
invoke MultiByteToWideChar,CP_ACP,0,ImageName,-1,offset FilenameW,MAX_PATH-1
invoke GdipCreateBitmapFromFile,offset FilenameW,offset pImage
test eax,eax
jnz Close_Gdiplus
invoke GdipBitmapLockBits,pImage,NULL,ImageLockModeRead,PixelFormat32bppARGB,offset GDIplusBitmapData
test eax,eax
jnz Close_Image
mov esi,pImageData
invoke ReleaseBitmapMemory,dword ptr [esi]
test eax,eax
jz UnlockBitmap
mov eax,GDIplusBitmapData.dwWidth
shl eax,2
mul GDIplusBitmapData.dwHeight
invoke VirtualAlloc,0,eax,MEM_COMMIT or MEM_RESERVE,PAGE_READWRITE
test eax,eax
jz UnlockBitmap
mov dword ptr [esi],eax ; save the virtual memory pointer
mov edi,eax ; pointer to the raw 32 bit bitmap data
CopyBitMap:
mov esi,GDIplusBitmapData.Scan0 ; pointer to the bitmap data
mov ecx,GDIplusBitmapData.dwHeight
Height_lp:
mov edx,GDIplusBitmapData.dwWidth
xor ebx,ebx
Width_lp:
mov eax,dword ptr [esi+ebx]
mov dword ptr [edi],eax
add ebx,4
add edi,4
dec edx
jnz Width_lp
add esi,GDIplusBitmapData.Stride
dec ecx
jnz Height_lp
mov ReturnMessage,D3D_OK
UnlockBitmap:
invoke GdipBitmapUnlockBits,pImage,offset GDIplusBitmapData
Close_Image:
invoke GdipDisposeImage,pImage
Close_Gdiplus:
mov ecx,GDIplusBitmapData.dwWidth
mov edx,GDIplusBitmapData.dwHeight
Exit_LoadImage:
mov eax,ReturnMessage
ret
LoadBitmapData endp
; function:
invoke LoadBitmapData,TEXT_("Lena320.png"),offset SourceBitmapData
I_pBits = the memory pointer to the raw 32bit ARGB bitmap data in memory
I_Stride = I_Width * 4
invoke GdipCreateBitmapFromScan0,I_Width,I_Height,I_Stride,PixelFormat32bppRGB,I_pBits,offset pImage
This is how to create a pImage from the raw 32bit ARGB bitmap data in memoryCode: [Select]I_pBits = the memory pointer to the raw 32bit ARGB bitmap data in memory
I_Stride = I_Width * 4
invoke GdipCreateBitmapFromScan0,I_Width,I_Height,I_Stride,PixelFormat32bppRGB,I_pBits,offset pImage
Now you can use "pImage" with one of the routines in GdiPlusColors.zip to save it in the format you like.
DEFINE_GUID IMGFMT_WEBP, 0x0B96B3CB7,0x0728,0x11d3,0x9d,0x7b,0x00,0x00,0xf8,0x1e,0xf3,0x2e);
DEFINE_GUID IMGFMT_HEIF, 0x0B96B3CB6,0x0728,0x11d3,0x9d,0x7b,0x00,0x00,0xf8,0x1e,0xf3,0x2e);
GUID_WICPixelFormatDontCare dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0
GUID_WICPixelFormat1bppIndexed dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 1
GUID_WICPixelFormat2bppIndexed dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 2
GUID_WICPixelFormat4bppIndexed dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 3
GUID_WICPixelFormat8bppIndexed dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 4
GUID_WICPixelFormatBlackWhite dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 5
GUID_WICPixelFormat2bppGray dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 6
GUID_WICPixelFormat4bppGray dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 7
GUID_WICPixelFormat8bppGray dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 8
GUID_WICPixelFormat16bppBGR555 dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 9
GUID_WICPixelFormat16bppBGR565 dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0Ah
GUID_WICPixelFormat16bppBGRA5551 dd 5EC7C2Bh
dw 0F1E6h
dw 4961h
db 0ADh, 46h, 0E1h, 0CCh, 81h, 0Ah, 87h, 0D2h
GUID_WICPixelFormat16bppGray dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0Bh
GUID_WICPixelFormat24bppBGR dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0Ch
GUID_WICPixelFormat24bppRGB dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0Dh
GUID_WICPixelFormat32bppBGR dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0Eh
GUID_WICPixelFormat32bppBGRA dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 0Fh
GUID_WICPixelFormat32bppPBGRA dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 10h
GUID_WICPixelFormat32bppGrayFloat dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 11h
GUID_WICPixelFormat32bppRGBA dd 0F5C7AD2Dh
dw 6A8Dh
dw 43DDh
db 0A7h, 0A8h, 0A2h, 99h, 35h, 26h, 1Ah, 0E9h
GUID_WICPixelFormat32bppPRGBA dd 3CC4A650h
dw 0A527h
dw 4D37h
db 0A9h, 16h, 31h, 42h, 0C7h, 0EBh, 0EDh, 0BAh
GUID_WICPixelFormat48bppRGB dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 15h
GUID_WICPixelFormat48bppBGR dd 0E605A384h
dw 0B468h
dw 46CEh
db 0BBh, 2Eh, 36h, 0F1h, 80h, 0E6h, 43h, 13h
GUID_WICPixelFormat64bppRGBA dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 16h
GUID_WICPixelFormat64bppBGRA dd 1562FF7Ch
dw 0D352h
dw 46F9h
db 97h, 9Eh, 42h, 97h, 6Bh, 79h, 22h, 46h
GUID_WICPixelFormat64bppPRGBA dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 17h
GUID_WICPixelFormat64bppPBGRA dd 8C518E8Eh
dw 0A4ECh
dw 468Bh
db 0AEh, 70h, 0C9h, 0A3h, 5Ah, 9Ch, 55h, 30h
GUID_WICPixelFormat16bppGrayFixedPoint dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 13h
GUID_WICPixelFormat48bppBGRFixedPoint dd 49CA140Eh
dw 0CAB6h
dw 493Bh
db 9Dh, 0DFh, 60h, 18h, 7Ch, 37h, 53h, 2Ah
GUID_WICPixelFormat32bppCMYK dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 1Ch
GUID_WICPixelFormat64bppRGBAFixedPoint dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 1Dh
GUID_WICPixelFormat64bppBGRAFixedPoint dd 356DE33Ch
dw 54D2h
dw 4A23h
db 0BBh, 4, 9Bh, 7Bh, 0F9h, 0B1h, 0D4h, 2Dh
GUID_WICPixelFormat64bppRGBFixedPoint dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 40h
GUID_WICPixelFormat64bppRGBAHalf dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 3Ah
GUID_WICPixelFormat64bppRGBHalf dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 42h
GUID_WICPixelFormat16bppGrayHalf dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 3Eh
GUID_WICPixelFormat32bppGrayFixedPoint dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 3Fh
GUID_WICPixelFormat64bppCMYK dd 6FDDC324h
dw 4E03h
dw 4BFEh
db 0B1h, 85h, 3Dh, 77h, 76h, 8Dh, 0C9h, 1Fh
GugaImageCodecInfo struc
CLSID dd ?
FormatID dd ?
CodecName dd ?
FormatDescription dd ?
FilenameExtension dd ?
MimeType dd ?
Version dd 4 dup(?)
Signature dd ?
SignatureReserved dd ?
pInstanceFunction dd ?
GugaImageCodecInfo ends
Marinus, what´s the format of the structure"GDIplusBitmapData" ? I can´t find it .
BitmapData struct
dwWidth dd ?
dwHeight dd ?
Stride dd ?
PixelFormat dd ?
Scan0 dd ?
Reserved dd ?
BitmapData ends
GDIplusBitmapData BitmapData <?>
; https://www.autohotkey.com/boards/viewtopic.php?style=17&t=66335
; https://git.cs.usask.ca/SumatraPDF/SumatraPDF/blob/master/src/utils/WebpReader.cpp
; https://chromium.googlesource.com/webm/libwebp/+/master/src/webp/decode.h
; https://developers.google.com/speed/webp/docs/api
Proc DecodeWebpFile:
Arguments @pFileName
Local @hFile, @FileSize, @pImgWidth, @pImgHeight, @pBits, @TmpOpenedFile
lea ebx D@TmpOpenedFile ; Buffer to hold the loaded webp in memory
lea ecx D@FileSize ; and store it´s size
call ReadOpenedFile 0, D@pFileName, ebx, ecx ; A function to open the webp file and read it´s contents.
; Just to get some info of the web p image. It´s not needed since the width and height will already be retrieved with the function below
; lea ebx D@pImgWidth
; lea ecx D@pImgHeight
; C_call 'libwebp.WebPGetInfo' eax, D@FileSize, ebx, ecx
; This function gets the Pixel data from the webp image and also it´s width and height
lea ebx D@pImgWidth
lea ecx D@pImgHeight
C_call 'libwebp.WebPDecodeBGRA' D@TmpOpenedFile, D@FileSize, ebx, ecx ; libwebp library is in cdecl calling convention. The stack needs to be adjusted, that´s why i´m using RosAsm "C_call" macro that simply add XXX bytes to the stack after the call.
mov D@pBits eax ; Copy the image pixels to export to other formats
call CreateImageFromMemory D@pBits, D@pImgWidth, D@pImgHeight, pImage ; Marinus routine to create a pImage from the raw 32bit ARGB bitmap data in memory
call SaveNonIndexedImage D$pImage, {B$ "GugaWebP.png", 0}, Image_PNG, PixelFormat32bppRGB
call SaveJpgQualityImage D$pImage, {B$ "GugaWebP.jpg", 0}, 40
C_call 'libwebp.WebPFree' D@pBits ; Must release the data loaded from WebPDecodeBGRA
call 'RosMem.VMemFree' D@TmpOpenedFile ; and released our allocated memory
EndP
______________
;;
I_pBits = the memory pointer to the raw 32bit ARGB bitmap data in memory
I_Stride = I_Width * 4
;;
Proc CreateImageFromMemory:
Arguments @ImgBits, @ImageWidth, @ImageHeight, @pImage
Local @Stride
mov eax D@ImageWidth | shl eax 2
call 'gdiplus.GdipCreateBitmapFromScan0' D@ImageWidth, D@ImageHeight, eax, PIXELFORMAT_32BPPARGB, D@ImgBits, D@pImage
EndP
____________________
; Just a function to open the file onto a memory buffer
[NumberOfReadBytes: 0]
Proc ReadOpenedFile:
Arguments @Adresseem, @Filename, @pFileData, @pFileLenght
Local @hFile, @hFileSize, @MemFile
Uses edi, ecx, edx
call 'KERNEL32.CreateFileA' D@Filename, &GENERIC_READ,
&FILE_SHARE_READ+&FILE_SHARE_WRITE, 0, &OPEN_EXISTING,
&FILE_ATTRIBUTE_NORMAL, 0
If eax = &INVALID_HANDLE_VALUE
call 'USER32.MessageBoxA' 0, {B$ "File cannot be opened!", 0}, {'Error', 0}, &MB_OK__&MB_ICONWARNING__&MB_SYSTEMMODAL
xor eax eax | ExitP
End_If
mov D@hFile eax
call 'KERNEL32.GetFileSize' eax, 0
mov D@hFileSize eax
On eax = 0, ExitP
mov edi D@pFileLenght
add eax 10 | mov edx eax
mov D$edi edx
; Allocate enough memory
mov D@MemFile 0 | lea eax D@MemFile
call 'RosMem.VMemAlloc' eax, edx
mov edi D@pFileData
mov D$edi eax
call 'KERNEL32.ReadFile' D@hFile, eax,
D@hFileSize, NumberOfReadBytes, 0
call 'KERNEL32.CloseHandle' D@hFile
EndP
webpUrl := "https://upload.wikimedia.org/wikipedia/commons/b/b2/Vulphere_WebP_OTAGROOVE_demonstration_2.webp"
libwebp32Url := "https://s3.amazonaws.com/resizer-dynamic-downloads/webp/0.5.2/x86/libwebp.dll"
libwebp64Url := "https://s3.amazonaws.com/resizer-dynamic-downloads/webp/0.5.2/x86_64/libwebp.dll"
filePath := A_ScriptDir . "\test.webp"
bitness := A_PtrSize*8
lib := A_ScriptDir . "\libwebp" . bitness . ".dll"
if !FileExist(filePath)
URLDownloadToFile, % webpUrl, % filePath
if !FileExist(lib)
URLDownloadToFile, % libwebp%bitness%Url, % lib
hBitmap := HBitmapFromWebP(lib, filePath, width, height)
Gui, Margin, 0, 0
Gui, Add, Pic, w600 h-1, HBITMAP:%hBitmap%
Gui, Show
Return
GuiClose() {
ExitApp
}
HBitmapFromWebP(libwebp, WebpFilePath, ByRef width, ByRef height) {
file := FileOpen(WebpFilePath, "r")
len := file.RawRead(buff, file.Length)
file.Close()
if !len
throw Exception("Failed to read the image file")
if !hLib := DllCall("LoadLibrary", Str, libwebp, Ptr)
throw Exception("Failed to load library. Error:" . A_LastError)
if !pBits := DllCall(libwebp . "\WebPDecodeBGRA", Ptr, &buff, Ptr, len, IntP, width, IntP, height) {
DllCall("FreeLibrary", Ptr, hLib)
throw Exception("Failed to decode the image file")
}
oGDIp := new GDIp
pBitmap := oGDIp.CreateBitmapFromScan0(width, height, pBits)
hBitmap := oGDIp.CreateHBITMAPFromBitmap(pBitmap)
DllCall(libwebp . "\WebPFree", Ptr, pBits)
oGDIp.DisposeImage(pBitmap)
DllCall("FreeLibrary", Ptr, hLib)
Return hBitmap
}
class GDIp {
__New() {
if !DllCall("GetModuleHandle", Str, "gdiplus", Ptr)
DllCall("LoadLibrary", Str, "gdiplus")
VarSetCapacity(si, A_PtrSize = 8 ? 24 : 16, 0), si := Chr(1)
DllCall("gdiplus\GdiplusStartup", UPtrP, pToken, Ptr, &si, Ptr, 0)
this.token := pToken
}
__Delete() {
DllCall("gdiplus\GdiplusShutdown", Ptr, this.token)
if hModule := DllCall("GetModuleHandle", Str, "gdiplus", Ptr)
DllCall("FreeLibrary", Ptr, hModule)
}
CreateBitmapFromScan0(Width, Height, pBits, PixelFormat := 0x26200A, Stride := "") {
if !Stride {
bpp := (PixelFormat & 0xFF00) >> 8
Stride := ((Width * bpp + 31) & ~31) >> 3
}
DllCall("gdiplus\GdipCreateBitmapFromScan0", Int, Width, Int, Height
, Int, Stride, Int, PixelFormat
, Ptr, pBits, PtrP, pBitmap)
Return pBitmap
}
CreateHBITMAPFromBitmap(pBitmap, Background=0xffffffff) {
DllCall("gdiplus\GdipCreateHBITMAPFromBitmap", Ptr, pBitmap, PtrP, hbm, Int, Background)
return hbm
}
DisposeImage(pBitmap) {
return DllCall("gdiplus\GdipDisposeImage", Ptr, pBitmap)
}
}
/* Copyright 2015 the SumatraPDF project authors (see AUTHORS file).
License: Simplified BSD (see COPYING.BSD) */
#include "BaseUtil.h"
#include "WebpReader.h"
#ifndef NO_LIBWEBP
#include <webp/decode.h>
using namespace Gdiplus;
namespace webp {
// checks whether this could be data for a WebP image
bool HasSignature(const char *data, size_t len)
{
return len > 12 && str::StartsWith(data, "RIFF") && str::StartsWith(data + 8, "WEBP");
}
Size SizeFromData(const char *data, size_t len)
{
Size size;
WebPGetInfo((const uint8_t *)data, len, &size.Width, &size.Height);
return size;
}
Bitmap *ImageFromData(const char *data, size_t len)
{
int w, h;
if (!WebPGetInfo((const uint8_t *)data, len, &w, &h))
return nullptr;
Bitmap bmp(w, h, PixelFormat32bppARGB);
Rect bmpRect(0, 0, w, h);
BitmapData bmpData;
Status ok = bmp.LockBits(&bmpRect, ImageLockModeWrite, PixelFormat32bppARGB, &bmpData);
if (ok != Ok)
return nullptr;
if (!WebPDecodeBGRAInto((const uint8_t *)data, len, (uint8_t *)bmpData.Scan0, bmpData.Stride * h, bmpData.Stride))
return nullptr;
bmp.UnlockBits(&bmpData);
// hack to avoid the use of ::new (because there won't be a corresponding ::delete)
return bmp.Clone(0, 0, w, h, PixelFormat32bppARGB);
}
}
#else
namespace webp {
bool HasSignature(const char *data, size_t len) { UNUSED(data); UNUSED(len); return false; }
Gdiplus::Size SizeFromData(const char *data, size_t len) { UNUSED(data); UNUSED(len); return Gdiplus::Size(); }
Gdiplus::Bitmap *ImageFromData(const char *data, size_t len) { UNUSED(data); UNUSED(len); return nullptr; }
}
#endif
Marinus, if you would like to add more equates to GdiPlus header. here are some more:
WebP images will not load with GdiPlus
Maybe you need to install a decoder/encoder for GdiPlus???
I say this because you provided the "IMAGE_WEBP 9" Image Type.
Hi guga,
We have to calculate with square roots for the magnitudes.
How much precision do we need?
We can calculate with high precision at the cost of +/- 39 clock cycles.
Or we can calculate the reciprocal square roots with 12 bit precision at +/- 8 clock cycles.
Or do it the rude way ( not very accurate ) and make the X Y gradients absolute values and add them together at +/- 3 clock cycles.
I myself would go for the reciprocal square roots, it's more than enough precision we need for the pixel bytes, which are actually rounded to natural numbers.
Just take a moment rest, then I have a chance to keep up with you. :biggrin:
Log and Exp routines are not available in SSE, we have to create them ourselves.
Did a test with the machine image from wikipedia's Sobel example.
The principal of the 16 pixel sse2 parallel sobel routine works, but there is some fine tuning needed.
There are misalignments on the Y axis, as you can see in the attachment.
Cool.learn from RCP** instructions the backward math:way of faster coding together with replace x DIV** coefficients to x MUL** 1/coeffcients with upto REAL8 precisicion
So, you're in love now with SSE... :bgrin:
Maybe it is possible to transfer it into a matrix with pre-calculated coefficients???