Blue noise dithering in APL

This is translated from an earlier article of mine written in Chinese. This article talked about the void-cluster method used to generate a dither matrix.

Table of Contents

1. Introduction
2. Putting up all together
3. Bonus

1. Introduction

Unlike Floyd or Atkinson dithering, the blue noise dithering is essentially ordered dithering with a special matrix, which can be expressed as:

  |dither←{(⊢>(⍺⌷⍨(⍴⍺∘⊣)|⊢)¨∘⍳∘⍴)⍵}                           ①
  |dither←{
  |    x y←⍴⍵
  |    a b←⍴⍺
5 |    ⍵>(,⍺)[(x y⍴y⍴⍳b)+⍉y x⍴x⍴b×⍳a]
  |}

①: Notice ⎕io←0.

This essentially means the left threshold matrix argument is spanned over the right matrix that represents an image, and values lower than the threshold is discarded.

    (2 2⍴2 1 3 0)⊂⍤{(⊢>(⍺⌷⍨(⍴⍺∘⊣)|⊢)¨∘⍳∘⍴)⍵}⊃⎕←⊂?5 5⍴4
  
┌─────────┐
│0 2 1 3 1│
│2 2 1 2 0│
│2 0 0 0 3│
│3 3 1 2 2│
│3 2 0 2 1│
└─────────┘
┌─────────┐
│0 1 0 1 0│
│0 1 0 1 0│
│0 0 0 0 1│
│0 1 0 1 0│
│1 1 0 1 0│
└─────────┘

Then we need a spherical Gaussian filter, $e^{- \frac{r^{2}}{2 σ^{2}}}$ . That $σ = 1.5$ gives ideal result.

Then it is so in APL:

gauss←{*-4.5÷⍨(⌽,1∘↓⍤1)(⊖⍪1∘↓)⍵ ⍵⍴(+/2*⍨⊢)¨⍸⍵ ⍵⍴1}

The essential void-cluster method is convolution over torus space. Which can then be described so with code:

cluster←{
    s←a+~2|a←⍴⍵
    g←gauss⊃⌈s÷2 ⋄ (¯1 1×a)↓(1 ¯1×a)↓({+/,g×⍵}⌺s)(⍪⍨⍪⊢)(,⍨,⊢)⍵
}

2. Putting up all together

The two functions, mkbp creates bit-pattern, mkda then creates dither array.

   |∇r←mkbp hm;m;l;v
   | m←2×hm
   | r←?m m⍴2
   |loop:
 5 | l←imax cluster r
   | (l⌷r)←0
   | v←imax void r
   | (v⌷r)←0
   | →(l≡v)/0
10 |∇
   | 
   |∇r←mkda hm;bp;pt;m;all;rank;loc
   | m←2×hm
   | r←m m⍴0
15 | bp←mkbp hm
   | pt←bp
   | rank←¯1++/,bp
   | :While rank≥0
   |     loc←imax cluster pt
20 |     ⍞←'.'
   |     (loc⌷pt)←0
   |     (loc⌷r)←rank
   |     rank-←1
   | :EndWhile
25 | pt←bp
   | rank←+/,bp
   | all←m*2
   | :While rank<all
   |     loc←imax void pt
30 |     ⍞←'*'
   |     (loc⌷pt)←1
   |     (loc⌷r)←rank
   |     rank+←1
   | :EndWhile
35 |∇

3. Bonus

A Common Lisp prototype I've written before the APL one:

   |(defun make-bp (hm)
   |  (let* ((m (* 2 hm))
   |        (bp (make-noise m))
   |        last
 5 |        void)
   |    (loop
   |       (setf last (find-location 1 bp))
   |       (setf (apply #'aref bp last) 0)
   |       (setf void (find-location 0 bp))
10 |       (setf (apply #'aref bp void) 1)
   |       (if (equal void last)
   |           (return bp)))))
   | 
   |(defun make-da (bp m hm)
15 |  (declare (optimize (speed 3)))
   |  (let ((da (make-array (list m m) :element-type 'fixnum))
   |        (hf (* hm m))
   |        (all (* m m)))
   |    (let ((pattern (alexandria:copy-array bp))
20 |          (rank (1- (ones bp m))))
   |      (loop while (>= rank 0)
   |            do (let ((loc (cluster pattern)))
   |                 (setf (apply #'aref pattern loc) 0)
   |                 (setf (apply #'aref da loc) rank)
25 |                 (decf rank))))
   |    (let ((pattern (alexandria:copy-array bp))
   |          (rank (ones bp m)))
   |      (loop while (< rank hf)
   |            do (let ((loc (void pattern)))
30 |                 (setf (apply #'aref pattern loc) 1)
   |                 (setf (apply #'aref da loc) rank)
   |                 (incf rank)))
   |      (loop while (< rank all)
   |            do (let ((loc (cluster pattern)))
35 |                 (setf (apply #'aref pattern loc) 1)
   |                 (setf (apply #'aref da loc) rank)
   |                 (incf rank))))
   |    da))