Image posterization and dithering (with optional correction for gamma)

[wilbert]

When you posterize or dither an sRGB image without correcting for gamma, it will be too bright.
The less levels you choose, the more noticeable this is.
When you search the internet for gamma correct or gamma aware dithering, you can find some more detailed explanation.

The code below allows to posterize or dither with or without correcting for gamma.
Most of the code uses integer math and lookup tables to speeds things up.
It should be cross platform compatible.

Code: Select all

Procedure Posterize(Image, Levels, Dither = 0, Grayscale=#False, GammaCorrection = #True)
  
  ; Posterize (2024-08-31)
  
  ; Dither:
  ; 0 = no dither
  ; 1 = Sierra Lite (error diffusion)
  ; 2 = Shiau-Fan (error diffusion)
  ; 3 = 16x16 blue noise matrix
  ; 4 = 16x16 bayer matrix
  ; 5 = 6x6 clustered dot
  ; 6 = 6x8 clustered dot
  ; 7 = 6x6 diagonal lines 1
  ; 8 = 6x6 diagonal lines 2
  
  ; Error diffusion matrices:
  ; Sierra Lite       Shiau-Fan
  ;        X  1/2                        X  1/2
  ; 1/4  1/4          1/16  1/16  1/8  1/4
  
  ; limit levels [2, 64]
  If Levels<2: Levels=2: ElseIf Levels>64: Levels=64: EndIf
  Protected.l c0, c1, c2, dx, dy, i, l, m0, m2, mx, my, n, p, t, x, y
  Protected.Ascii *c, *dm, *e, *l
  If IsImage(Image) And StartDrawing(ImageOutput(Image))
    Levels-1
    ; fill lookup tables
    Protected Dim LUT0a.l(255)
    Protected Dim LUT0b.a(255); level index in LUT1
    Protected Dim LUT1a.l(Levels)
    Protected Dim LUT1b.l(Levels); delta to next level
    Protected Dim LUT1c.l(Levels); halfway to next level
    Protected Dim LUT1d.a(Levels); output value for level
    For i = 0 To 255
      If GammaCorrection
        LUT0a(i) = 1e9*Pow(i/255.0, 2.23)
      Else
        LUT0a(i) = $400000*i
      EndIf
    Next
    LUT1a(Levels) = LUT0a(255)
    LUT1c(Levels) = $7fffffff; there is no halfway to next level
    LUT1d(Levels) = 255
    For i = Levels-1 To 0 Step -1
      l = 255*i/Levels
      LUT1a(i) = LUT0a(l)
      LUT1b(i) = LUT1a(i+1)-LUT1a(i)
      LUT1c(i) = (LUT1a(i+1)+LUT1a(i))>>1
      LUT1d(i) = l
    Next
    For i = 0 To 255
      l = Levels*i/255
      If Dither
        LUT0b(i) = l; level index
      Else
        LUT0b(i) = LUT1d(l-(LUT1c(l)-LUT0a(i))>>31); closest color value
      EndIf
    Next
    ; update the image pixels
    mx = OutputWidth()-1: my=OutputHeight()-1
    n = OutputDepth()>>3-2
    If Dither=1 Or Dither=2
      ; error diffusion buffers
      Protected Dim _c0.l(mx+3)
      Protected Dim _c1.l(mx+3)
      Protected Dim _c2.l(mx+3)
    ElseIf Dither
      Protected Dim dm.q(31)
      Select Dither
        Case 4:
          ; 16x16 Bayer data
          dx=16: dy=16
          dm(00)=$a8288808a0208000: dm(01)=$aa2a8a0aa2228202: dm(02)=$68e848c860e040c0: dm(03)=$6aea4aca62e242c2
          dm(04)=$9818b8389010b030: dm(05)=$9a1aba3a9212b232: dm(06)=$58d878f850d070f0: dm(07)=$5ada7afa52d272f2
          dm(08)=$a4248404ac2c8c0c: dm(09)=$a6268606ae2e8e0e: dm(10)=$64e444c46cec4ccc: dm(11)=$66e646c66eee4ece
          dm(12)=$9414b4349c1cbc3c: dm(13)=$9616b6369e1ebe3e: dm(14)=$54d474f45cdc7cfc: dm(15)=$56d676f65ede7efe
          dm(16)=$ab2b8b0ba3238303: dm(17)=$a9298909a1218101: dm(18)=$6beb4bcb63e343c3: dm(19)=$69e949c961e141c1
          dm(20)=$9b1bbb3b9313b333: dm(21)=$9919b9399111b131: dm(22)=$5bdb7bfb53d373f3: dm(23)=$59d979f951d171f1
          dm(24)=$a7278707af2f8f0f: dm(25)=$a5258505ad2d8d0d: dm(26)=$67e747c76fef4fcf: dm(27)=$65e545c56ded4dcd
          dm(28)=$9717b7379f1fbf3f: dm(29)=$9515b5359d1dbd3d: dm(30)=$57d777f75fdf7fff: dm(31)=$55d575f55ddd7dfd
        Case 5:
          ; 6x6 clustered dot (converted from image magick)
          dx=6: dy=6
          dm(00)=$f8dc23156a87b1bf: dm(01)=$7895eacd32075ca3: dm(02)=$87b1bf23156a404e: dm(03)=$4e78a3f8dc32075c
          dm(04)=$95eacd40
        Case 6:
          ; 6x8 clustered dot
          dx=6: dy=8
          dm(00)=$1045a5af855a507a: dm(01)=$c525053ae4efba1b: dm(02)=$9acf8f653070dafa: dm(03)=$efba5a507aa5af85
          dm(04)=$3adafac51b1045e4: dm(05)=$6530709acf8f2505
        Case 7:
          ; 6x6 diagonal lines 1
          dx=6: dy=6
          dm(00)=$4aca27fba651d17c: dm(01)=$19ed98437520f49f: dm(02)=$3cbc6712e691c36e: dm(03)=$59048a35b5600bdf
          dm(04)=$d8832eae
        Case 8:
          ; 6x6 diagonal lines 2
          dx=6: dy=6
          dm(00)=$ed98277cd1fba651: dm(01)=$0b60b5df43196ec3: dm(02)=$d8832e0459ae8a35: dm(03)=$4a20bce6913c1267
          dm(04)=$75caf49f  
        Default:
          ; 16x16 blue noise taken data from LDR_LLL1_42.png by Christoph Peters (CC0 license)
          dx=16: dy=16
          dm(00)=$f78c1a26f289e7ae: dm(01)=$5264cc9130a469c1: dm(02)=$0862d5a97b5e1c71: dm(03)=$05a6187a550e49ce
          dm(04)=$b272e4560199bfd8: dm(05)=$96f825bcedd99736: dm(06)=$e99f2f41b6fe3b2b: dm(07)=$8047d13e6dac851d
          dm(08)=$510fc283d2134f66: dm(09)=$e2b55b9d002a60fb: dm(10)=$7994236beb90aac8: dm(11)=$0b8b16f576deca40
          dm(12)=$a8e1cb5c337720f3: dm(13)=$7034c34db08e14bb: dm(14)=$04f44a0da5db4558: dm(15)=$a2e8821ee5315768
          dm(16)=$3a28b38afabd0395: dm(17)=$27cd639a43d3a186: dm(18)=$da986e541b8165d7: dm(19)=$b93d06ff6f0ceec5
          dm(20)=$197ceccf3f2eadf1: dm(21)=$177d53beaf7f225f: dm(22)=$b74607a0c4e67348: dm(23)=$8da7dc2c92374ef9
          dm(24)=$8f32d6612493095a: dm(25)=$35f66715ead0a36c: dm(26)=$c611ab8450fcb4d4: dm(27)=$c0219ec7754202e3
          dm(28)=$59f074df126a9b10: dm(29)=$78874c0a5d88b129: dm(30)=$3c9c4bba38c9442d: dm(31)=$ddef39b8fde01f7e
      EndSelect
    EndIf
    If DrawingBufferPixelFormat() & (#PB_PixelFormat_24Bits_BGR|#PB_PixelFormat_32Bits_BGR)
      m0=29: m2=77; BGR
    Else
      m0=77: m2=29; RGB
    EndIf
    If DrawingBufferPixelFormat() & #PB_PixelFormat_ReversedY
      *l = DrawingBuffer()+DrawingBufferPitch()*my
      p = -DrawingBufferPitch()
    Else
      *l = DrawingBuffer()
      p = DrawingBufferPitch()
    EndIf
    For y=0 To my
      *c=*l
      If Grayscale
        If Dither<=0
          ; posterize only (gray)
          For x=0 To mx
            l=*c\a*m0: *c+1: l+*c\a*150: *c+1: l=(l+*c\a*m2+128)>>8: *c-2: l=LUT0b(l)
            *c\a=l: *c+1: *c\a=l: *c+1: *c\a=l: *c+n
          Next
        ElseIf Dither=1
          ; Sierra Lite error diffusion dither (gray)
          c0=0
          For x=0 To mx
            l=*c\a*m0: *c+1: l+*c\a*150: *c+1: l=(l+*c\a*m2+128)>>8: *c-2
            i=LUT0b(l): l=LUT0a(l)+c0>>1+(_c0(x)+_c0(x+1))>>2
            i-(LUT1c(i)-l)>>31: c0=l-LUT1a(i): _c0(x)=c0: l=LUT1d(i)
            *c\a=l: *c+1: *c\a=l: *c+1: *c\a=l: *c+n
          Next
        ElseIf Dither=2
          ; Shiau-Fan error diffusion dither (gray)
          c0=0
          For x=0 To mx
            l=*c\a*m0: *c+1: l+*c\a*150: *c+1: l=(l+*c\a*m2+128)>>8: *c-2
            i=LUT0b(l): l=LUT0a(l)+c0>>1+_c0(x)>>2+_c0(x+1)>>3+(_c0(x+2)+_c0(x+3))>>4
            i-(LUT1c(i)-l)>>31: c0=l-LUT1a(i): _c0(x)=c0: l=LUT1d(i)
            *c\a=l: *c+1: *c\a=l: *c+1: *c\a=l: *c+n
          Next
        Else
          ; ordered dither (gray)
          *dm=@dm()+(y%dy)*dx: *e=*dm+dx-2
          For x=0 To mx
            t=*dm\a<<1+1: *dm=*dm+1-(dx&((*e-*dm)>>31)); threshold
            l=*c\a*m0: *c+1: l+*c\a*150: *c+1: l=(l+*c\a*m2+128)>>8: *c-2
            i=LUT0b(l): i-(t*((LUT1b(i))>>9)+LUT1a(i)-LUT0a(l))>>31: l=LUT1d(i)
            *c\a=l: *c+1: *c\a=l: *c+1: *c\a=l: *c+n
          Next
        EndIf
      Else
        If Dither<=0
          ; posterize only (color)
          For x=0 To mx
            *c\a=LUT0b(*c\a): *c+1: *c\a=LUT0b(*c\a): *c+1: *c\a=LUT0b(*c\a): *c+n
          Next
        ElseIf Dither=1
          ; Sierra Lite error diffusion dither (color)
          c0=0: c1=0: c2=0
          For x=0 To mx
            l=LUT0a(*c\a)+c0>>1+(_c0(x)+_c0(x+1))>>2; c0
            i=LUT0b(*c\a): i-(LUT1c(i)-l)>>31: c0=l-LUT1a(i): _c0(x)=c0
            *c\a=LUT1d(i): *c+1     
            l=LUT0a(*c\a)+c1>>1+(_c1(x)+_c1(x+1))>>2; c1
            i=LUT0b(*c\a): i-(LUT1c(i)-l)>>31: c1=l-LUT1a(i): _c1(x)=c1
            *c\a=LUT1d(i): *c+1
            l=LUT0a(*c\a)+c2>>1+(_c2(x)+_c2(x+1))>>2; c2
            i=LUT0b(*c\a): i-(LUT1c(i)-l)>>31: c2=l-LUT1a(i): _c2(x)=c2
            *c\a=LUT1d(i): *c+n        
          Next
        ElseIf Dither=2
          ; Shiau-Fan error diffusion dither (color)
          c0=0: c1=0: c2=0
          For x=0 To mx
            l=LUT0a(*c\a)+c0>>1+_c0(x)>>2+_c0(x+1)>>3+(_c0(x+2)+_c0(x+3))>>4; c0
            i=LUT0b(*c\a): i-(LUT1c(i)-l)>>31: c0=l-LUT1a(i): _c0(x)=c0
            *c\a=LUT1d(i): *c+1     
            l=LUT0a(*c\a)+c1>>1+_c1(x)>>2+_c1(x+1)>>3+(_c1(x+2)+_c1(x+3))>>4; c1
            i=LUT0b(*c\a): i-(LUT1c(i)-l)>>31: c1=l-LUT1a(i): _c1(x)=c1
            *c\a=LUT1d(i): *c+1
            l=LUT0a(*c\a)+c2>>1+_c2(x)>>2+_c2(x+1)>>3+(_c2(x+2)+_c2(x+3))>>4; c2
            i=LUT0b(*c\a): i-(LUT1c(i)-l)>>31: c2=l-LUT1a(i): _c2(x)=c2
            *c\a=LUT1d(i): *c+n        
          Next
        Else
          ; ordered dither (color)
          *dm=@dm()+(y%dy)*dx: *e=*dm+dx-2
          For x=0 To mx
            t=*dm\a<<1+1: *dm=*dm+1-(dx&((*e-*dm)>>31)); threshold
            i=LUT0b(*c\a): i-(t*((LUT1b(i))>>9)+LUT1a(i)-LUT0a(*c\a))>>31; c0
            *c\a=LUT1d(i): *c+1
            i=LUT0b(*c\a): i-(t*((LUT1b(i))>>9)+LUT1a(i)-LUT0a(*c\a))>>31; c1
            *c\a=LUT1d(i): *c+1
            i=LUT0b(*c\a): i-(t*((LUT1b(i))>>9)+LUT1a(i)-LUT0a(*c\a))>>31; c2
            *c\a=LUT1d(i): *c+n
          Next
        EndIf
      EndIf
      *l+p
    Next
    StopDrawing()
    ProcedureReturn #True
  Else
    ProcedureReturn #False
  EndIf
EndProcedure

Here's also an example how to use it.

Code: Select all

UseJPEGImageDecoder()
UsePNGImageDecoder()

If OpenWindow(0, 0, 0, 990, 580, "Image posterization (and dithering)", #PB_Window_SystemMenu | #PB_Window_MinimizeGadget | #PB_Window_ScreenCentered)
  ImageGadget(0, 10, 10, 480, 480, 0, #PB_Image_Border)
  ImageGadget(1, 500, 10, 480, 480, 0, #PB_Image_Border)
  ButtonGadget(2, 20, 520, 120, 30, "Load image")
  TrackBarGadget(3, 180, 520, 220, 30, 2, 32, #PB_TrackBar_Ticks)
  SetGadgetState(3, 4)
  TextGadget(4, 410, 520, 80, 30, "4 Levels")
  ComboBoxGadget(5, 500, 520, 160, 30)
  AddGadgetItem(5, -1, "Posterize only")
  AddGadgetItem(5, -1, "Sierra Lite")
  AddGadgetItem(5, -1, "Shiau-Fan")
  AddGadgetItem(5, -1, "16x16 blue noise")
  AddGadgetItem(5, -1, "16x16 bayer")
  AddGadgetItem(5, -1, "6x6 clustered dot")
  AddGadgetItem(5, -1, "6x8 clustered dot")
  AddGadgetItem(5, -1, "6x6 diagonal lines 1")
  AddGadgetItem(5, -1, "6x6 diagonal lines 2")
  SetGadgetState(5, 0)
  CheckBoxGadget(6, 710, 520, 90, 30, "Grayscale")
  CheckBoxGadget(7, 810, 520, 160, 30, "Gamma correction")
  SetGadgetState(7, #True)  
  
  Repeat
    Event = WaitWindowEvent()
    If Event = #PB_Event_Gadget
      Select EventGadget()
        Case 2; load image
          File.s = OpenFileRequester("Choose image file", GetCurrentDirectory(), "", 0)
          If File And LoadImage(0, File)
            scale.d = DesktopScaledX(480)/ImageWidth(0)
            scaleY.d = DesktopScaledY(480)/ImageHeight(0)
            If scaleY<scale: scale=scaleY: EndIf
            If scale<1
              ResizeImage(0, ImageWidth(0)*scale, ImageHeight(0)*scale)   
            EndIf
            SetGadgetState(0, ImageID(0))
          EndIf
        Case 3; levels
          SetGadgetText(4, Str(GetGadgetState(3))+" Levels")
      EndSelect
      If IsImage(0)
        ; update result image
        CopyImage(0,1)
        Posterize(1, GetGadgetState(3), GetGadgetState(5), GetGadgetState(6), GetGadgetState(7))
        SetGadgetState(1, ImageID(1))
      EndIf
    EndIf
  Until Event = #PB_Event_CloseWindow
EndIf

[wilbert]

I updated the code form my previous post to include some additional dithering options.

[Caronte3D]

Nice!
Thanks!

[Fred]

That's cool stuff

[Seymour Clufley]

Very useful, Wilbert. Thanks for sharing this.

[loulou2522]

Is that's programm allow to have a better image for OCR ?

[wilbert]

Is that's programm allow to have a better image for OCR ?

I'm not an expert on OCR but I don't think it will help with that.

[wilbert]

Another update ...

Changes:
- Grayscale is a procedure argument now. It made more sense to me to do it this way. It also works now for posterize only and error diffusion.
- Added a Shiau-Fan error diffusion algorithm. Its matrix is quite similar to that of Sierra Lite. Sierra Lite and Shiau-Fan are some of the easiest to implement dithering algorithms (easier as Floyd–Steinberg I think) and both require little working memory.
- Added a second diagonal lines dither.

A little note about the Sierra Lite and Shiau-Fan algorithms.
When looking at the Sierra Lite matrix

Code: Select all

 Sierra Lite
        X  1/2
 1/4  1/4

my initial approach was to write 1/4 to the location below and update the one to the left of that by adding 1/4 so when processing the next line, it can be read as one value.
At some point I realized I could also look at it another way.
The matrix also means that for every pixel (x,y), the previous error to add, comes 1/2 from the previous pixel (x-1,y), 1/4 from the pixel at (x, y-1) and 1/4 from the pixel at (x+1, y-1) so instead of adding the error values together (updating) when writing them, I just write them once (the full error value) and read them back from multiple positions. This worked faster.

[matalog]

Hi wilbert,

This is great, thanks for sharing it.

Could you explain to me please, how do you achieve any more than 2 colours output in the diffusion dithering?

This is my own little b/w Floyd Steinberg dithering method, but I would like to make it output 3 or 4 colours. For now I will only be using shades of grey.

Can you possibly suggest how I could amend this to output more than 2 colours?

Code: Select all

StartDrawing(ImageOutput(0))
For y=0 To #height-1
  For x=0 To #width-1
    pixels(x,y)=Red(Point(x,y))
  Next
Next
StopDrawing()


StartDrawing(ImageOutput(1))
For y=1 To #height-2
  For x=1 To #width-2
    col=pixels(x,y)
    newcol=Round(col/255,#PB_Round_Nearest)*255
    
        
    If col<88 								; EDIT	I was able to add this code and it allows something like 3 colours.
      newcol=0
    EndIf
    If col>87 And col<192
      newcol=128
    EndIf
    If col>191 
      newcol=255
      endif
    
    
    
    pixels(x,y)=newcol
    quanterror=col-newcol
    
    Plot(x,y,RGB(newcol,newcol,newcol))
    pixels(x+1,y)=pixels(x+1,y)+ quanterror * 7/16
    pixels(x-1,y+1)=pixels(x-1,y+1)+ quanterror * 3/16
    pixels(x,y+1)=pixels(x,y+1)+ quanterror * 5/16
    pixels(x+1,y+1)=pixels(x+1,y+1)+ quanterror * 1/16
    
 ;Debug pixels(x,y)

  Next
Next

[threedslider]

Nice effects ! It is a good job on that ! Very well

[wilbert]

Could you explain to me please, how do you achieve any more than 2 colours output in the diffusion dithering?

You already answered it yourself in your edit. It's just creating more output levels.
A very simple way to get 4 levels would be something like

Code: Select all

newcol = (col>>6)*85

A little remark.
Currently you are only using the red channel of the input image. If you only use one channel, it works better to use green instead of red.

Nice effects ! It is a good job on that ! Very well

Thanks

[matalog]

I have a new question, how would I get a dithering method to only choose one of 4 colours, the closest one to the original photo? The 4 colours are RGBY.

I have tried a few naive methods, but they seem to land on red far too often, and they aren't dithered yet.

[wilbert]

I have a new question, how would I get a dithering method to only choose one of 4 colours, the closest one to the original photo? The 4 colours are RGBY.

To get the closest color you need some kind of color distance formula (you can find it on the forum) to calculate the distance between your four colors and the source color and then pick the color with the closest distance.
To dither you keep track of the errors for the red, green and blue channels separately.
I'm not sure what you mean by RGBY but if you mean pure red, green, blue and yellow, you will not be able to get a decent result. You will need black and preferably white as well.

[matalog]

I have a new question, how would I get a dithering method to only choose one of 4 colours, the closest one to the original photo? The 4 colours are RGBY.

To get the closest color you need some kind of color distance formula (you can find it on the forum) to calculate the distance between your four colors and the source color and then pick the color with the closest distance.
To dither you keep track of the errors for the red, green and blue channels separately.
I'm not sure what you mean by RGBY but if you mean pure red, green, blue and yellow, you will not be able to get a decent result. You will need black and preferably white as well.

Yes, I mean the only colours allowed to be output will be Red, Green, Blue and Yellow. These are the colours available to use on an old computer I am makimg images for.

I'm with you on the colour distance, I didn't know about that formula, thanks.

Regarding keeping track of the RGB channels seperately, will that not leave the possibility of each channel writing a value, and therefore the possibility of white appearing, which isn't allowed (yellow would be used instead).

I cannot have colours mix in any pixel, the only pixel to be set with any values from 2 channels is yellow, and I can set both red and green to 255 manually.

Edit: Having now had time to try this out, I see why I need the 3 channels, and it seems to be working mostly, and looks terrible .

[wilbert]

Yes, I mean the only colours allowed to be output will be Red, Green, Blue and Yellow. These are the colours available to use on an old computer I am makimg images for.

I'm very curious what kind of computer that would be, that it has only red, green, blue, yellow and no black
And what is your exact goal?
I'd like to give it a try myself if I better understand what you want to do.