vol-interpolation.lisp

   1 (in-package :vol)
   2
   3 ;; pbourke linearly interpolating points within a box
   4 (def-generator (interpolate (rank type))
   5   (let ((params (ecase rank
   6                   (1 `(xx))
   7                   (2 `(yy xx))
   8                   (3 `(zz yy xx)))))
   9     `(defun ,name (vol ,@params)
  10        "Linear interpolation of a value."
  11       (declare (inline)
  12                ((simple-array ,long-type ,rank) vol)
  13                (single-float ,@params)
  14                (values ,long-type &optional))
  15       (let ((zero ,(coerce 0 long-type))
  16             (one ,(coerce 1 long-type)))
  17         ,(ecase rank
  18                 (1 `(multiple-value-bind (i x)
  19                         (floor xx)
  20                       (labels ((ref (a)
  21                                  (let ((x (+ i a)))
  22                                    (if (array-in-bounds-p vol x)
  23                                        (aref vol x)
  24                                        zero))))
  25                         (let ((o (ref 0))
  26                               (i (ref 1))
  27                               (a (- one x)))
  28                           (+ (* o a)
  29                              (* i x))))))
  30                 (2 `(multiple-value-bind (i x)
  31                         (floor xx)
  32                       (multiple-value-bind (j y)
  33                           (floor yy)
  34                         (labels ((ref (a b)
  35                                    (let ((x (+ i a))
  36                                          (y (+ j b)))
  37                                      (if (array-in-bounds-p vol y x)
  38                                          (aref vol y x)
  39                                          zero))))
  40                           (let ((oo (ref 0 0))
  41                                 (io (ref 1 0))
  42                                 (oi (ref 0 1))
  43                                 (ii (ref 1 1))
  44                                 (a (- one x))
  45                                 (b (- one y)))
  46                             (+ (* oo a b)
  47                                (* io x b)
  48                                (* oi a y)
  49                                (* ii x y)))))))
  50                 (3 `(multiple-value-bind (i x)
  51                         (floor xx)
  52                       (multiple-value-bind (j y)
  53                           (floor yy)
  54                         (multiple-value-bind (k z)
  55                             (floor zz)
  56                           (labels ((ref (a b c)
  57                                      (let ((x (+ i a))
  58                                            (y (+ j b))
  59                                            (z (+ k c)))
  60                                        (if (array-in-bounds-p vol z y x)
  61                                            (aref vol z y x)
  62                                            zero))))
  63                             (let ((ooo (ref 0 0 0))
  64                                   (ioo (ref 1 0 0))
  65                                   (oio (ref 0 1 0))
  66                                   (ooi (ref 0 0 1))
  67                                   (iio (ref 1 1 0))
  68                                   (ioi (ref 1 0 1))
  69                                   (oii (ref 0 1 1))
  70                                   (iii (ref 1 1 1))
  71                                   (a (- one x))
  72                                   (b (- one y))
  73                                   (c (- one z)))
  74                               (+ (* ooo a b c)
  75                                  (* ioo x b c)
  76                                  (* oio a y c)
  77                                  (* ooi a b z)
  78                                  (* iio x y c)
  79                                  (* ioi x b z)
  80                                  (* oii a y z)
  81                                  (* iii x y z)))))))))))))
  82 #+nil
  83 (def-interpolate-rank-type 1 sf)
  84 #+nil
  85 (let ((a (make-array 3
  86                      :element-type 'single-float
  87                      :initial-contents '(1.0 2.0 3.0))))
  88  (interpolate-1-sf a 1.5))
  89
  90 (defmacro def-interpolate-functions (ranks types)
  91   (let ((functions nil)
  92         (cases nil))
  93     (loop for rank in ranks do
  94          (loop for type in types do
  95               (push `(def-interpolate-rank-type ,rank ,type)
  96                     functions)))
  97     (loop for rank in ranks do
  98                        (loop for type in types do
  99                             (let ((long-type (get-long-type type)))
 100                               (push `((simple-array ,long-type ,rank)
 101                                       (,(format-symbol
 102                                          "interpolate-~a-~a" rank type)
 103                                         a ,@(subseq '(x y z) 0 rank)))
 104                                     cases))))
 105     `(progn ,@functions
 106             (defun interpolate (a x &optional y z)
 107               (declare ((or null single-float) x y z))
 108               (etypecase a
 109                 ,@cases
 110          (t (error "type not supported.")))))))
 111
 112 (def-interpolate-functions (1 2 3) (sf df csf cdf))
 113
 114 #+nil
 115 (let* ((n 3)
 116        (a (make-array (list n n n)
 117                       :element-type 'single-float))
 118        (a1 (sb-ext:array-storage-vector a)))
 119   (dotimes (i (length a1))
 120     (setf (aref a1 i) (* 1s0 i)))
 121   (interpolate-3-sf a 1.2s0 1.3s0 1.2s0)
 122   (interpolate a 1.2s0 1.3s0 1.2s0))
 123
 124
 125 (def-generator (resample-half (rank type))
 126   `(defun ,name (vol)
 127      (declare ((simple-array ,long-type ,rank) vol)
 128               (values (simple-array ,long-type ,rank) &optional))
 129      (let ((dims (array-dimensions vol)))
 130        ,(ecase
 131          rank
 132          (1 `(destructuring-bind (x)
 133                  dims
 134                (let* ((xx (floor x 2))
 135                       (small (make-array (list xx)
 136                                         :element-type (quote ,long-type))))
 137                  (do-region ((i) (xx))
 138                    (setf (aref small i)
 139                          (aref vol (* i 2))))
 140                  small)))
 141          (2 `(destructuring-bind (y x)
 142                  dims
 143                (let* ((xx (floor x 2))
 144                       (yy (floor y 2))
 145                       (small (make-array (list yy xx)
 146                                          :element-type (quote ,long-type))))
 147                 (do-region ((j i) (yy xx))
 148                   (setf (aref small j i)
 149                         (aref vol (* j 2) (* i 2))))
 150                 small)))
 151          (3 `(destructuring-bind (z y x)
 152                  dims
 153                (let* ((xx (floor x 2))
 154                       (yy (floor y 2))
 155                       (zz (floor z 2))
 156                       (small (make-array (list zz yy xx)
 157                                          :element-type (quote ,long-type))))
 158                 (do-region ((k j i) (zz yy xx))
 159                   (setf (aref small k j i)
 160                         (aref vol (* k 2) (* j 2) (* i 2))))
 161                 small)))))))
 162 #+nil
 163 (def-resample-half-rank-type 1 ub8)
 164
 165 (defmacro def-resample-half-functions (ranks types)
 166   (let ((result nil))
 167     (loop for rank in ranks do
 168          (loop for type in types do
 169              (push `(def-resample-half-rank-type ,rank ,type)
 170                    result)))
 171     `(progn ,@result)))
 172
 173
 174 (def-resample-half-functions (1 2 3) (ub8 sf df cdf csf))
 175
 176 #+nil
 177 (resample-half-1-ub8
 178  (let ((l '(1 2 3 4 5 6)))
 179    (make-array (length l) :element-type '(unsigned-byte 8)
 180                :initial-contents l)))
 181
 182
 183
 184 ;; The following function is meant to resample an angular PSF that has
 185 ;; been calculated with too much resolution. The important thing is
 186 ;; that the central sample stays central. Consider a 1D array with n
 187 ;; elements and stepsize dx: The length of the array will be
 188 ;; X=n*dx. There are c=center(n) samples right of the center and b
 189 ;; samples (with b+c=n) left of the center. The new step size dx'
 190 ;; suggests that c'=(floor c*dx dx') samples are needed right of the
 191 ;; center and analogously b' samples left of the center. The size of
 192 ;; the new array is n'=b'+c'. The float coordinate x' for the new
 193 ;; sample is (i+b')*dx' where i runs from -b' below c'. The central
 194 ;; sample is i=0. This following function implements this scheme for 3
 195 ;; dimensions.
 196 (def-generator (resample (rank type))
 197   (let ((par-deltas (ecase rank (1 `(dx)) (2 `(dx dy)) (3 `(dx dy dz))))
 198         (par-deltas^ (ecase rank (1 `(dx^)) (2 `(dx^ dy^)) (3 `(dx^ dy^ dz^)))))
 199  `(defun ,name (vol ,@par-deltas ,@par-deltas^)
 200     "Resample by linear interpolation. Ensure that the
 201   former center stays in the center."
 202     (declare ((simple-array ,long-type ,rank) vol)
 203              (single-float ,@par-deltas ,@par-deltas^)
 204              (values (simple-array ,long-type ,rank) &optional))
 205     (labels ((center (n) ;; give the index of the central sample
 206                (declare (fixnum n)
 207                         (values fixnum &optional))
 208                (let ((c (if (evenp n)
 209                             (ceiling n 2)
 210                             (floor n 2))))
 211                  c))
 212              (dims (n d d^) ;; calculate new array size from old one and deltas
 213                (declare (fixnum n)
 214                         (single-float d d^)
 215                         (values fixnum fixnum fixnum &optional))
 216                (let* ((c (center n))
 217                       (b (- n c))
 218                       (c^ (floor (* c d) d^))
 219                       (b^ (floor (* b d) d^))
 220                       (n^ (+ c^ b^)))
 221                  (values b^ c^ n^))))
 222       (macrolet ((coord (i dir) ;; convert index into float coordinate
 223                    `(* (/ ,(format-symbol "D~a^" dir)
 224                            ,(format-symbol "D~a" dir))
 225                         (+ ,i ,(format-symbol "B~a" dir)))))
 226         ,(ecase rank
 227                 (1
 228                  `(destructuring-bind (x)
 229                       (array-dimensions vol)
 230                     (multiple-value-bind (bx cx nx)
 231                         (dims x dx dx^)
 232                       (let ((result (make-array (list nx)
 233                                                 :element-type (quote ,long-type))))
 234                         (do-region ((i) (cx) ((- bx)))
 235                           (let ((xx (coord i x)))
 236                            (setf (aref result (+ bx i))
 237                                  (,(format-symbol "interpolate-~a-~a" rank type) vol xx))))
 238                         result))))
 239                 (2
 240                 `(destructuring-bind (y x)
 241                      (array-dimensions vol)
 242                    (multiple-value-bind (bx cx nx)
 243                        (dims x dx dx^)
 244                      (multiple-value-bind (by cy ny)
 245                          (dims y dy dy^)
 246                        (let ((result (make-array (list ny nx)
 247                                                  :element-type (quote ,long-type))))
 248                          (do-region ((j i) (cy cx) ((- by) (- bx)))
 249                            (let ((xx (coord i x))
 250                                  (yy (coord j y)))
 251                              (setf (aref result (+ by j) (+ bx i))
 252                                    (,(format-symbol "interpolate-~a-~a" rank type) vol xx yy))))
 253                          result)))))
 254                (3
 255                 `(destructuring-bind (z y x)
 256                      (array-dimensions vol)
 257                    (multiple-value-bind (bx cx nx)
 258                        (dims x dx dx^)
 259                      (multiple-value-bind (by cy ny)
 260                          (dims y dy dy^)
 261                        (multiple-value-bind (bz cz nz)
 262                            (dims z dz dz^)
 263                          (let ((result (make-array (list nz ny nx)
 264                                                    :element-type (quote ,long-type))))
 265                            (do-region ((k j i) (cz cy cx) ((- bz) (- by) (- bx)))
 266                              (let ((xx (coord i x))
 267                                    (yy (coord j y))
 268                                    (zz (coord k z)))
 269                                (setf (aref result (+ bz k) (+ by j) (+ bx i))
 270                                      (,(format-symbol "interpolate-~a-~a" rank type) vol xx yy zz))))
 271                            result))))))))))))
 272 #+nil
 273 (def-resample-rank-type 1 sf)
 274 #+nil
 275 (let ((a (make-array 3 :element-type 'single-float :initial-contents '(1s0 2s0 3s0))))
 276  (resample-1-sf a 1.0 .5))
 277
 278 (defmacro def-resample-functions (ranks types)
 279   (let ((result nil))
 280     (loop for rank in ranks do
 281          (loop for type in types do
 282               (push `(def-resample-rank-type ,rank ,type)
 283                     result)))
 284     `(progn ,@result)))
 285
 286 (def-resample-functions (1 2 3) (sf df csf cdf))
 287
 288 #+nil
 289 (time
 290  (let* ((s0 (convert3-ub8/cdf-complex (draw-sphere-ub8 12d0 41 239 232)))
 291         (dr 160d0)
 292         (dz 300d0)
 293         (ddr 200d0)
 294         (ddz 1000d0)
 295         (s1 (resample3-cdf s0 dr dr dz ddr ddr ddz))
 296         (s2 (resample3-cdf s1 ddr ddr ddz dr dr dz)))
 297    (save-stack-ub8 "/home/martin/tmp/s0/" (normalize3-cdf/ub8-realpart s0))
 298    (save-stack-ub8 "/home/martin/tmp/s1/" (normalize3-cdf/ub8-realpart s1))
 299    (save-stack-ub8 "/home/martin/tmp/s2/" (normalize3-cdf/ub8-realpart s2))
 300    (save-stack-ub8 "/home/martin/tmp/s2-s0" (normalize3-cdf/ub8-realpart (.- s0 s2)))))
 301
 302
 303 (def-generator (cross-section-xz (type))
 304   `(defun ,name (a &optional (y (floor (array-dimension a 1) 2)))
 305      (declare ((simple-array ,long-type 3) a)
 306               (fixnum y)
 307               (values (simple-array ,long-type 2)
 308                       &optional))
 309      (destructuring-bind (z y-size x)
 310          (array-dimensions a)
 311        (unless (<= 0 y (1- y-size))
 312          (error "Y is out of bounds."))
 313        (let ((b (make-array (list z x)
 314                             :element-type (quote ,long-type))))
 315          (do-region ((j i) (z x))
 316            (setf (aref b j i)
 317                  (aref a j y i)))
 318          b))))
 319
 320 #+nil
 321 (def-cross-section-xz-type sf)
 322
 323 #+nil
 324 (let ((a (make-array (list 2 2 2) :element-type 'single-float)))
 325  (cross-section-xz-sf a))
 326
 327 (defmacro def-cross-section-xz-functions (types)
 328   `(progn
 329      ,@(loop for type in types collect
 330             `(def-cross-section-xz-type ,type))
 331      (store-new-function (format-symbol "cross-section-xz"))
 332      (defun cross-section-xz (a &optional (y (floor (array-dimension a 1) 2)))
 333        (declare (fixnum y))
 334        (etypecase a
 335          ,@(loop for type in types collect
 336                 (let ((long-type (get-long-type type)))
 337                        `((simple-array ,long-type 3) (,(format-symbol
 338                                                          "cross-section-xz-~a" type) a y))))
 339          (t (error "type not supported."))))))
 340
 341 (def-cross-section-xz-functions (ub8 sf df csf cdf))
 342
 343 #+nil
 344 (let ((a (make-array (list 2 2 2) :element-type 'single-float)))
 345  (cross-section-xz a))
 346
 347 (def-generator (decimate-xy (type))
 348  `(defun ,name (dx vol)
 349    "Increase transversal sampling distance by odd integer factor DX."
 350    (declare (fixnum dx)
 351             ((simple-array ,long-type 3) vol)
 352             (values (simple-array ,long-type 3) &optional))
 353    (unless (eq (mod dx 2) 1)
 354      (error "Factor DX has to be odd."))
 355    (destructuring-bind (z y x)
 356        (array-dimensions vol)
 357      (let* ((dx2 (* dx dx))
 358             (nx (floor x dx))
 359             (ny (floor y dx))
 360             (result (make-array (list z ny nx)
 361                                 :element-type (quote ,long-type))))
 362        (declare (integer nx ny))
 363        (do-region ((k j i) (nx ny z))
 364          (let ((sum (coerce 0 (quote ,long-type))))
 365            (do-region ((jj ii) (dx dx))
 366              (incf sum (aref vol
 367                              k
 368                              (+ (* dx j) jj)
 369                              (+ (* dx i) ii))))
 370            (setf (aref result k j i) ,(if (eq type 'ub8)
 371                                           `(floor sum dx2)
 372                                           `(/ sum dx2)))))
 373        result))))
 374
 375 (defmacro def-decimate-xz-functions (types)
 376   `(progn
 377      ,@(loop for type in types collect
 378             `(def-decimate-xy-type ,type))))
 379
 380 (def-decimate-xz-functions (ub8 sf df csf cdf))
 381
 382 ;; for dx=5:
 383 ;; 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3
 384 ;; x x x q x
 385 ;;           x o x x x
 386 ;;                     x x x x x
 387 ;;                               x x x x x
 388 ;;                                         x x x x ?
 389 ;; dxh=2, n=24
 390 ;; output size floor(n,dx)=4
 391 ;; position of q: ii=3, i=0: dx*i+ii=5*0+3=3
 392 ;; position of o: ii=1, i=1: dx*i+ii=5+1=6
 393
 394 #+nil
 395 (decimate-xy-ub8 5 (make-array (list 2 10 10)
 396                                :element-type '(unsigned-byte 8)))