vol.lisp

   1 (defpackage :vol
   2   (:use :cl :sb-alien :sb-c-call)
   3   (:export
   4    #:fftshift2
   5    #:ft2
   6    #:ift2
   7    #:convolve2-circ
   8
   9    #:fftshift3
  10    #:ft3
  11    #:ift3
  12    #:read-pgm
  13    #:write-pgm
  14    #:histogram
  15    #:get-linear-array
  16    #:read-stack
  17    #:square
  18    #:linear-regression
  19    #:clamp
  20    #:interp1
  21    #:interpolate2
  22
  23    #:+forward+
  24    #:+backward+
  25    #:+estimate+
  26
  27    #:do-rectangle
  28    #:do-box
  29    #:with-slice
  30
  31    #:save-stack
  32    #:save-stack-ub8
  33    #:.*
  34    #:.+
  35    #:s*
  36    #:.*2
  37    #:.+2
  38    #:s*2
  39    #:convolve3-circ
  40    #:convolve3
  41
  42    #:convert-vol
  43    #:convert-img
  44    #:resample-half
  45    #:cross-section-xz
  46    #:normalize-img
  47    #:normalize-vol))
  48
  49 ;; for i in `cat vol.lisp|grep defconst|cut -d " " -f 2`;do echo \#:$i ;done
  50
  51 (in-package :vol)
  52
  53 #+nil (declaim (optimize (speed 3) (debug 1) (safety 1)))
  54 (declaim (optimize (speed 2) (debug 3) (safety 3)))
  55
  56 (load-shared-object "/usr/lib/libfftw3.so.3")
  57
  58 (define-alien-routine ("fftw_execute" execute)
  59     void
  60   (plan (* int)))
  61
  62 (defconstant +forward+ 1)
  63 (defconstant +backward+ -1)
  64 (defconstant +estimate+ (ash 1 6))
  65
  66
  67 (define-alien-routine ("fftw_plan_dft_3d" plan-dft-3d)
  68     (* int)
  69   (n0 int)
  70   (n1 int)
  71   (n2 int)
  72   (in (* double-float))
  73   (out (* double-float))
  74   (sign int)
  75   (flags unsigned-int))
  76
  77 (define-alien-routine ("fftw_plan_dft_2d" plan-dft-2d)
  78     (* int)
  79   (n0 int)
  80   (n1 int)
  81   (in (* double-float))
  82   (out (* double-float))
  83   (sign int)
  84   (flags unsigned-int))
  85
  86 ;; C-standard "row-major" order, so that the last dimension has the
  87 ;; fastest-varying index in the array.
  88
  89
  90
  91 (defmacro do-rectangle ((j i ymin ymax xmin xmax) &body body)
  92   "Loop through 2d points in ymin .. ymax-1. e.g. call with 0 y 0 x to
  93 visit every point."
  94   `(loop for ,j from ,ymin below ,ymax do
  95         (loop for ,i from ,xmin below ,xmax do
  96              (progn ,@body))))
  97
  98 (defmacro do-box ((k j i zmin zmax ymin ymax xmin xmax) &body body)
  99   "Loop through 3d points e.g. call with 0 z 0 y 0 x to visit every
 100 point."
 101   `(loop for ,k from ,zmin below ,zmax do
 102         (loop for ,j from ,ymin below ,ymax do
 103              (loop for ,i from ,xmin below ,xmax do
 104                   (progn ,@body)))))
 105
 106 (defun fftshift1 (in)
 107   (declare ((simple-array (complex double-float) 1) in)
 108            (values (simple-array (complex double-float) 1) &optional))
 109   (let* ((n (length in))
 110          (nh (floor n 2))
 111          (out (make-array n
 112                           :element-type '(complex double-float))))
 113     (dotimes (i (length in))
 114       (let ((ii (mod (+ i nh) n)))
 115         (setf (aref out ii) (aref in i))))
 116     out))
 117 #+nil
 118 (let* ((ls '(1 2 3 4 5 6 7 8 9))
 119       (a (make-array (length ls)
 120                      :element-type '(complex double-float)
 121                      :initial-contents (mapcar
 122                                         #'(lambda (z) (coerce z
 123                                                          '(complex double-float)))
 124                                         ls))))
 125   (fftshift1 a))
 126
 127
 128
 129 (defun fftshift2 (in)
 130   (declare ((simple-array (complex double-float) 2) in)
 131            (values (simple-array (complex double-float) 2) &optional))
 132   (let ((out (make-array (array-dimensions in)
 133                          :element-type '(complex double-float))))
 134    (destructuring-bind (w1 w0)
 135        (array-dimensions in)
 136      (let ((wh0 (floor w0 2))
 137            (wh1 (floor w1 2)))
 138       (do-rectangle (j i 0 w1 0 w0)
 139         (let* ((ii (mod (+ i wh0) w0))
 140                (jj (mod (+ j wh1) w1)))
 141           (setf (aref out j i)
 142                 (aref in jj ii))))))
 143    out))
 144
 145 (defun ft2 (in &key (forward t))
 146   (declare ((simple-array (complex double-float) 2) in)
 147            (boolean forward)
 148            (values (simple-array (complex double-float) 2) &optional))
 149   (let ((dims (array-dimensions in)))
 150     (destructuring-bind (y x)
 151         dims
 152       (let* ((out (make-array dims :element-type '(complex double-float))))
 153         (sb-sys:with-pinned-objects (in out)
 154           (let ((p (plan-dft-2d y x
 155                                 (sb-sys:vector-sap
 156                                  (sb-ext:array-storage-vector in))
 157                                 (sb-sys:vector-sap
 158                                  (sb-ext:array-storage-vector out))
 159                                 (if forward
 160                                     +forward+
 161                                    +backward+)
 162                                 +estimate+)))
 163             (execute p)))
 164         (when forward ;; normalize if forward
 165           (let ((1/n (/ 1d0 (* x y))))
 166             (do-rectangle (j i 0 y 0 x)
 167               (setf (aref out j i) (* 1/n (aref out j i))))))
 168         out))))
 169
 170 (defmacro ift2 (in)
 171   `(ft2 ,in :forward nil))
 172
 173
 174 ;; was originally fftshift3 /home/martin/usb/y2009/1123/1.lisp
 175 ;; now it is better and work for odd dimensions
 176 (defun fftshift3 (in)
 177   (declare ((simple-array (complex double-float) 3) in)
 178            (values (simple-array (complex double-float) 3) &optional))
 179   (let ((out (make-array (array-dimensions in)
 180                          :element-type '(complex double-float))))
 181     (destructuring-bind (w2 w1 w0)
 182         (array-dimensions in)
 183       (let ((wh0 (floor w0 2))
 184             (wh1 (floor w1 2))
 185             (wh2 (floor w2 2)))
 186         (do-box (k j i 0 w2 0 w1 0 w0)
 187           (let ((ii (mod (+ i wh0) w0))
 188                 (jj (mod (+ j wh1) w1))
 189                 (kk (mod (+ k wh2) w2)))
 190             (setf (aref out k j i)
 191                   (aref in kk jj ii))))))
 192     out))
 193
 194
 195 (declaim (ftype (function ((simple-array (complex double-float) (* * *))
 196                            &key (:forward boolean))
 197                           (values (simple-array (complex double-float) (* * *))
 198                                   &optional))
 199                 ft3))
 200 (defun ft3 (in &key (forward t))
 201   (let ((dims (array-dimensions in)))
 202    (destructuring-bind (z y x)
 203        dims
 204      (let* ((out (make-array dims :element-type '(complex double-float))))
 205        (sb-sys:with-pinned-objects (in out)
 206          (let ((p (plan-dft-3d z y x
 207                                (sb-sys:vector-sap
 208                                 (sb-ext:array-storage-vector in))
 209                                (sb-sys:vector-sap
 210                                 (sb-ext:array-storage-vector out))
 211                                (if forward
 212                                    +forward+
 213                                    +backward+)
 214                                +estimate+)))
 215            (execute p)))
 216        (when forward ;; normalize if forward
 217          (let ((1/n (/ 1d0 (* x y z))))
 218           (do-box (k j i 0 z 0 y 0 x)
 219             (setf (aref out k j i) (* 1/n (aref out k j i))))))
 220        out))))
 221
 222 (defmacro ift3 (in)
 223   `(ft3 ,in :forward nil))
 224
 225 (defun convolve2-circ (vola volb)
 226   (declare ((simple-array (complex double-float) 2) vola volb)
 227            (values (simple-array (complex double-float) 2) &optional))
 228   (let* ((da (array-dimensions vola))
 229          (db (array-dimensions volb))
 230          (compare-ab (map 'list #'(lambda (x y) (eq x y)) da db)))
 231     (when (some #'null compare-ab)
 232       (error "convolve3-circ expects both input arrays to have the same dimensions.")))
 233   (ift2 (.*2 (ft2 vola) (ft2 volb))))
 234
 235 (declaim
 236  (type (function
 237         (string)
 238         (values (simple-array (unsigned-byte 8) (* *))
 239                 &optional))
 240        read-pgm))
 241 (defun read-pgm (filename)
 242   (with-open-file (s filename)
 243     (unless (equal (symbol-name (read s)) "P5")
 244       (error "no PGM file"))
 245     (let* ((w (read s))
 246            (h (read s))
 247            (grays (read s))
 248            (pos (file-position s))
 249            (data (make-array
 250                   (list h w)
 251                   :element-type '(unsigned-byte 8)))
 252            (data-1d (make-array
 253                      (* h w)
 254                      :element-type '(unsigned-byte 8)
 255                      :displaced-to data)))
 256       (declare ((simple-array (unsigned-byte 8) (* *)) data)
 257                ((array (unsigned-byte 8) (*)) data-1d)
 258                ((integer 0 65535) grays w h))
 259       (unless (= grays 255)
 260         (error "image has wrong bitdepth"))
 261       (with-open-file (s filename
 262                          :element-type '(unsigned-byte 8))
 263         (file-position s pos)
 264         (read-sequence data-1d s))
 265       data)))
 266
 267 (declaim
 268  (ftype (function
 269          ((array (unsigned-byte 8) 2)
 270           string)
 271          (values null &optional))
 272         write-pgm))
 273 (defun write-pgm (img filename)
 274   (destructuring-bind (h w)
 275       (array-dimensions img)
 276     (declare ((integer 0 65535) w h))
 277     (with-open-file (s filename
 278                        :direction :output
 279                        :if-exists :supersede
 280                        :if-does-not-exist :create)
 281       (declare (stream s))
 282       (format s "P5~%~D ~D~%255~%" w h))
 283     (with-open-file (s filename
 284                        :element-type '(unsigned-byte 8)
 285                        :direction :output
 286                        :if-exists :append)
 287       (let ((data-1d (make-array
 288                       (* h w)
 289                       :element-type '(unsigned-byte 8)
 290                       :displaced-to img)))
 291         (write-sequence data-1d s)))
 292     nil))
 293
 294 (declaim
 295  (ftype (function
 296          ((array (unsigned-byte 8) (*)) &optional fixnum)
 297          (values (simple-array fixnum (*)) (unsigned-byte 8) (unsigned-byte 8) (unsigned-byte 8)
 298                  &optional))
 299         histogram))
 300 (defun histogram (img &optional (bins 30))
 301   (let* ((maxval (aref img 0))
 302          (minval maxval)
 303          (w (length img)))
 304     (dotimes (i w)
 305       (let ((v (aref img i)))
 306         (when (< v minval)
 307           (setf minval v))
 308         (when (< maxval v)
 309           (setf maxval v))))
 310     (let* ((len (1+ (- maxval minval)))
 311            (result (make-array len :element-type 'fixnum))
 312            (binsm (min (- maxval minval) bins)))
 313       (when (eq maxval minval)
 314         #+nil (error "data is too boring.")
 315         (return-from histogram (values result minval maxval binsm)))
 316       (dotimes (i w)
 317         (incf (aref result (floor (* binsm
 318                                      (- (aref img i)
 319                                         minval))
 320                                   (- maxval minval)))))
 321       (values result minval maxval binsm))))
 322
 323 (declaim (ftype (function ((array * (* *)))
 324                           (values (simple-array * (*)) &optional))
 325                 get-linear-array))
 326 (defun get-linear-array (img)
 327   (sb-ext:array-storage-vector img)
 328   #+nil (make-array (* (array-dimension img 0)
 329                  (array-dimension img 1))
 330               :element-type (array-element-type img)
 331               :displaced-to img))
 332
 333
 334 (declaim (ftype (function (string)
 335                           (values (simple-array (unsigned-byte 8) (* * *))
 336                                   &optional))
 337                 read-stack))
 338 (defun read-stack (fn)
 339   (let* ((files (directory fn))
 340          (slices (length files))
 341          (a (read-pgm (first files))))
 342     (destructuring-bind (h w)
 343         (array-dimensions a)
 344       (let* ((result (make-array (list slices h w)
 345                                  :element-type '(unsigned-byte 8))))
 346         (dotimes (k slices)
 347           (let* ((a (read-pgm (elt files k))))
 348             (dotimes (j h)
 349               (dotimes (i w)
 350                 (setf (aref result k j i)
 351                       (aref a j i))))))
 352         result))))
 353
 354 (defmacro with-slice ((slice-array array slice-nr) &body body)
 355   "Returns SLICE-NRth slice of ARRAY as the 2D SLICE-ARRAY."
 356   (let ((x (gensym))
 357         (y (gensym))
 358         (z (gensym)))
 359     `(destructuring-bind (,z ,y ,x)
 360         (array-dimensions ,array)
 361       (when (or (< ,slice-nr 0) (<= ,z ,slice-nr))
 362         (error "slice-nr=~d out of range [0,~d]" ,slice-nr (1- ,z)))
 363       (let* ((,slice-array (make-array (list ,y ,x)
 364                                        :element-type '(unsigned-byte 8)
 365                                        :displaced-to ,array
 366                                        :displaced-index-offset (* ,slice-nr ,x ,y))))
 367         ,@body))))
 368
 369 (declaim (ftype (function (double-float) (values double-float &optional))
 370                 square))
 371 (defun square (x)
 372   (* x x))
 373
 374 (declaim (ftype (function ((array double-float *)
 375                            &optional
 376                            (array double-float *))
 377                           (values double-float double-float
 378                                   double-float double-float &optional))
 379                 linear-regression))
 380 ;; chernov/book.pdf p. 20
 381 (defun linear-regression (y &optional
 382                           (x (let* ((n (length y)))
 383                                (make-array
 384                                 n
 385                                 :element-type 'double-float
 386                                 :initial-contents
 387                                 (loop for i below n collect (* 1d0 i))))))
 388   "Linear regression of the values in Y with the function y=a*x+b. If
 389 X isn't supplied its assumed to be 0,1,2, ... . Returned are the
 390 fitting parameters A and B and their errors DELTA_A and DELTA_B."
 391   (let* ((n (length y))
 392          (xmean (/ (loop for xi across x sum xi) n))
 393          (ymean (/ (loop for xi across y sum xi) n))
 394          (sxx (loop for xi across x sum (square (- xi xmean))))
 395         #+nil (syy (loop for xi across y sum (square (- xi ymean))))
 396          (sxy (loop for i below n sum (* (- (aref x i) xmean)
 397                                          (- (aref y i) ymean))))
 398          (bhat (/ sxy sxx))
 399          (ahat (- ymean (* bhat xmean)))
 400          (var (/ (loop for i below n sum (square (- (aref y i) ahat
 401                                                        (* bhat (aref x i)))))
 402                     (- n 2)))
 403          (vara (* var (+ (/ (square xmean)
 404                             sxx)
 405                          (/ n))))
 406          (varb (/ var sxx)))
 407     (values ahat bhat (sqrt vara) (sqrt varb))))
 408
 409 #+nil
 410 (linear-regression (let* ((ll (list 1d0 2.01d0 3d0 4d0)))
 411                      (make-array (length ll)
 412                       :element-type 'double-float
 413                       :initial-contents ll)))
 414
 415 (declaim (ftype (function (integer)
 416                           (values (unsigned-byte 8) &optional))
 417                 clamp))
 418 (defun clamp (a)
 419   (when (< a 0)
 420     (return-from clamp 0))
 421   (when (< 255 a)
 422     (return-from clamp 255))
 423   a)
 424
 425
 426 (declaim (ftype (function ((simple-array (unsigned-byte 8) 2)
 427                            fixnum fixnum)
 428                           (values (unsigned-byte 8) &optional))
 429                 aref2-zero-ub8))
 430 (defun aref2-zero-ub8 (a j i)
 431   "Like AREF but return zero if subscripts J and I point outside of
 432 the array bounds."
 433   (if (array-in-bounds-p a j i)
 434       (aref a j i)
 435       0))
 436
 437 (declaim (ftype (function ((simple-array double-float 2)
 438                            fixnum fixnum)
 439                           (values double-float &optional))
 440                 aref2-zero-df))
 441 (defun aref2-zero-df (a j i)
 442   "Like AREF but return zero if subscripts J and I point outside of
 443 the array bounds."
 444   (if (array-in-bounds-p a j i)
 445       (aref a j i)
 446       0d0))
 447
 448 (declaim (ftype (function ((simple-array (complex double-float) 2)
 449                            fixnum fixnum)
 450                           (values (complex double-float) &optional))
 451                 aref2-zero-cdf))
 452 (defun aref2-zero-cdf (a j i)
 453   "Like AREF but return zero if subscripts J and I point outside of
 454 the array bounds."
 455   (if (array-in-bounds-p a j i)
 456       (aref a j i)
 457       (complex 0d0)))
 458
 459
 460 (declaim (ftype (function ((unsigned-byte 8)
 461                            (unsigned-byte 8)
 462                            double-float)
 463                           (values double-float &optional))
 464                 interp1-ub8)
 465          (inline interp1-ub8))
 466 (defun interp1-ub8 (a b xi)
 467     "Interpolate between values A and B. Returns A for xi=0 and B for
 468   xi=1."
 469   (+ (* (- 1d0 xi) a) (* xi b)))
 470
 471 (declaim (ftype (function (double-float double-float double-float)
 472                           (values double-float &optional))
 473                 interp1-df)
 474          (inline interp1-df))
 475 (defun interp1-df (a b xi)
 476     "Interpolate between values A and B. Returns A for xi=0 and B for
 477   xi=1."
 478   (+ (* (- 1d0 xi) a) (* xi b)))
 479
 480 (declaim (ftype (function ((complex double-float)
 481                            (complex double-float) double-float)
 482                           (values (complex double-float) &optional))
 483                 interp1-cdf)
 484          (inline interp1-cdf))
 485 (defun interp1-cdf (a b xi)
 486     "Interpolate between values A and B. Returns A for xi=0 and B for
 487   xi=1."
 488   (+ (* (- 1d0 xi) a) (* xi b)))
 489
 490 (declaim (ftype (function ((or (unsigned-byte 8)
 491                                double-float
 492                                (complex double-float))
 493                            (or (unsigned-byte 8)
 494                                double-float
 495                                (complex double-float))
 496                            double-float)
 497                           (values (or double-float
 498                                       (complex double-float)) &optional))
 499                 interp1)
 500          (inline interp1))
 501 (defun interp1 (a b xi)
 502   (etypecase a
 503     ((unsigned-byte 8) (interp1-ub8 a b xi))
 504     (double-float (interp1-df a b xi))
 505     ((complex double-float) (interp1-cdf a b xi))))
 506
 507 (declaim (ftype (function ((simple-array (unsigned-byte 8) 2)
 508                            double-float double-float)
 509                           (values double-float &optional))
 510                 interpolate2-ub8))
 511 (defun interpolate2-ub8 (img x y)
 512   "Bilinear interpolation on an image IMG. The coordinates X and Y can
 513 be floating point values. If they point outside of IMG 0 is returned."
 514   (multiple-value-bind (i ix)
 515       (floor x)
 516     (multiple-value-bind (j jx)
 517         (floor y)
 518       ;; values on grid points, top left is i,j and stored in a
 519       ;; |
 520       ;;-a-b
 521       ;; |
 522       ;; c d
 523       (let ((a (aref img i j))
 524             (b (aref2-zero-ub8 img (1+ i) j))
 525             (c (aref2-zero-ub8 img i (1+ j)))
 526             (d (aref2-zero-ub8 img (1+ i) (1+ j))))
 527         ;; now interpolate verticals
 528         ;; a  b
 529         ;; |  |
 530         ;; e  f
 531         ;; |  |
 532         ;; c  d
 533         (let ((e (interp1-ub8 a c jx))
 534               (f (interp1-ub8 b d jx)))
 535           ;; and now horizontal
 536           ;; e - g - f
 537           (let ((g (interp1-df e f ix)))
 538             g))))))
 539
 540 (declaim (ftype (function ((simple-array double-float 2)
 541                            double-float double-float)
 542                           (values double-float &optional))
 543                 interpolate2-df))
 544 (defun interpolate2-df (img x y)
 545   "Bilinear interpolation on an image IMG. The coordinates X and Y can
 546 be floating point values. If they point outside of IMG 0 is returned."
 547   (multiple-value-bind (i ix)
 548       (floor x)
 549     (multiple-value-bind (j jx)
 550         (floor y)
 551       (let ((a (aref img i j))
 552             (b (aref2-zero-df img (1+ i) j))
 553             (c (aref2-zero-df img i (1+ j)))
 554             (d (aref2-zero-df img (1+ i) (1+ j))))
 555         (let ((e (interp1-df a c jx))
 556               (f (interp1-df b d jx)))
 557           (let ((g (interp1-df e f ix)))
 558             g))))))
 559
 560 (declaim (ftype (function ((simple-array (complex double-float) 2)
 561                            double-float double-float)
 562                           (values (complex double-float) &optional))
 563                 interpolate2-cdf))
 564 (defun interpolate2-cdf (img x y)
 565   "Bilinear interpolation on an image IMG. The coordinates X and Y can
 566 be floating point values. If they point outside of IMG 0 is returned."
 567   (multiple-value-bind (i ix)
 568       (floor x)
 569     (multiple-value-bind (j jx)
 570         (floor y)
 571       (let ((a (aref img i j))
 572             (b (aref2-zero-cdf img (1+ i) j))
 573             (c (aref2-zero-cdf img i (1+ j)))
 574             (d (aref2-zero-cdf img (1+ i) (1+ j))))
 575         (let ((e (interp1-cdf a c jx))
 576               (f (interp1-cdf b d jx)))
 577           (let ((g (interp1-cdf e f ix)))
 578             g))))))
 579
 580
 581 (declaim (ftype (function ((simple-array * 2) double-float double-float)
 582                           (values (or double-float
 583                                       (complex double-float)) &optional))
 584                 interpolate2))
 585 (defun interpolate2 (img x y)
 586   "Bilinear interpolation on an image IMG. The coordinates X and Y can
 587 be floating point values. If they point outside of IMG 0 is returned."
 588   (etypecase img
 589     ((simple-array (unsigned-byte 8) 2) (interpolate2-ub8 img x y))
 590     ((simple-array double-float 2) (interpolate2-df img x y))
 591     ((simple-array (complex double-float) 2) (interpolate2-cdf img x y))))
 592
 593 #+nil
 594 (let ((a (make-array (list 2 2) :element-type '(unsigned-byte 8)
 595                      :initial-contents '((1 2) (2 3)))))
 596   (interpolate2 a .5d0 .2d0))
 597
 598 (declaim (ftype (function (string (simple-array (complex double-float) 3)
 599                                   &key (:function function))
 600                           (values null &optional))
 601                 save-stack))
 602 (defun save-stack (fn vol &key (function #'realpart))
 603   ;; add a slash / if there isn't one
 604   (ensure-directories-exist (if (eq (1- (length fn))
 605                                             (position #\/ fn :from-end t))
 606                                         fn
 607                                         (format nil "~a/" fn)))
 608   (destructuring-bind (z y x)
 609       (array-dimensions vol)
 610     (let ((b (make-array (list y x)
 611                          :element-type '(unsigned-byte 8))))
 612       (dotimes (k z)
 613         (do-rectangle (j i 0 y 0 x)
 614           (setf (aref b j i)
 615                 (clamp (floor (* 255 (funcall function (aref vol k j i)))))))
 616         (write-pgm b (format nil "~a/~3,'0d.pgm" fn k)))))
 617   nil)
 618
 619 (declaim (ftype (function (string (simple-array (unsigned-byte 8) 3))
 620                           (values null &optional))
 621                 save-stack-ub8))
 622 (defun save-stack-ub8 (fn vol)
 623   (ensure-directories-exist (if (eq (1- (length fn))
 624                                     (position #\/ fn :from-end t))
 625                                 fn
 626                                 (format nil "~a/" fn)))
 627   (destructuring-bind (z y x)
 628       (array-dimensions vol)
 629     (let ((b (make-array (list y x)
 630                          :element-type '(unsigned-byte 8))))
 631       (dotimes (k z)
 632         (do-rectangle (j i 0 y 0 x)
 633           (setf (aref b j i)
 634                 (aref vol k j i)))
 635         (write-pgm b (format nil "~a/~3,'0d.pgm" fn k)))))
 636   nil)
 637
 638 (defun .*2 (vola volb)
 639   (declare ((simple-array (complex double-float) 2) vola volb)
 640            (values (simple-array (complex double-float) 2) &optional))
 641   (let ((result (make-array (array-dimensions vola)
 642                             :element-type (array-element-type vola))))
 643    (destructuring-bind (y x)
 644        (array-dimensions vola)
 645      (do-rectangle (j i 0 y 0 x)
 646        (setf (aref result j i)
 647              (* (aref vola j i)
 648                 (aref volb j i)))))
 649    result))
 650
 651 (defun .+2 (vola volb)
 652   (declare ((simple-array (complex double-float) 2) vola volb)
 653            (values (simple-array (complex double-float) 2) &optional))
 654   (let ((result (make-array (array-dimensions vola)
 655                             :element-type (array-element-type vola))))
 656    (destructuring-bind (y x)
 657        (array-dimensions vola)
 658      (do-rectangle (j i 0 y 0 x)
 659        (setf (aref result j i)
 660              (+ (aref vola j i)
 661                 (aref volb j i)))))
 662    result))
 663
 664 (defun .* (vola volb)
 665   (declare ((simple-array (complex double-float) 3) vola volb)
 666            (values (simple-array (complex double-float) 3) &optional))
 667   (let ((result (make-array (array-dimensions vola)
 668                             :element-type (array-element-type vola))))
 669    (destructuring-bind (z y x)
 670        (array-dimensions vola)
 671      (do-box (k j i 0 z 0 y 0 x)
 672        (setf (aref result k j i)
 673              (* (aref vola k j i)
 674                 (aref volb k j i)))))
 675    result))
 676
 677 (defun .+ (vola volb)
 678   (declare ((simple-array (complex double-float) 3) vola volb)
 679            (values (simple-array (complex double-float) 3) &optional))
 680   (let ((result (make-array (array-dimensions vola)
 681                             :element-type (array-element-type vola))))
 682    (destructuring-bind (z y x)
 683        (array-dimensions vola)
 684      (do-box (k j i 0 z 0 y 0 x)
 685        (setf (aref result k j i)
 686              (+ (aref vola k j i)
 687                 (aref volb k j i)))))
 688    result))
 689
 690 (declaim (ftype (function (double-float (simple-array (complex double-float) 3))
 691                           (values (simple-array (complex double-float) 3) &optional))
 692                 s*))
 693 (defun s* (s vol)
 694   (let* ((a (sb-ext:array-storage-vector vol))
 695          (n (length a)))
 696     (dotimes (i n)
 697       (setf (aref a i) (* s (aref a i)))))
 698   vol)
 699
 700 (defun s*2 (s vol)
 701   (declare (double-float s)
 702            ((simple-array (complex double-float) 2) vol)
 703            (values (simple-array (complex double-float) 2) &optional))
 704   (let* ((a (sb-ext:array-storage-vector vol))
 705          (n (length a)))
 706     (dotimes (i n)
 707       (setf (aref a i) (* s (aref a i)))))
 708   vol)
 709
 710
 711 (declaim (ftype (function ((simple-array (complex double-float) 3)
 712                            (simple-array (complex double-float) 3))
 713                           (values (simple-array (complex double-float) 3) &optional))
 714                 convolve3-circ convolve3))
 715 (defun convolve3-circ (vola volb)
 716   (let* ((da (array-dimensions vola))
 717          (db (array-dimensions volb))
 718          (compare-ab (map 'list #'(lambda (x y) (eq x y)) da db)))
 719     (when (some #'null compare-ab)
 720       (error "convolve3-circ expects both input arrays to have the same dimensions.")))
 721   (ift3 (.* (ft3 vola) (ft3 volb))))
 722
 723 ;; volb is the kernel
 724 (defun convolve3 (vola volb)
 725   (destructuring-bind (za ya xa)
 726       (array-dimensions vola)
 727     (destructuring-bind (zb yb xb)
 728         (array-dimensions volb)
 729       (let* ((biga (make-array (list (+ za zb)
 730                                      (+ ya yb)
 731                                      (+ xa xb))
 732                                :element-type '(complex double-float)))
 733              (bigb (make-array (array-dimensions biga)
 734                                :element-type '(complex double-float)))
 735              (fzb (floor zb 2))
 736              (fyb (floor yb 2))
 737              (fxb (floor xb 2))
 738              (fza (floor za 2))
 739              (fya (floor ya 2))
 740              (fxa (floor xa 2)))
 741         (do-box (k j i 0 za 0 ya 0 xa)
 742           (setf (aref biga (+ k fzb) (+ j fyb) (+ i fxb))
 743                 (aref vola k j i)))
 744         (do-box (k j i 0 zb 0 yb 0 xb)
 745           (setf (aref bigb (+ k fza) (+ j fya) (+ i fxa))
 746                 (aref volb k j i)))
 747         (let* ((conv (convolve3-circ biga (fftshift3 bigb)))
 748                (result (make-array (array-dimensions vola)
 749                                    :element-type '(complex double-float))))
 750           (do-box (k j i 0 za 0 ya 0 xa)
 751             (setf (aref result k j i)
 752                   (aref conv (+ k fzb) (+ j fyb) (+ i fxb))))
 753           result)))))
 754
 755 #+nil
 756 (let ((a (make-array (list 100 200 300)
 757                      :element-type '(complex double-float)))
 758       (b (make-array (list 10 200 30)
 759                      :element-type '(complex double-float))))
 760   (convolve3 a b)
 761   nil)
 762
 763
 764 (declaim (ftype (function ((simple-array (complex double-float) 3))
 765                           (values (simple-array (unsigned-byte 8) 3)
 766                                   &optional))
 767                 convert-vol))
 768 (defun convert-vol (vol)
 769   (destructuring-bind (z y x)
 770       (array-dimensions vol)
 771    (let ((result (make-array (array-dimensions vol)
 772                              :element-type '(unsigned-byte 8))))
 773      (do-box (k j i 0 z 0 y 0 x)
 774        (setf (aref result k j i)
 775              (clamp (floor (* 255 (abs (aref vol k j i)))))))
 776      result)))
 777
 778 (declaim (ftype (function ((simple-array (complex double-float) 2)
 779                            &optional function)
 780                           (values (simple-array (unsigned-byte 8) 2)
 781                                   &optional))
 782                 convert-img))
 783 (defun convert-img (img &optional (function #'abs))
 784   (destructuring-bind (y x)
 785       (array-dimensions img)
 786    (let ((result (make-array (array-dimensions img)
 787                              :element-type '(unsigned-byte 8))))
 788      (do-rectangle (j i 0 y 0 x)
 789        (setf (aref result j i)
 790              (clamp (floor (funcall function (aref img j i))))))
 791      result)))
 792
 793 (declaim (ftype (function ((simple-array (complex double-float) 3))
 794                           (values (simple-array (complex double-float) 3)
 795                                   &optional))
 796                 resample-half))
 797 (defun resample-half (vol)
 798   (destructuring-bind (z y x)
 799       (array-dimensions vol)
 800     (let* ((xx (floor x 2))
 801            (yy (floor y 2))
 802            (zz (floor z 2))
 803            (small (make-array (list zz yy xx)
 804                               :element-type '(complex double-float))))
 805       (do-box (k j i 0 zz 0 xx 0 xx)
 806        (setf (aref small k j i)
 807              (aref vol (* k 2) (* j 2) (* i 2))))
 808       small)))
 809
 810
 811 (declaim (ftype (function ((simple-array (complex double-float) 3)
 812                            &optional fixnum)
 813                           (values (simple-array (complex double-float) 2)
 814                                   &optional))
 815                 cross-section-xz))
 816 (defun cross-section-xz (a &optional (y (floor (array-dimension a 1) 2)))
 817   (destructuring-bind (z y-size x)
 818       (array-dimensions a)
 819     (unless (<= 0 y (1- y-size))
 820       (error "Y is out of bounds."))
 821     (let ((b (make-array (list z x)
 822                          :element-type `(complex double-float))))
 823       (do-rectangle (j i 0 z 0 x)
 824         (setf (aref b j i)
 825               (aref a j y i)))
 826       b)))
 827
 828 (declaim (ftype (function ((simple-array (complex double-float) 2)
 829                            &key (:function function))
 830                           (values (simple-array (unsigned-byte 8) 2)
 831                                   &optional))
 832                 normalize-img))
 833 (defun normalize-img (a &key (function #'abs))
 834   (destructuring-bind (y x)
 835       (array-dimensions a)
 836     (let ((b (make-array (list y x)
 837                          :element-type 'double-float)))
 838       (do-rectangle (j i 0 y 0 x)
 839         (setf (aref b j i)
 840               (funcall function (aref a j i))))
 841       (let* ((b1 (sb-ext:array-storage-vector b))
 842              (ma (reduce #'max b1))
 843              (mi (reduce #'min b1))
 844              (result (make-array (list y x)
 845                                  :element-type '(unsigned-byte 8))))
 846         (when (< (abs (- ma mi)) 1d-12)
 847           (return-from normalize-img result)
 848           #+nil(error "image contains constant data, can't normalize."))
 849         (let ((s (/ 255d0 (- ma mi))))
 850          (do-rectangle (j i 0 y 0 x)
 851            (setf (aref result j i)
 852                  (floor (* s (- (aref b j i) mi))))))
 853         result))))
 854
 855 (declaim (ftype (function ((simple-array (complex double-float) 3)
 856                            &key (:function function))
 857                           (values (simple-array (unsigned-byte 8) 3)
 858                                   &optional))
 859                 normalize-vol))
 860 (defun normalize-vol (a &key (function #'abs))
 861   (destructuring-bind (z y x)
 862       (array-dimensions a)
 863     (let ((b (make-array (list z y x)
 864                          :element-type 'double-float)))
 865       (do-box (k j i 0 z 0 y 0 x)
 866         (setf (aref b k j i)
 867               (funcall function (aref a k j i))))
 868       (let* ((b1 (sb-ext:array-storage-vector b))
 869              (ma (reduce #'max b1))
 870              (mi (reduce #'min b1))
 871              (result (make-array (list z y x)
 872                                  :element-type '(unsigned-byte 8)))
 873              (s (/ 255d0 (- ma mi))))
 874         (do-box (k j i 0 z 0 y 0 x)
 875           (setf (aref result k j i)
 876                 (floor (* s (- (aref b k j i) mi)))))
 877         result))))
 878