Image color quantization procedure (FAD)

[El_Choni]

Code updated for 5.20+

Hi,

I've been coding for some time a ColorDepth lib to change image color depth, including dither. Currently, I have this quantization procedure which uses Fast Adaptive Dissection (an enhancement of Median Cut) to quantize colours. But I won't be able to release it until a month or two, too much work going on, so I've decided to post the PB source so everybody can benefit and report problems or suggestions. Speed and memory usage are not optimized yet, but the procedure seems to work. The usage of the procedure is simple:

Code: Select all

SetDepth(ImageID, PixelDepth)

So, if you have a 24 bit image and you want it to have a maximum of 256 colours, you would call it this way:

Code: Select all

Structure BoundingBox
  Error.f
  Low.w
  Up.w
  Direction.b
EndStructure

#UsedDepth = 31

Procedure SetDepth(ImageID, PixelDepth)
  result = 0
  GetObject_(ImageID, SizeOf(BITMAP), bm.BITMAP)
  bmih.BITMAPINFOHEADER
  bmih\biSize = SizeOf(BITMAPINFOHEADER)
  bmih\biWidth = bm\bmWidth
  bmih\biHeight = bm\bmHeight
  bmih\biPlanes = 1
  bmih\biBitCount = 32
  bmih\biCompression = #BI_RGB
  ImageBits = AllocateMemory(0, bmih\biWidth*bmih\biHeight*4)
  WindowID = WindowID(0)
  hDC = GetDC_(WindowID)
  GetDIBits_(hDC, ImageID, 0, bmih\biHeight, ImageBits, bmih, 0)
  ReleaseDC_(WindowID, hDC)
  BMPSeeker = ImageBits
  MaxColours = Pow(2, PixelDepth)
  Dim Histogram($8000-1) ; 5/6/4? Otra combinacion? Basada en diferencias en la imagen?
  For i=0 To bmih\biHeight-1                        ; Tener en cuenta al partir
    For t=0 To bmih\biWidth-1
      RGBColour = PeekL(BMPSeeker)
      r = (RGBColour&$F8)<<7
      g = (RGBColour&$F800)>>6
      b = (RGBColour&$F80000)>>19
      Histogram(r|g|b)+1
      BMPSeeker+4
    Next t
  Next i
  HColours = $8000-1
  Dim Colours.w(HColours-1, 2)
  i = 0
  For r=0 To #UsedDepth
    For g=0 To #UsedDepth
      For b=0 To #UsedDepth
        Colour = ((r<<10)|(g<<5)|b)&$7FFF
        If Histogram(Colour)
          Colours(i, 0) = Colour
          i+1
        EndIf
      Next b
    Next g
  Next r
  ImageColours = i-1
  If ImageColours>MaxColours-1
    i = 0
    For g=0 To #UsedDepth
      For b=0 To #UsedDepth
        For r=0 To #UsedDepth
          Colour = ((r<<10)|(g<<5)|b)&$7FFF
          If Histogram(Colour)
            Colours(i, 1) = Colour
            i+1
          EndIf
        Next r
      Next b
    Next g
    i = 0
    For b=0 To #UsedDepth
      For r=0 To #UsedDepth
        For g=0 To #UsedDepth
          Colour = ((r<<10)|(g<<5)|b)&$7FFF
          If Histogram(Colour)
            Colours(i, 2) = Colour
            i+1
          EndIf
        Next g
      Next r
    Next b
    Dim Region.BoundingBox(MaxColours-1)
    CurrentRegion = 0
    Region(CurrentRegion)\Low = 0
    Region(CurrentRegion)\Up = ImageColours
    LastHigh = 0
    LastLow = 0
    CurrentList = 0
    Lastdiff.f = 0
    For i=0 To 2
      Low = (Colours(Region(CurrentRegion)\Low, i)>>((2-i)*5))&$1F
      High = (Colours(Region(CurrentRegion)\Up, i)>>((2-i)*5))&$1F
      Select i
        Case 0
          diff.f = (High-Low)*0.59
        Case 1
          diff.f = (High-Low)*0.30
        Case 2
          diff.f = (High-Low)*0.11
      EndSelect
      If diff>Lastdiff
        CurrentList = i
        Lastdiff = diff
      EndIf
    Next i
    NewRegion = 1
    Region(CurrentRegion)\Direction = %11
    Dim leftvec.f(ImageColours)
    Dim rightvec.f(ImageColours)
    While NewRegion<MaxColours
      L = Region(CurrentRegion)\Low
      U = Region(CurrentRegion)\Up
      If Region(CurrentRegion)\Direction&%1
        LowR = #UsedDepth
        LowG = #UsedDepth
        LowB = #UsedDepth
        HighR = 0
        HighG = 0
        HighB = 0
        Popularity = 0
        For i=L To U-1
          Colour = Colours(i, CurrentList)
          r = (Colour&$7C00)>>10
          g = (Colour&$3E0)>>5
          b = Colour&$1F
          If r<LowR
            LowR = r
          EndIf
          If r>HighR
            HighR = r
          EndIf
          If g<LowG
            LowG = g
          EndIf
          If g>HighG
            HighG = g
          EndIf
          If b<LowB
            LowB = b
          EndIf
          If b>HighB
            HighB = b
          EndIf
          Popularity+Histogram(Colour)
          SoR.f = (HighR-LowR)*0.59
          SoG.f = (HighG-LowG)*0.30
          SoB.f = (HighB-LowB)*0.11
          leftvec(i) = Popularity*(Sqr(Pow(SoR, 2)+Pow(SoG, 2))+Pow(SoB, 2))
        Next i
      EndIf
      If Region(CurrentRegion)\Direction&%10
        LowR = #UsedDepth
        LowG = #UsedDepth
        LowB = #UsedDepth
        HighR = 0
        HighG = 0
        HighB = 0
        Popularity = 0
        For i=U To L+1 Step -1
          Colour = Colours(i, CurrentList)
          r = (Colour&$7C00)>>10
          g = (Colour&$3E0)>>5
          b = Colour&$1F
          If r<LowR
            LowR = r
          EndIf
          If r>HighR
            HighR = r
          EndIf
          If g<LowG
            LowG = g
          EndIf
          If g>HighG
            HighG = g
          EndIf
          If b<LowB
            LowB = b
          EndIf
          If b>HighB
            HighB = b
          EndIf
          Popularity+Histogram(Colour)
          SoR.f = (HighR-LowR)*0.59
          SoG.f = (HighG-LowG)*0.30
          SoB.f = (HighB-LowB)*0.11
          rightvec(i) = Popularity*(Sqr(Pow(SoR, 2)+Pow(SoG, 2))+Pow(SoB, 2))
        Next i
      EndIf
      PreviousError = leftvec(L)+rightvec(L+1)
      SplitPlace = L
      For i=L To U-1
        If leftvec(i)+rightvec(i+1)<PreviousError
          PreviousError = leftvec(i)+rightvec(i+1)
          SplitPlace = i
        EndIf
      Next i
      If SplitPlace>L
        Region(CurrentRegion)\Error = leftvec(SplitPlace)
      Else
        Region(CurrentRegion)\Error = 0
      EndIf
      Region(CurrentRegion)\Up = SplitPlace
      Region(CurrentRegion)\Direction = %10
      If SplitPlace<U-1
        Region(NewRegion)\Error = rightvec(SplitPlace+1)
      EndIf
      Region(NewRegion)\Low = SplitPlace+1
      Region(NewRegion)\Up = U
      Region(NewRegion)\Direction = %1
      LowC = 0
      HighC = 0
      For i=0 To 2
        ThisL = CurrentList+i
        If ThisL>2:ThisL-3:EndIf
        LowC = (LowC<<5)|((Colours(L, CurrentList)>>((2-ThisL)*5))&$1F)
        HighC = (HighC<<5)|((Colours(SplitPlace, CurrentList)>>((2-ThisL)*5))&$1F)
      Next i
      Dim belowList.w(U-L)
      Dim aboveList.w(U-L)
      List2 = CurrentList+1
      If List2>2:List2-3:EndIf
      List3 = List2+1
      If List3>2:List3-3:EndIf
      For d=1 To 2
        aListIndex = 0
        bListIndex = 0
        RCLI = CurrentList+d
        If RCLI>2:RCLI-3:EndIf
        For i=L To U
          OutBounds = 1
          Colour = Colours(i, RCLI)
          PColour = 0
          For z=0 To 2
            ThisL = CurrentList+z
            If ThisL>2:ThisL-3:EndIf
            PColour = (PColour<<5)|((Colour>>((2-ThisL)*5))&$1F)
          Next z
          If PColour>=LowC And PColour<=HighC
            OutBounds = 0
          EndIf
          If OutBounds
            aboveList(aListIndex) = Colour
            aListIndex+1
          Else
            belowList(bListIndex) = Colour
            bListIndex+1
          EndIf
        Next i
        For i=0 To bListIndex-1
          Colours(i+L, RCLI) = belowList(i)
        Next i
        For i=0 To aListIndex-1
          Colours(i+L+bListIndex, RCLI) = aboveList(i)
        Next i
      Next d
      For i=0 To NewRegion
        If Region(i)\Error>Region(CurrentRegion)\Error
          CurrentRegion = i
        EndIf
      Next i
      Lastdiff.f = 0
      For i=0 To 2
        Low = (Colours(Region(CurrentRegion)\Low, i)>>((2-i)*5))&$1F
        High = (Colours(Region(CurrentRegion)\Up, i)>>((2-i)*5))&$1F
        Select i
          Case 0
            diff.f = (High-Low)*0.59
          Case 1
            diff.f = (High-Low)*0.30
          Case 2
            diff.f = (High-Low)*0.11
        EndSelect
        If diff>Lastdiff
          CurrentList = i
          Lastdiff = diff
        EndIf
      Next i
      NewRegion+1
    Wend
    Dim leftvec.f(0)
    Dim rightvec.f(0)
    Dim belowList.w(0)
    Dim aboveList.w(0)
    Dim Palette(MaxColours-1)
    PaletteIndex = 0
    While PaletteIndex<MaxColours
      AverageR.f = 0
      AverageG.f = 0
      AverageB.f = 0
      ColoursNumber = 0
      For i=Region(PaletteIndex)\Low To Region(PaletteIndex)\Up
        Colour = Colours(i, 0)&$7FFF
        r = (Colour&$7C00)>>10
        g = (Colour&$3E0)>>5
        b = Colour&$1F
        AverageR+(r*Histogram(Colour))
        AverageG+(g*Histogram(Colour))
        AverageB+(b*Histogram(Colour))
        ColoursNumber+Histogram(Colour)
        Histogram(Colour) = PaletteIndex
      Next i
      If ColoursNumber
        AvR = Int(AverageR*8/ColoursNumber)
        AvG = Int(AverageG*8/ColoursNumber)
        AvB = Int(AverageB*8/ColoursNumber)
        Palette(PaletteIndex) = ((AvB&$FF)<<16)|((AvG&$FF)<<8)|(AvR&$FF)
      EndIf
      PaletteIndex+1
    Wend
    Dim Colours.w(0, 0)
    Dim Region.BoundingBox(0)
    BMPSeeker = ImageBits
    For i=0 To bmih\biHeight-1
      For t=0 To bmih\biWidth-1
        RGBColour = PeekL(BMPSeeker)
        r = (RGBColour&$F8)<<7
        g = (RGBColour&$F800)>>6
        b = (RGBColour&$F80000)>>19
        ; Dither aquí
        PokeL(BMPSeeker, Palette(Histogram(r|g|b)))
        BMPSeeker+4
      Next t
    Next i
    Dim Histogram(0)
    hDC = GetDC_(WindowID)
    SetDIBits_(hDC, ImageID, 0, bmih\biHeight, ImageBits, bmih, 0)
    ReleaseDC_(WindowID, hDC)
    RedrawWindow_(WindowID, 0, 0, 7)
    result = 1
  ElseIf ImageColours=MaxColours
    BMPSeeker = ImageBits
    For i=0 To bmih\biHeight-1
      For t=0 To bmih\biWidth-1
        RGBColour = PeekL(BMPSeeker)
        RGBColour!%1110000011100000111
        ; Dither aquí
        PokeL(BMPSeeker, RGBColour)
        BMPSeeker+4
      Next t
    Next i
    hDC = GetDC_(WindowID)
    SetDIBits_(hDC, ImageID, 0, bmih\biHeight, ImageBits, bmih, 0)
    ReleaseDC_(WindowID, hDC)
    RedrawWindow_(WindowID, 0, 0, 7)
    result = 1
  EndIf
  FreeMemory(0)
  ProcedureReturn result
EndProcedure

hWindow = OpenWindow(0, 0, 0, 320, 256, "Load picture example", #PB_Window_SystemMenu|#PB_Window_MinimizeGadget|#PB_Window_ScreenCentered|#PB_Window_Invisible)
If hWindow
  FileName$ = OpenFileRequester("Open image", "", "All supported formats|*.bmp;*.ico;*.cur|BMP image (*.bmp)|*.bmp|Icon file (*.ico)|*.ico|Cursor file (*.cur)|*.cur", 0)
  If FileName$
    CurrentDir$ = GetPathPart(FileName$)
    ImageID = LoadImage(0, FileName$)
    
    If ImageID
      width = ImageWidth(0)
      height = ImageHeight(0)
      ButtonImageGadget(0, 0, 0, width, height, ImageID)
      ResizeWindow(0, #PB_Ignore, #PB_Ignore, width, height)
      HideWindow(0, 0)
      AddKeyboardShortcut(0, #PB_Shortcut_Space, 0)
      SetForegroundWindow_(hWindow)
      Repeat
        EventID = WaitWindowEvent()
        If EventID=#PB_Event_Gadget
          FileName$ = OpenFileRequester("Open image", "", "All supported formats|*.bmp;*.ico;*.cur|BMP image (*.bmp)|*.bmp|Icon file (*.ico)|*.ico|Cursor file (*.cur)|*.cur", 0)
          If FileName$
            CurrentDir$ = GetPathPart(FileName$)
            FreeImage(0)
            ImageID = LoadImage(0, FileName$)
            If ImageID
              width = ImageWidth(0)
              height = ImageHeight(0)
              ResizeWindow(0, #PB_Ignore, #PB_Ignore, width, height)
              ResizeGadget(0, 0, 0, width, height)
              SetGadgetState(0, ImageID)
            EndIf
          EndIf
        ElseIf EventID=#PB_Event_Menu
          SetDepth(ImageID(0), 8)
        EndIf
      Until EventID=#PB_Event_CloseWindow
    EndIf
  EndIf
EndIf
End

[Hi-Toro]

Good stuff -- very fast too! Thanks for sharing it

[El_Choni]

Thank you for the feedbak

[Fred]

Just tested it and it works great..

[blueb]

Funny... :roll:

When I try to run this program, I get an error message on line #212

Code: Select all

        LowC = (LowC<<5)|((Colours(L, CurrentList)>>((2-ThisL)*5))&$1F) 

that states: Too complex expression.. please split it!

Running XP Pro SP1

Blueb

[Fred]

You got the Beta1 ?

[El_Choni]

I have no problems with 3.70 (not beta), but tf it says 'split it', then split it

Code: Select all

tmp = (Colours(L, CurrentList)>>((2-ThisL)*5))&$1F
LowC = (LowC<<5)|tmp

[LCD]

Any chance to add error diffusion? I tried myself, but if works perfectly only with Grey scale images (or splited RGB Channels, which cause grey to be transformed into random coloured pixels). Doing a error diffusion on quantized palette is still a miracle for me.

[El_Choni]

I'm working on that, but my first attempts look awful (much more awful than the non-dithered version). I have several dithering procedures that work very well when translating 24 bit to monochrome, so I guess I'll be able to adapt them. Some day (if you look at the source, you'll see that the places to do the dithering are marked).

[LCD]

Well, mine 24 bit to monochrome (B/W and grey scale) works well, but colour.... I even tried to work directly on 24 bit space in place of splited RGB channels, but this did not work too. And yes, I saw the dithering remark... Anyway, the procedure does a excellent job.

[El_Choni]

Maybe sharing our ideas we can come to a solution. Let's take Floyd-Steinberg as an example. As you know, this is the FS pattern:

Code: Select all

  X 7
3 5 1

Let's say we have a 256 colour quantized palette from a 24 bit image. First, we look at the corresponding palette colour for the current pixel, we substitute it and keep separated error values for R, G and B.

Now we add the weighted error to the R, G and B pixels that will hold the error diffusion.

And now we go to the next pixel; the problem now is that this pixel has been modified and is not directly associated to a colour in the palette, so we must find the closest match in the palette. I do it this way:

Code: Select all

CurrentColour = PeekL(CurrentPixel)
r = CurrentColour&$FF
g = (CurrentColour&$FF00)>>8
b = (CurrentColour&$FF0000)>>16
LastError.f = $FFFFFF
For m=0 to 255
  rP = Palette(m)&$FF
  gP = (Palette(m)&$FF00)>>8
  bP = (Palette(m)&$FF0000)>>16
   ; this calculates the distance between two colours in the RGB space
   ; there should be another Sqr around the whole expression, but this is Ok for
   ; doing comparisons, like in this case
  Error.f = Sqr(Pow(r-rP, 2)+Pow(g-gP, 2))+Pow(b-bP, 2)
  If Error<LastError
    LastError = Error
    ClosestMatch = m
  EndIf
Next m
rP = Palette(ClosestMatch)&$FF
gP = (Palette(ClosestMatch)&$FF00)>>8
bP = (Palette(ClosestMatch)&$FF0000)>>16
ErrorR = r-rP
ErrorG = g-gP
ErrorB = b-bP
PokeL(CurrentPixel, Palette(ClosestMatch))
; And do error diffusion

I think this should work, but the results are awful, as I said. Not exactly incorrect, just ugly. Any ideas you can share about this?

[LCD]

I´m at work now, but then I will check the new code, maybe I can find something...
Edit: Should add: x/16 to the floyd steinberg matrix
And Pow() requires two parameters if I remember it correctly.

[El_Choni]

Corrected Pow, thanks. Yes, I assume that every value in the matrix is divided by the sum of all the values, 16 in the case of Floyd-Steinberg.

[LCD]

Yeah, that´s better... Better than my results, because my try with own code gave me even much stranger results. Your code can be improved, I believe.

I pointed the value 16 in FS, because I write the matrixes in this way:
x/16 *7
*3 *5 *1

Or simple and fast error diffusion:
x/2 *1
*1

And there are other matrixes like JaJuNi, Burkes, and other with different values. But that´s right, I just need to add all matrix values together to get the division value. On the other side I assume that some other coders do not know what this all means, but they want to help to improve the code .

[LCD]

Now after some experiments I think, you are on the right way. The method you are using, is the pytagoras distance calculation, and I did exactly the same (okay, I optimised the speed a bit)

Code: Select all

      r1=err(x,y)
      g1=erg(x,y)
      b1=erb(x,y)
      minall=16777215
      col=0
      For a=0 To maxcolors
    dist=Sqr(((palr(a)-r1)*(palr(a)-r1))+((palg(a)-g1)*(palg(a)-g1))+((palb(a)-b1)*(palb(a)-b1)))
        If dist<minall
          col=a
          minall=dist
        EndIf
      Next a
      scr(x,y)=col
      error diffusion here...

it works great, so I don't know why your results should be bad. There is even a faster way of distance calculation, which sometimes works better, just one line should be replaced:

Code: Select all

dist=Abs(palr(a)-r1)+Abs(palg(a)-g1)+Abs(palb(a)-b1)

Maybe your Palette quantisation was not perfect... The comments are not in my language, and I'm not sure how do you compute the histogram, do you use Octree? Octree is one of the easiest methods to create a suitable palette for color reduction from 24 but images, without the need to sort a histogram field of 16777216 entrys.