Blue noise dithering

               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.


          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←{(⊢>(⍺⌷⍨(⍴⍺∘⊣)|⊢)¨∘⍳∘⍴)⍵}            [1]
            |dither←{
            |    x y←⍴⍵
            |    a b←⍴⍺
          5 |    ⍵>(,⍺)[(x y⍴y⍴⍳b)+⍉y x⍴x⍴b×⍳a]
            |}

           [1] Notice ⎕io←0.


               This essentially means the left threshold matrix  argu‐
          ment  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^(−sq(r)/2sq(σ)) that σ=1.5 which 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))