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