Fast median algorithm

[Siekmanski]

Maby we could use this as the Laplace edge detection operator????

[ 0  -1  0 ]
[-1   4 -1 ]
[ 0  -1  0 ]

Not sure though.......

[guga]

Hi Guys....about the fast median algo, after a loong delay....i found a small bug.

It is not calculating the median when we have 0 as a valid valid. I mean, suppose we have a sequence of 0, 0,2,2, 0, 2   or 0, 0,1,1, 0, 1. I both cases, the median value should result in 0.5 or 1 (in the 2nd example). Insetad it keeps resulting 0

i'm trying to fix the algo, but couldn´t do it.

The original algo Siekmanski did was:

Code Select Expand

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

I did a small variation (In RosAsm syntax) where it is used to take a byte at a time on a sequence of pixels. For example. Imagine we have a sequence of 6 images and i want to calculate the median of the 1st byte on each one of them. Suppose that each image has only 3 pixels (Therefore, 12 bytes long (3*4, because 4 is the amount of bytes used in RGBA pixel data). Like this:

[Image1: Pixel1, Pixel2, Pixel3
Image2: Pixel1, Pixel2, Pixel3
Image3: Pixel1, Pixel2, Pixel3
Image4: Pixel1, Pixel2, Pixel3
Image5: Pixel1, Pixel2, Pixel3
Image6: Pixel1, Pixel2, Pixel3]

So, i want to calculate the median of Pixel 1 from Image1 to Image6. So suppose they have these values:
Pixel = RGA, and i want the Median of Red channel, for example

[Image1: B$ 0, 158, 198, 255,   43, 158, 198, 255,   4, 158, 198, 255
Image2: B$ 0, 13, 255, 255,   14, 13, 255, 255,   14, 13, 255, 255
Image3: B$ 1, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255
Image4: B$ 1, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255
Image5: B$ 0, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255
Image6: B$ 1, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255]

So, i want to calculate the median on the 1st byte on all of them (in red). So i have a sequence of
0, 0, 1, 1, 0, 1

The correct result should be a median of 0.5

I converted Siekmanski algorithm to output the results on a Real8 variable. (Remembering hat the input is a byte (integer)

Code Select Expand

Macros used:
[Do | D0:]
[Repeat_Until_Zero | dec #1 | jne D0<]

[Loop_Until | #=3 | cmp #1 #3 | jn#2 D0<]

[Test_Do | Q6:]
[Test_Until | #=2 | test #1 #2 | jz Q6<]



Code Select Expand


;;
Input - Pointer to the start of a sequence of images
ChunckSize - The number of images used to calculate the median at each position in the input
SkipPos  - The size of each image (In bytes)
pOutput - Pointer to a Real8 variable to hold the median value from the sequence of pixels at XX position in Input.

Remarks: Since "Input" refers to the starting address of each pixel, to calculate the median of the Red channel of pixel1 in all images, all we need to do is set as input the 1st byte (red) of pixel 1. Example. Say we use edi as input

mov edi nput
Pixel1 Red channel = 1st byte in edi
Pixel1 Green channel = 2nd byte in edi . So to calculate the median of all Green in Image1 to Image6 we need only to add 1 to edi. So, it will point exactly to the green byte related to the 1st pixel
Pixel1 Blue channel = 3rd byte in edi So to calculate the median of all Blue in Image1 to Image6 we need only to add 1 to edi. So, it will point exactly to the green byte related to the 1st pixel


Example of Usage:


[GugaOutput: R$ 0]

[testingNew: B$ 0, 158, 198, 255,   43, 158, 198, 255,   4, 158, 198, 255
testingNew2: B$ 0, 13, 255, 255,   14, 13, 255, 255,   14, 13, 255, 255
testingNew3: B$ 1, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255
testingNew4: B$ 1, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255
testingNew5: B$ 0, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255
testingNew6: B$ 1, 83, 5, 255,     19, 83, 5, 255,     19, 83, 5, 255]

to calculate the median of Red data in the 6 images located all in the 1st Pixel
mov esi testingNew
call CalculateMedianByte esi, 6, 12, GugaOutput

to calculate the media of green we simply do
mov esi testingNew | inc esi ; advance esi only 1 byte to we reach always the Blue pixel of the desired position on all images.
call CalculateMedianByte esi, 6, 12, GugaOutput

to calculate the red of the 2nd pixel we simply add 4 bytes to esi. to it reach the related Red Byte on the 2nd pixel of image 1 (and therefore, will also reach the same byte from Image2 to Image6 - as long all of the images are in sequence and has the same size)


;;

[MedianBuffer: D$ 0 #(256*2)] ; minimum buffer size must be as large as the biggest number in the data set
[Float_Half: R$ 0.5]

Proc CalculateMedianByte:
    Arguments @Input, @ChunckSize, @SkipPos, @pOutput
    Local @Result, @Ishalf
    Uses ebx, esi, edi, ecx, edx

    call 'RosMem.FastZeroMem' MedianBuffer, (256*4)
    mov D@Ishalf 0
    mov esi D@Input
    mov edi MedianBuffer ; Zero the MedianBuffer before using this routine.
    mov ecx D@ChunckSize

    ; PopulateLP Diagonally from left to right
    ; if skippos = 0, we increment data not from 4 to 4, but from the size of each image. So the next byte to collect correspond to the next image)
    mov ebx D@SkipPos
    If ebx = 0
        inc ebx
    End_If

    Do
        movzx eax B$esi
        inc D$edi+eax*4   ; Populate the Median buffer with the values as index numbers
        add esi ebx       ; and keep track for duplicates in one go. And skip positions if needed
    Repeat_Until_Zero ecx

    xor esi esi         ;  Uneven marker
    mov ecx D@ChunckSize
    mov ebx ecx
    shr ebx 1           ; The middle of the PopulationSize.
    and ecx 1 | jne @UnEven ; Check even or uneven.
    dec ebx             ; Make even
    inc esi             ; Set marker true = even

; UnEven
@UnEven:

    xor edx edx     ; Median counter
    mov eax 0-1     ; Number counter.

    ; FindMedianLP
    Do
        inc eax
        mov ecx D$edi+eax*4 ; Count of single and/or duplicate numbers. ( zero is not a number )
        add edx ecx         ; Add it to the Median counter.
    Loop_Until edx > ebx    ; Compare if we reached the middle of the PopulationSize and Repeat until we reached the middle.

    dec esi | js @Ready        ; Check even or uneven. Ready if uneven

    ; Check if next number is a duplicate. ( no averaging needed )
    ; Ready if next number is the same value as the previous number
    .If ecx <= 1
        mov edx eax     ; Save the first middle number to calculate the average of 2 middle numbers
        ; here we are always at the end of the left middle. Ex if size = 100, we are at pos 49th
        ; FindSecondMiddleNumberLP. Funciona uando o numero for par
        Test_Do
            inc eax
            mov ecx D$edi+eax*4
        Test_Until ecx ecx
        add eax edx
        mov D@Ishalf 1
        ;shr eax 1           ; And get the average. No longer needed since the result is in Real8)
    .End_If

@Ready:

    mov edi D@pOutput
    mov D@Result eax
    fild D@Result
    If D@Ishalf = 1 ; And get the average.
        fmul R$Float_Half
    End_If
    fstp R$edi

EndP

So, 1st problem

How to fix it ?


2nd Problem

Once it is fixed how do we calculate also the MAD (median absolute distance) ? The median absolute deviation is given as:
https://www.statskingdom.com/median-absolute-deviation-calculator.html

So, once we calculate the median and taking onto account that the result must be a Real8 (for both median and MAD), how to calculate MAD faster too ?

The function cab be extended to reuse the table of integers as in "FindMedianLP" and later making it export a seconds outpu as the MAD. SO the funciotion can contains 2 outputs. One for median and other for MAD