vol-misc.lisp

   1 (in-package :vol)
   2
   3 (defun clamp (a)
   4   (declare (integer a)
   5            (values (unsigned-byte 8) &optional))
   6   (when (< a 0)
   7     (return-from clamp 0))
   8   (when (< 255 a)
   9     (return-from clamp 255))
  10   a)
  11
  12 (defun count-non-zero-ub8 (vol)
  13   (declare ((simple-array (unsigned-byte 8) 3) vol)
  14            (values fixnum &optional))
  15   (let* ((sum 0)
  16          (vol1 (sb-ext:array-storage-vector vol)))
  17     (dotimes (i (length vol1))
  18       (when (< 0 (aref vol1 i))
  19         (incf sum)))
  20     sum))
  21
  22 (defun histogram (img &optional (bins 30))
  23   (declare ((array (unsigned-byte 8) 1) img)
  24            (fixnum bins)
  25            (values (simple-array fixnum 1) (unsigned-byte 8) (unsigned-byte 8) (unsigned-byte 8)
  26                  &optional))
  27   (let* ((maxval (aref img 0))
  28          (minval maxval)
  29          (w (length img)))
  30     (dotimes (i w)
  31       (let ((v (aref img i)))
  32         (when (< v minval)
  33           (setf minval v))
  34         (when (< maxval v)
  35           (setf maxval v))))
  36     (let* ((len (1+ (- maxval minval)))
  37            (result (make-array len :element-type 'fixnum))
  38            (binsm (min (- maxval minval) bins)))
  39       (when (eq maxval minval)
  40         #+nil (error "data is too boring.")
  41         (return-from histogram (values result minval maxval binsm)))
  42       (dotimes (i w)
  43         (incf (aref result (floor (* binsm
  44                                      (- (aref img i)
  45                                         minval))
  46                                   (- maxval minval)))))
  47       (values result minval maxval binsm))))
  48
  49 (defun square (x)
  50   (declare (my-float x)
  51            (values my-float &optional))
  52   (* x x))
  53
  54 ;; chernov/book.pdf p. 20
  55 (defun linear-regression (y &optional
  56                           (x (let* ((n (length y)))
  57                                (make-array
  58                                 n
  59                                 :element-type 'my-float
  60                                 :initial-contents
  61                                 (loop for i below n collect (* one i))))))
  62   "Linear regression of the values in Y with the function y=a*x+b. If
  63 X isn't supplied its assumed to be 0,1,2, ... . Returned are the
  64 fitting parameters A and B and their errors DELTA_A and DELTA_B."
  65   (declare ((array my-float 1) y x)
  66            (values my-float my-float my-float my-float &optional))
  67   (let* ((n (length y))
  68          (xmean (/ (loop for xi across x sum xi) n))
  69          (ymean (/ (loop for xi across y sum xi) n))
  70          (sxx (loop for xi across x sum (square (- xi xmean))))
  71         #+nil (syy (loop for xi across y sum (square (- xi ymean))))
  72          (sxy (loop for i below n sum (* (- (aref x i) xmean)
  73                                          (- (aref y i) ymean))))
  74          (bhat (/ sxy sxx))
  75          (ahat (- ymean (* bhat xmean)))
  76          (var (/ (loop for i below n sum (square (- (aref y i) ahat
  77                                                        (* bhat (aref x i)))))
  78                     (- n 2)))
  79          (vara (* var (+ (/ (square xmean)
  80                             sxx)
  81                          (/ n))))
  82          (varb (/ var sxx)))
  83     (values ahat bhat (sqrt vara) (sqrt varb))))
  84
  85 #+nil
  86 (linear-regression (let* ((ll (list one 2.0one 3d0 4d0)))
  87                      (make-array (length ll)
  88                       :element-type 'my-float
  89                       :initial-contents ll)))