大津の閾値判別法

大津の閾値判別法

(ql:quickload :iterate)
(ql:quickload :lispbuilder-sdl)
(ql:quickload :lispbuilder-sdl-gfx)
(use-package :iterate)

(defclass <data> ()
  ((max :initform (error "max must specified") :reader data-max :initarg :max)
   (min :initform (error "min must specified") :reader data-min :initarg :min)
   (data :initform (error "data must specified") :reader data :initarg :data)))

(defun my-random (min max)
  (+ min (random (1+ (- max min)))))

(defun make-random-data (len min max)
  (loop repeat len
       for datum = (my-random min max)
       maximize datum into data-max
       minimize datum into data-min
       collect datum into data
       finally (return (make-instance '<data>
				      :min data-min
				      :max data-max
				      :data data))))

(defmethod make-hist-gram ((data <data>))
  (let ((hist-gram (make-hash-table)))
    (loop for datum in (data data)
       do
	 (incf (gethash datum hist-gram 0)))
    hist-gram))

(defun calc-avarage (data)
  (/ (loop for datum in data
	  sum datum)
     (length data)))

(defmethod split-by ((data <data>) (threshold number))
  (loop for datum in (data data)
       if (<= datum threshold)
       collect datum into class1
       else
       collect datum into class2
     finally (return
	       (cons class1 class2))))

(defmethod calc-threshold ((data <data>))
  (iter:iter (for th from (data-min data) below (data-max data))
	     (finding th maximizing
		      (let* ((cls(split-by data th))
                             (c1 (car cls)) (c2 (cdr cls)))
                        (* (length c1)
                           (length c2)
                           (expt (- (calc-avarage c1)
                                    (calc-avarage c2)) 2))))))

(defmethod threshold-viewer ((data <data>))
  (let* ((hist-gram (make-hist-gram data))
         (threshold (calc-threshold data))
         (width (+ (data-max data) 50))
         (height (+ (loop for datum in (data data)
                          maximizing (gethash datum hist-gram))
                    25)))
    (sdl:with-init ()
        (sdl:window width height)
        (sdl:clear-display sdl:*white*)
        ;;Draw data
        (loop for datum in (data data)
              do
              (sdl:draw-line-* 
                (+ 25 datum) height
                (+ 25 datum) (- height (gethash datum hist-gram))
                :color sdl:*black*)) 
        ;;show threshold
        (sdl:draw-line-* 
          threshold 0
          threshold height
          :color sdl:*red*)
        (sdl:update-display)
        (sdl:with-events ()
          (:quit-event () t)))))

CL-USER(3): (calc-threshold (make-random-data 10000 0 255))

128