2D FFT image quantization

[idle]

An experiment in 2D image FFT quantization between brain farts (as in while I can still think)

Code: Select all

Structure complex
  Re.d
  Im.d
EndStructure

Structure arcomplex
  ar.complex[0]
EndStructure

Procedure _stockham(*x.arcomplex, n.i, flag.i, n2.i, *y.arcomplex)
 
  Protected *y_orig.arcomplex
  Protected *tmp.arcomplex
  
  Protected i.i, j.i, k.i, k2.i, Ls.i, r.i, jrs.i
  Protected half, m, m2
  Protected wr.d, wi.d, tr.d, ti.d
  
  *y_orig = *y
  half = n >> 1
  r = half
  Ls = 1
  
  While(r >= n2)
    *tmp = *x
    *x = *y
    *y = *tmp
    m = 0
    m2 = half
    j = 0
    While j < Ls
      wr = Cos(#PI*j/Ls)
      wi = -flag * Sin(#PI*j/Ls)
      jrs = j*(r+r)
      k = jrs
      While k < jrs+r
        k2 = k + r
        tr =  wr * *y\ar[k2]\Re - wi * *y\ar[k2]\Im
        ti =  wr * *y\ar[k2]\Im + wi * *y\ar[k2]\Re
        *x\ar[m]\Re = *y\ar[k]\Re + tr
        *x\ar[m]\Im = *y\ar[k]\Im + ti
        *x\ar[m2]\Re = *y\ar[k]\Re - tr
        *x\ar[m2]\Im = *y\ar[k]\Im - ti
        m = m + 1
        m2 = m2 + 1
        k = k + 1
      Wend
      j = j + 1
    Wend
    r = r >> 1
    Ls = Ls << 1
  Wend
  
  If (*y <> *y_orig)
    For i = 0 To n -1
      *y\ar[i]\Re = *x\ar[i]\Re
      *y\ar[i]\Im = *x\ar[i]\Im
    Next
  EndIf
  
EndProcedure

Procedure _cooley_tukey(*x.arcomplex, n.i, flag.i, n2.i)
  
  Protected c.complex
  Protected i.i, j.i, k.i, m.i, p.i, n1.i
  Protected Ls.i, ks.i, ms.i, jm.i, dk.i
  Protected wr.d, wi.d, tr.d, ti.d
  Protected tempRe.d, tempIm.d
  
  n1 = n / n2
  k=0
  While k < n1
    j = 0
    m = k
    p = 1
    While p < n1
      j << 1
      j + (m & 1)
      m >> 1
      p << 1
    Wend
    
    If j > k
      i=0
      While i < n2
        tempRe = *x\ar[k*n2+i]\Re
        tempIm = *x\ar[k*n2+i]\Im
        *x\ar[k*n2+i]\Re = *x\ar[j*n2+i]\Re
        *x\ar[k*n2+i]\Im = *x\ar[j*n2+i]\Im
        *x\ar[j*n2+i]\Re = tempRe
        *x\ar[j*n2+i]\Im = tempIm
        i = i + 1
      Wend
    EndIf
    k = k + 1
  Wend
  
  p = 1
  While p < n1
    Ls = p
    p << 1
    jm = 0
    dk = p * n2
    j = 0
    While j < Ls
      wr = Cos((#PI*j/Ls))
      wi = -flag * Sin((#PI*j/Ls))
      k = jm
      While k < n
        ks = k + Ls * n2
        i = 0
        While i < n2
          m = k + i
          ms = ks + i
          tr =  (wr * *x\ar[ms]\Re) - (wi * *x\ar[ms]\Im)
          ti =  (wr * *x\ar[ms]\Im) + (wi * *x\ar[ms]\Re)
          
          tempRe = *x\ar[m]\Re
          tempIm = *x\ar[m]\Im
          
          *x\ar[ms]\Re = tempRe - tr
          *x\ar[ms]\Im = tempIm - ti
          
          *x\ar[m]\Re = tempRe + tr
          *x\ar[m]\Im = tempIm + ti
          
          i = i + 1
        Wend
        k = k + dk
      Wend
      jm = jm + n2
      j = j + 1
    Wend
  Wend
EndProcedure

Procedure fft2D(*x.arcomplex, n1.i, n2.i, flag.i=1)
  Protected *y.arcomplex
  Protected i, n
  
  n = n1*n2
  *y = AllocateMemory(n2*SizeOf(complex))
  i=0
  While i < n
    _stockham(@*x\ar[i], n2, flag, 1, *y)
    i = i + n2
  Wend
  FreeMemory(*y)
  _cooley_tukey(@*x\ar[0], n, flag, n2)
EndProcedure

#block_size = 16
Global img_width = 256
Global img_height = 256
Global image_id.i, output_image_id.i
Global quantization_level.i = 10
Global lquantization_level

Procedure Quantize(*block.arcomplex, size.i, quant_level.d)
 
  Protected i.i
  Protected quant_re.d, quant_im.d
  
  For i = 0 To size - 1
    quant_re = *block\ar[i]\Re / quant_level
    quant_im = *block\ar[i]\Im / quant_level
    *block\ar[i]\Re = Int(quant_re)
    *block\ar[i]\Im = Int(quant_im)
  Next i
EndProcedure

Procedure DeQuantize(*block.arcomplex, size.i, quant_level.d)
  
  Protected i.i
  For i = 0 To size - 1
    *block\ar[i]\Re * quant_level
    *block\ar[i]\Im * quant_level
  Next i
EndProcedure

Procedure UnpackBlockData(*buffer, *Y_block.arcomplex, *Cb_block.arcomplex, *Cr_block.arcomplex, size.i)
  
  Protected i.i, offset.i
  
  If *buffer
   
    For i = 0 To size - 1
      *Y_block\ar[i]\Re = PeekD(*buffer + offset)
      offset + 8
      *Y_block\ar[i]\Im = PeekD(*buffer + offset)
      offset + 8
    Next i
       
    For i = 0 To size - 1
      *Cb_block\ar[i]\Re = PeekD(*buffer + offset)
      offset + 8
      *Cb_block\ar[i]\Im = PeekD(*buffer + offset)
      offset + 8
    Next i
       
    For i = 0 To size - 1
      *Cr_block\ar[i]\Re = PeekD(*buffer + offset)
      offset + 8
      *Cr_block\ar[i]\Im = PeekD(*buffer + offset)
      offset + 8
    Next i
  EndIf
EndProcedure

Procedure PackBlockData(*Y_block.arcomplex, *Cb_block.arcomplex, *Cr_block.arcomplex, size.i)
 
  Protected *buffer, i.i, offset.i
   
  *buffer = AllocateMemory(size * SizeOf(complex) * 3)
  
  If *buffer
   
    For i = 0 To size - 1
      PokeD(*buffer + offset, *Y_block\ar[i]\Re)
      offset + 8 
      PokeD(*buffer + offset, *Y_block\ar[i]\Im)
      offset + 8
    Next i
       
    For i = 0 To size - 1
      PokeD(*buffer + offset, *Cb_block\ar[i]\Re)
      offset + 8
      PokeD(*buffer + offset, *Cb_block\ar[i]\Im)
      offset + 8
    Next i
       
    For i = 0 To size - 1
      PokeD(*buffer + offset, *Cr_block\ar[i]\Re)
      offset + 8
      PokeD(*buffer + offset, *Cr_block\ar[i]\Im)
      offset + 8
    Next i
  EndIf
  
  ProcedureReturn *buffer
EndProcedure

UseZipPacker()

Procedure ProcessImage()
  
  Protected x.i, y.i, i.i, j.i, block_x.i, block_y.i
  Protected R.i, G.i, B.i, Yd.d, Cb.d, Cr.d
  Protected *Y_block.arcomplex, *Cb_block.arcomplex, *Cr_block.arcomplex
  Protected pixel_color.i
  Protected block_count.i, total_size.i
  
  block_count = (img_width / #block_size) * (img_height / #block_size)
    
  original_data_size = block_count * (#block_size * #block_size * SizeOf(complex) * 3)
   
  *Y_block = AllocateMemory(#block_size * #block_size * SizeOf(complex))
  *Cb_block = AllocateMemory(#block_size * #block_size * SizeOf(complex))
  *Cr_block = AllocateMemory(#block_size * #block_size * SizeOf(complex))
    
  Protected *source_buffer = AllocateMemory(original_data_size)
  Protected *compressed_buffer = AllocateMemory(original_data_size)
  Protected offset.i
  
  If ImageID(output_image_id) 
    FreeImage(output_image_id) 
  EndIf 
  
  output_image_id = CreateImage(#PB_Any,ImageWidth(image_id),ImageHeight(image_id)) 
    
  StartDrawing(ImageOutput(image_id)) 
   
  For block_y = 0 To img_height - #block_size Step #block_size
    For block_x = 0 To img_width - #block_size Step #block_size
          
      For j = 0 To #block_size - 1
        For i = 0 To #block_size - 1
          x = block_x + i
          y = block_y + j
          pixel_color = Point(x, y)
          R = Red(pixel_color)
          G = Green(pixel_color)
          B = Blue(pixel_color)
          
          Yd = 0.299 * R + 0.587 * G + 0.114 * B
          Cb = -0.1687 * R - 0.3313 * G + 0.5 * B + 128
          Cr = 0.5 * R - 0.4187 * G - 0.0813 * B + 128
          
          *Y_block\ar[j * #block_size + i]\Re = Yd - 128
          *Cb_block\ar[j * #block_size + i]\Re = Cb - 128
          *Cr_block\ar[j * #block_size + i]\Re = Cr - 128
        Next i
      Next j
           
      fft2D(*Y_block, #block_size, #block_size, 1)
      fft2D(*Cb_block, #block_size, #block_size, 1)
      fft2D(*Cr_block, #block_size, #block_size, 1)
     
      Quantize(*Y_block, #block_size * #block_size, quantization_level)
      Quantize(*Cb_block, #block_size * #block_size, quantization_level * 2)
      Quantize(*Cr_block, #block_size * #block_size, quantization_level * 2)
            
      Protected *block_buffer = PackBlockData(*Y_block, *Cb_block, *Cr_block, #block_size * #block_size)
      CopyMemory(*block_buffer, *source_buffer + offset, #block_size * #block_size * SizeOf(complex) * 3)
      offset + (#block_size * #block_size * SizeOf(complex) * 3)
      FreeMemory(*block_buffer)
      
    Next block_x
  Next block_y
  
  StopDrawing()
  
  
  compressed_data_size = CompressMemory(*source_buffer, original_data_size,*compressed_buffer,original_data_size, #PB_PackerPlugin_Zip)
   
  Debug Str(compressed_data_size) + " " + Str(original_data_size / compressed_data_size)
      
  ;-decode it 
  
  If compressed_data_size 
  
    Protected *decompressed_buffer
    
    If UncompressMemory(*compressed_buffer,compressed_data_size,*source_buffer,original_data_size,#PB_PackerPlugin_Zip) = original_data_size
      *decompressed_buffer = *source_buffer
    EndIf 
    
    If *decompressed_buffer
      StartDrawing(ImageOutput(output_image_id))
      offset = 0
      
      For block_y = 0 To img_height - #block_size Step #block_size
        For block_x = 0 To img_width - #block_size Step #block_size
                    
          UnpackBlockData(*decompressed_buffer + offset, *Y_block, *Cb_block, *Cr_block, #block_size * #block_size)
          offset + (#block_size * #block_size * SizeOf(complex) * 3)
          
          DeQuantize(*Y_block, #block_size * #block_size, quantization_level)
          DeQuantize(*Cb_block, #block_size * #block_size, quantization_level * 2)
          DeQuantize(*Cr_block, #block_size * #block_size, quantization_level * 2)
          
          fft2D(*Y_block, #block_size, #block_size, -1)
          fft2D(*Cb_block, #block_size, #block_size, -1)
          fft2D(*Cr_block, #block_size, #block_size, -1)
                
          Protected norm_factor.d = #block_size * #block_size
          For i = 0 To #block_size * #block_size - 1
            *Y_block\ar[i]\Re / norm_factor
            *Cb_block\ar[i]\Re / norm_factor
            *Cr_block\ar[i]\Re / norm_factor
          Next i
                 
          For j = 0 To #block_size - 1
            For i = 0 To #block_size - 1
              x = block_x + i
              y = block_y + j
              
              Yd = *Y_block\ar[j * #block_size + i]\Re + 128
              Cb = *Cb_block\ar[j * #block_size + i]\Re + 128
              Cr = *Cr_block\ar[j * #block_size + i]\Re + 128
              
              R = Yd + 1.402 * (Cr - 128)
              G = Yd - 0.344136 * (Cb - 128) - 0.714136 * (Cr - 128)
              B = Yd + 1.772 * (Cb - 128)
                          
              If R < 0 : R = 0 : EndIf : If R > 255 : R = 255 : EndIf
              If G < 0 : G = 0 : EndIf : If G > 255 : G = 255 : EndIf
              If B < 0 : B = 0 : EndIf : If B > 255 : B = 255 : EndIf
              
              Plot(x, y, RGB(R, G, B))
            Next i
          Next j
        Next block_x
      Next block_y
      
      DrawText(10,10,Str(compressed_data_size) + " compression ratio " + Str(original_data_size / compressed_data_size))
      
      StopDrawing()
      FreeMemory(*decompressed_buffer)
    EndIf
  EndIf
  
  FreeMemory(*compressed_buffer)
  
  FreeMemory(*Y_block)
  FreeMemory(*Cb_block)
  FreeMemory(*Cr_block)
EndProcedure

Procedure DrawUI(window.i)
 
  Protected event.i, eventGadget.i
  
  If IsWindow(window)
    StartDrawing(WindowOutput(window))
        
    DrawImage(ImageID(image_id), 10, 30)
    DrawingMode(#PB_2DDrawing_Transparent)
    DrawText(10, 10, "Original Image", #Black)
        
    DrawImage(ImageID(output_image_id), img_width + 30, 30)
    DrawText(img_width + 30, 10, "Processed Image", #Black)
    
    StopDrawing()
  EndIf
EndProcedure

UsePNGImageDecoder()

If OpenWindow(0, 0, 0, img_width * 2 + 40, img_height + 100, "FFT Codec", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
  
  
  image_id = LoadImage(#PB_Any, #PB_Compiler_Home + "\Examples\3D\Data\Textures\spheremap.png")
  
  ResizeImage(image_ID,256,256)
  
  If image_id = 0
    MessageRequester("Error", "Check image directory.")
    End
  EndIf
  
  output_image_id = CreateImage(#PB_Any, img_width, img_height)
    
  TrackBarGadget(0, 10, img_height + 40, img_width * 2 + 20, 25, 5, 256)
  SetGadgetState(0, quantization_level)
  lql = quantization_level
  ProcessImage()
  DrawUI(0)
  
  Repeat
    event = WaitWindowEvent()
    Select event
      Case #PB_Event_CloseWindow
        Break
        
      Case #PB_Event_Gadget
        eventGadget = EventGadget()
        If eventGadget = 0
                   
          quantization_level = GetGadgetState(0)
          If lquantization_level <> quantization_level  
            ProcessImage()
            lquantization_level = quantization_level 
          EndIf 
          DrawUI(0)
          
        EndIf
        
    EndSelect
  ForEver
  
EndIf

To do
adaptive quantization to reduce block effects
and a block dictionary to do stream compression

[threedslider]

@idle : Nice work on that

But warnings my antivirus says as threads is a virus so don't worry I deactivated and it works

Keep it up !

[RASHAD]

Very nice work idle
I have a feeling that you are very close to to do AI enlarge image instead of compress it using PB
Thanks for sharing

[idle]

Very nice work idle
I have a feeling that you are very close to to do AI enlarge image instead of compress it using PB
Thanks for sharing

That probably wouldn't be to hard for an AI to do, make it bigger.
I actually used Gemini to produce it and while it didn't have to many issues it still took numerous corrections explanations and some coercion after it decided it didn't want to continue.
I supplied it with the FFT code and walked it through what I wanted it to do and to be honest what would have taken me a day or two to wring out my addled mind and debug took less than hour.
Unfortunately I can't really get it to help me any more as the rest of what I'm doing is novel and AI doesn't do novel.

I will eventually make a native lossy image and movie module with frame callbacks so you can render to what ever target you want and maybe add audio as well which will undoubtedly hurt my head. Though as in everything I do, if there's no interest or pat on the head to say good man, I will get bored and do something else.