Merge branch 'master' of github.com:slyrus/opticl

README.md

@@ -268,3 +268,7 @@ expository example code into
 Some examples of using opticl code can be found here:
 
 [http://www.cyrusharmon.org/static/opticl-examples/opticl-examples.xhtml](http://www.cyrusharmon.org/static/opticl-examples/opticl-examples.xhtml)
+
+# Contributors
+
+Thanks to Ivan Chernetsky for contributing code thresholding grayscale images

cluster.lisp

@@ -111,7 +111,7 @@
          (flet ((recompute-means ()
                   (declare (type 8-bit-rgb-image image)
                            (type (simple-array fixnum (* *)) z)
-                           (type (simple-array fixnum *) counts)
+                           (type (simple-array fixnum (*)) counts)
                            (optimize (speed 3)))
                   ;; clear out the old values
                   (dotimes (q k)

opticl.asd

@@ -1,6 +1,7 @@
 
 (asdf:defsystem :opticl
   :name "opticl"
+  :description "A library for representing and processing images"
   :author "Cyrus Harmon <ch-lisp@bobobeach.com>"
   :licence "BSD"
   :serial t
@@ -26,4 +27,5 @@
    (:file "pnm")
    (:file "gif")
    (:file "io")
-   (:file "cluster")))
+   (:file "cluster")
+   (:file "thresholding")))

opticl.lisp

@@ -61,8 +61,8 @@
                                           (when element-type
                                             `(:element-type ,element-type)))))
               (defun ,ctor-function (height width &key
-                                     (initial-element nil initial-element-p)
-                                     (initial-contents nil initial-contents-p))
+                                                    (initial-element nil initial-element-p)
+                                                    (initial-contents nil initial-contents-p))
                 (apply #'make-array (append (list height width)
                                             (when ,(and channels
                                                         (> channels 1))
@@ -76,14 +76,14 @@
      (frobber ()
        `(progn
           ,@(loop for image-spec in *image-types*
-               collect 
-                 (destructuring-bind (name &key channels element-type)
-                     image-spec
-                   `(frob-image ,name
-                                ,@(if channels
-                                      `(:channels ,channels))
-                                ,@(if element-type
-                                      `(:element-type ,element-type))))))))
+                  collect 
+                  (destructuring-bind (name &key channels element-type)
+                      image-spec
+                    `(frob-image ,name
+                                 ,@(if channels
+                                       `(:channels ,channels))
+                                 ,@(if element-type
+                                       `(:element-type ,element-type))))))))
   (frobber))
 
 ;;; support functions/constants for the pixel setf-expander need to
@@ -117,6 +117,11 @@
   (defconstant +max-image-channels+ 4))
 
 (define-setf-expander pixel (image-var y x &environment env)
+  "Sets the (possibly multiple) image intensity value(s) at position
+y, x to the provided value(s). For example, to set pixel [0,0] in an
+rgb-image to R20, G40, B60, one would do (setf (pixel img 0 0) (values
+20 40 60)). With proper type declarations for images, use of this
+macro should yield non-consing setting of image intensity data. "
   (multiple-value-bind (dummies vals newval setter getter)
       (get-setf-expansion image-var env)
     (declare (ignore newval setter))
@@ -139,12 +144,12 @@
                           `(,@vals ,y ,x)
                           stores
                           `(progn (setf ,@(loop for i from 0
-                                             for store in stores
-                                             collect `(aref ,getter ,temp-y ,temp-x ,i)
-                                             collect store))
+                                                for store in stores
+                                                collect `(aref ,getter ,temp-y ,temp-x ,i)
+                                                collect store))
                                   (values ,@stores))
                           `(values ,@(loop for i from 0 below (length stores)
-                                        collect `(aref ,getter ,temp-y ,temp-x ,i)))))))
+                                           collect `(aref ,getter ,temp-y ,temp-x ,i)))))))
           (let ((syms (map-into (make-list +max-image-channels+) #'gensym)))
             (let ((temp-y (gensym))
                   (temp-x (gensym)))
@@ -173,7 +178,7 @@
                                   (setf (aref ,getter ,temp-y ,temp-x 2) ,(elt syms 2))
                                   (setf (aref ,getter ,temp-y ,temp-x 3) ,(elt syms 3))))
                                 (t (loop for i below d
-                                      collect (setf (aref ,getter ,temp-y ,temp-x i) (elt (list ,@syms) i)))))))
+                                         collect (setf (aref ,getter ,temp-y ,temp-x i) (elt (list ,@syms) i)))))))
                          (2 (setf (aref ,getter ,temp-y ,temp-x) ,(elt syms 0))))
                       `(ecase (array-rank ,getter)
                          (3
@@ -199,17 +204,20 @@
                                 (aref ,getter ,temp-y ,temp-x 3)))
                               (t (values-list
                                   (loop for i below d
-                                     collect (aref ,getter ,temp-y ,temp-x i)))))))
+                                        collect (aref ,getter ,temp-y ,temp-x i)))))))
                          (2 (aref ,getter ,temp-y ,temp-x))))))))))
 
 (defmacro pixel (image-var y x &environment env)
+  "Returns the (possibly multiple) image intensity value(s) at
+position y, x. With proper type declarations for images, use of this
+macro should yield non-consing access to image intensity data. "
   (let ((image-dimensions (%get-image-dimensions image-var env)))
     (if image-dimensions
         (progn
           (ecase (length image-dimensions)
             (2 `(aref ,image-var ,y ,x))
             (3 `(values ,@(loop for i below (third image-dimensions)
-                             collect `(aref ,image-var ,y ,x ,i))))))
+                                collect `(aref ,image-var ,y ,x ,i))))))
         `(ecase (array-rank ,image-var)
            (2 (aref ,image-var ,y ,x))
            (3 (ecase (array-dimension ,image-var 2)
@@ -254,63 +262,60 @@ function does that.")
            (,xmax-var (array-dimension ,img 1))
            (,channels ,(when (or (not image-dimensions)
                                  (> (length image-dimensions) 2))
-                             `(when (= (array-rank ,img) 3)
-                                (array-dimension ,img 2)))))
+                         `(when (= (array-rank ,img) 3)
+                            (array-dimension ,img 2)))))
        (declare (ignorable ,channels)
                 (type fixnum ,ymax-var)
                 (type fixnum ,xmax-var))
        ,@body)))
 
 (defmacro do-pixels ((i-var j-var) image &body body)
   (alexandria:with-gensyms (height width)
-    `(with-image-bounds (,height ,width) ,image
+    `(with-image-bounds (,height ,width)
+         ,image
        (loop for ,i-var fixnum below ,height
-          do (loop for ,j-var fixnum below ,width
-                do ,@body)))))
+             do (loop for ,j-var fixnum below ,width
+                      do ,@body)))))
 
 (defmacro set-pixels ((i-var j-var) image &body body)
   (alexandria:with-gensyms (height width)
-    `(with-image-bounds (,height ,width) ,image
+    `(with-image-bounds (,height ,width)
+         ,image
        (loop for ,i-var fixnum below ,height
-          do (loop for ,j-var fixnum below ,width
-                do (setf (pixel ,image ,i-var ,j-var)
-                         (progn
-                           ,@body)))))))
+             do (loop for ,j-var fixnum below ,width
+                      do (setf (pixel ,image ,i-var ,j-var)
+                               (progn ,@body)))))))
 
 (defmacro do-region-pixels ((i-var j-var y1 x1 y2 x2) image &body body)
   (declare (ignorable image))
   `(loop for ,i-var fixnum from ,y1 below ,y2
-      do (loop for ,j-var fixnum from ,x1 below ,x2
-            do ,@body)))
+         do (loop for ,j-var fixnum from ,x1 below ,x2
+                  do ,@body)))
 
 (defmacro set-region-pixels ((i-var j-var y1 x1 y2 x2) image &body body)
   (declare (ignorable image))
   `(loop for ,i-var fixnum from ,y1 below ,y2
-      do (loop for ,j-var fixnum from ,x1 below ,x2
-            do (setf (pixel ,image ,i-var ,j-var)
-                     (progn
-                       ,@body)))))
+         do (loop for ,j-var fixnum from ,x1 below ,x2
+                  do (setf (pixel ,image ,i-var ,j-var)
+                           (progn
+                             ,@body)))))
 
 (defun clear-image (image)
   (with-image-bounds (height width channels)
       image
     (declare (ignore height width))
     (if channels
         (ecase channels
-          (2 (set-pixels (i j) image
-               (values 0 0)))
-          (3 (set-pixels (i j) image
-               (values 0 0 0)))
-          (4 (set-pixels (i j) image
-               (values 0 0 0 0))))
+          (2 (set-pixels (i j) image (values 0 0)))
+          (3 (set-pixels (i j) image (values 0 0 0)))
+          (4 (set-pixels (i j) image (values 0 0 0 0))))
         (set-pixels (i j) image 0)))
   image)
 
-(defun copy-array (src &key
-                   (element-type (array-element-type src))
-                   (fill-pointer (and (array-has-fill-pointer-p src)
-                                      (fill-pointer src)))
-                   (adjustable (adjustable-array-p src)))
+(defun copy-array (src &key (element-type (array-element-type src))
+                            (fill-pointer (and (array-has-fill-pointer-p src)
+                                               (fill-pointer src)))
+                            (adjustable (adjustable-array-p src)))
   "Returns an undisplaced copy of ARRAY, with same fill-pointer and
 adjustability (if any) as the original, unless overridden by the keyword
 arguments."

package.lisp

@@ -136,6 +136,10 @@
    ;; k-means clustering
    #:k-means-cluster-image-pixels
 
+   ;; thresholding
+   #:threshold-image
+   #:min-error-threshold-image
+
    ;;
    ;; deprectated!!! do not use!!!
    #:convert-image-to-grayscale

thresholding.lisp

@@ -0,0 +1,115 @@
+
+;;;
+;;; Thresholding code contributed buy Ivan Chernetsky in 2011.
+;;; See COPYRIGHT file for details.
+
+(in-package :opticl)
+
+(defun threshold-image (image threshold)
+  "Performs simple thresholding of grayscale image and returns a
+binarized image of type 1-bit-gray-image. Before thresholding
+threshold is coerced to type of image's elements.
+
+An error of type type-error is signaled if image is not of
+gray-image type."
+  (etypecase image
+    (gray-image
+     (with-image-bounds (height width) image
+       (let ((binary-image (make-1-bit-gray-image height width :initial-element 0))
+             (threshold (coerce threshold (array-element-type image))))
+         (declare (type gray-image image)
+                  (type 1-bit-gray-image binary-image))
+         (do-pixels (i j) image
+           (when (>= (pixel image i j) threshold)
+             (setf (pixel binary-image i j) 1)))
+         binary-image)))))
+
+;; Below is an implementation of minimum error thresholding
+
+(defconstant +8-bit-values-count+ 256
+  "Number of different values an element of 8-bit-gray-image
+image can be of.")
+
+(defconstant +8-bit-max-value+ 255
+  "Maximum value an element of 8-bit-gray-image image can be.")
+
+(defun compute-histogram (image)
+  "Computes a normalized histogram of 8-bit-gray-image image,
+i.e. an estimate of the probability density function of gray
+levels, and returns it."
+  (declare (type 8-bit-gray-image image))
+  (let ((histogram (make-array +8-bit-values-count+
+                               :element-type 'double-float
+                               :initial-element 0d0)))
+    (do-pixels (i j) image
+      (incf (aref histogram (pixel image i j))))
+    (with-image-bounds (height width) image
+      (let ((size (* height width)))
+        (loop for i below +8-bit-values-count+
+           do (setf (aref histogram i)
+                    (/ (aref histogram i) size)))))
+    histogram))
+
+(defconstant +j-fn-max-value+ most-positive-double-float
+  "A constant that is used as a default value of J(T) function,
+if its value cannot be calculated for a given T.")
+
+(defun sample-j-fn (histogram)
+  "Returns an array of values of J(T) for T in [0; 255].
+The elements at indexex 0 and 255 are set to +j-fn-max-value+."
+  (declare (type (simple-array double-float (256)) histogram))
+  (let* ((histlen (length histogram))
+         (samples (make-array histlen :element-type 'double-float
+                              :initial-element +j-fn-max-value+))
+         (apriori-1 (aref histogram 0)))
+    (declare (type double-float apriori-1))
+    (loop for idx from 1 upto (- histlen 2) do
+         (progn
+           (incf apriori-1 (aref histogram idx))
+           (if (zerop apriori-1)
+               (setf (aref samples idx) +j-fn-max-value+)
+               (let* ((apriori-2 (- 1d0 apriori-1))
+                      (mean-1 (/ (loop for i upto idx
+                                    summing (* (aref histogram i) i))
+                                 apriori-1))
+                      (mean-2 (/ (loop for i from (1+ idx) to (- histlen 1)
+                                    summing (* (aref histogram i) i))
+                                 apriori-2))
+                      (variance-1 (/ (loop for i upto idx
+                                        summing (* (expt (- i mean-1) 2)
+                                                   (aref histogram i)))
+                                     apriori-1))
+                      (variance-2 (/ (loop for i from (1+ idx) to (- histlen 1)
+                                        summing (* (expt (- i mean-2) 2)
+                                                   (aref histogram i)))
+                                     apriori-2)))
+                 (declare (type double-float apriori-2 mean-1 mean-2 variance-1 variance-2))
+                 (if (or (zerop variance-1) (zerop variance-2))
+                     (setf (aref samples idx) +j-fn-max-value+)
+                     (let ((deviation-1 (sqrt variance-1))
+                           (deviation-2 (sqrt variance-2)))
+                       (declare (type double-float deviation-1 deviation-2))
+                       (setf (aref samples idx)
+                             (+ 1
+                                (* 2 (+ (* apriori-1 (log deviation-1))
+                                        (* apriori-2 (log deviation-2))))
+                                (* -2 (+ (* apriori-1 (log apriori-1))
+                                         (* apriori-2 (log apriori-2))))))))))))
+
+    samples))
+
+(defun min-error-threshold-image (image)
+  "Binarize 8-bit-gray-image image with an automatically guessed
+threshold. Returns multiple values: 1-bit-gray-image image and a guessed
+threshold.
+
+For further details, please refer 'Minumum error thresholding' by
+J. Kittler and J. Illingworth."
+  (etypecase image
+    (8-bit-gray-image
+     (let* ((histogram (compute-histogram image))
+            (j-fn-samples (sample-j-fn histogram))
+            (min-sample (reduce #'min j-fn-samples))
+            (threshold (position min-sample j-fn-samples)))
+       (values (the 1-bit-gray-image (threshold-image image (truncate threshold)))
+               threshold)))))