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.

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 yy⍴⍳b)+⍉y xxb×⍳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 r 2 2 σ 2 . That σ = 1.5 gives ideal result.

Then it is so in APL:

  |gauss{*-4.5÷(⌽,11)(⊖⍪1) (+/2*)¨ 1}

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

  |cluster{
  |    sa+~2|a
  |    ggauss⊃⌈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.

   |rmkbp hm;m;l;v
   | m2×hm
   | r?m m2
   |loop:
 5 | limax cluster r
   | (lr)0
   | vimax void r
   | (vr)0
   | (lv)/0
10 |
   | 
   |rmkda hm;bp;pt;m;all;rank;loc
   | m2×hm
   | rm m0
15 | bpmkbp hm
   | ptbp
   | rank¯1++/,bp
   | :While rank0
   |     locimax cluster pt
20 |     '.'
   |     (locpt)0
   |     (locr)rank
   |     rank-1
   | :EndWhile
25 | ptbp
   | rank+/,bp
   | allm*2
   | :While rank<all
   |     locimax void pt
30 |     '*'
   |     (locpt)1
   |     (locr)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))