deconvolve.lisp

   1 #.(progn (require :asdf)
   2          (require :vector)
   3          (require :vol)
   4          (require :psf))
   5
   6 (defpackage :deconvolve
   7   (:use :cl :vector :vol))
   8 (in-package :deconvolve)
   9
  10 (defun draw-sphere (vol radius)
  11   (declare ((simple-array (complex double-float) 3) vol)
  12            (double-float radius)
  13            (values (simple-array (complex double-float) 3) &optional))
  14   (destructuring-bind (z y x)
  15      (array-dimensions vol)
  16    (let ((xh (floor x 2))
  17          (yh (floor y 2))
  18          (zh (floor z 2))
  19          (radius2 (* radius radius)))
  20      (do-box (k j i 0 z 0 y 0 x)
  21        (let* ((a (- i xh))
  22               (b (- j yh))
  23               (c (- k zh))
  24               (r2 (+ (* a a) (* b b) (* c c))))
  25          (when (< r2 radius2)
  26            (setf (aref vol k j i)
  27                  (complex 1d0)))))))
  28   vol)
  29
  30 (defun add-noise (vol avg)
  31   (declare ((simple-array (complex double-float) 3) vol)
  32            (double-float avg)
  33            (values (simple-array (complex double-float) 3) &optional))
  34   (destructuring-bind (z y x)
  35       (array-dimensions vol)
  36     (do-box (k j i 0 z 0 y 0 x)
  37       (setf (aref vol k j i)
  38             (+ (aref vol k j i)
  39                (random avg)))))
  40   vol)
  41
  42 #+nil
  43 (time
  44  (let* ((dx 1d0)
  45         (dz dx)
  46         (x 32)
  47         (y 32)
  48         (z 32)
  49         (psf (psf:intensity-psf z y x (* dx x) (* dz z)
  50                                 :integrand-evaluations 100))
  51         (sphere (make-array (list 32 32 32)
  52                             :element-type '(complex double-float))))
  53
  54    (defparameter *sphere* sphere)
  55    (draw-sphere sphere 4.2d0)
  56    (defparameter *psf* psf)
  57
  58
  59    (write-pgm (normalize-img (cross-section-xz psf))
  60               "/home/martin/tmp/deconv000-psf.pgm")
  61    (write-pgm (normalize-img (cross-section-xz sphere))
  62               "/home/martin/tmp/deconv001-sphere.pgm")
  63
  64    (let ((sphere-x-psf (convolve3 sphere psf)))
  65    ;;  (add-noise sphere-x-psf 100d0)
  66      (defparameter *sphere-x-psf* sphere-x-psf)
  67
  68      (write-pgm (normalize-img (cross-section-xz sphere-x-psf))
  69                 "/home/martin/tmp/deconv002-sphere-x-psf.pgm"))))
  70
  71 (defun divide-minus-1 (result mi ei)
  72   (declare ((simple-array (complex double-float) 3) result mi ei)
  73            (values (simple-array (complex double-float) 3) &optional))
  74   (destructuring-bind (z y x)
  75       (array-dimensions mi)
  76     (do-box (k j i 0 z 0 y 0 x)
  77       (setf (aref result k j i)
  78             (complex (- (/ (realpart (aref mi k j i))
  79                            (realpart (aref ei k j i)))
  80                         1d0)))))
  81   result)
  82
  83 (defun multiply-add (pl q fl)
  84   (declare ((simple-array (complex double-float) 3) fl pl)
  85            (double-float q)
  86            (values (simple-array (complex double-float) 3) &optional))
  87   (destructuring-bind (z y x)
  88       (array-dimensions pl)
  89     (do-box (k j i 0 z 0 y 0 x)
  90       (incf (aref pl k j i)
  91             (complex (* q
  92                         (realpart (aref fl k j i))
  93                         (realpart (aref pl k j i)))))))
  94   pl)
  95
  96 #+nil
  97 (time
  98  (let* ((p *sphere* #+nil (make-array (array-dimensions *sphere-x-psf*)
  99                        :element-type '(complex double-float)
 100                        :initial-element (complex 1d0)))
 101         (qi (make-array (array-dimensions *sphere-x-psf*)
 102                         :element-type '(complex double-float)
 103                         :initial-element (complex 1d0)))
 104         (qs '(1 2 1 4 1 8 1 4 1 16 1 2 1 16 1 2 1
 105               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 106               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 107               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 108               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 109               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 110               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 111               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 112               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 113               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 114               16 1 2 1 16 1 2 1 16 1 2 1 16 1 2 1
 115               1 1 1 1 1 1
 116               ))
 117         (count 2))
 118    (loop for i in qs do
 119         (incf count)
 120         (let ((ei (convolve3 p *psf*)))
 121           (divide-minus-1 qi *sphere-x-psf* ei)
 122           (write-pgm (normalize-img (cross-section-xz qi))
 123                      (format nil "/home/martin/tmp/deconv~3,'0d-qi.pgm" count))
 124           (let ((fl (convolve3 qi *psf*)))
 125             (incf count)
 126             (write-pgm (normalize-img (cross-section-xz fl))
 127                        (format nil "/home/martin/tmp/deconv~3,'0d-fl.pgm" count))
 128             (multiply-add p 1d0 fl)
 129             (incf count)
 130             (write-pgm (normalize-img (cross-section-xz p))
 131                    (format nil "/home/martin/tmp/deconv~3,'0d-pl.pgm" count)))))))