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