source: project/release/4/cell-geometry/trunk/cell-geometry.scm @ 28431

Last change on this file since 28431 was 28431, checked in by Ivan Raikov, 8 years ago

cell-geometry updates and beginnings of test suite

File size: 20.9 KB
Line 
1;;
2;; Representation of geometric objects based on boundaries and points.
3;;
4;;
5;; Copyright 2012-2013 Ivan Raikov and the Okinawa Institute of
6;; Science and Technology.
7;;
8;; This program is free software: you can redistribute it and/or
9;; modify it under the terms of the GNU General Public License as
10;; published by the Free Software Foundation, either version 3 of the
11;; License, or (at your option) any later version.
12;;
13;; This program is distributed in the hope that it will be useful, but
14;; WITHOUT ANY WARRANTY; without even the implied warranty of
15;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16;; General Public License for more details.
17;;
18;; A full copy of the GPL license can be found at
19;; <http://www.gnu.org/licenses/>.
20;;
21
22;; TODO: genpointset: a set of genpoints with bb
23;;       pointset: a set of points with bb
24
25(module cell-geometry
26
27        ( 
28         verbose 
29
30         bounds? make-bounds bounds-empty bounds-add bounds-translate
31         layer-boundary? Bounds BoundsXZ BoundsYZ
32         genpoint? make-genpoint genpoint-coords genpoint-parent-index genpoint-parent-distance
33         originpoint? make-originpoint originpoint-coords originpoint-index
34         cell-type cell-index cell-origin cell-compartments
35         
36         cells-compartments->kd-tree
37         cells-origins->kd-tree
38
39         list->point-tree.bb
40         linear-curve->genpoints.bb
41
42         pointset/uniform-random
43         genpointset/axisline
44         genpointset/angular-perturbation
45
46         load-pointset-from-file
47         )
48
49        (import scheme chicken)
50        (require-extension typeclass kd-tree datatype matchable random-mtzig )
51        (require-library data-structures files posix irregex mathh extras bvsp-spline parametric-curve)
52        (include "mathh-constants.scm")
53
54        (import srfi-1 srfi-4 srfi-13 srfi-63 foreign
55                (only irregex string->irregex irregex-match)
56                (only files make-pathname)
57                (only posix glob)
58                (only data-structures ->string alist-ref compose identity string-split)
59                (only extras fprintf read-line read-lines pp)
60                (only mathh cosh tanh log10)
61                (only kd-tree <KdTree> <Point> Point3d KdTree3d make-point point?
62                      kd-tree? kd-tree-for-each kd-tree-fold-right kd-tree-map
63                      kd-tree->list kd-tree->list* kd-tree-empty? 
64                      kd-tree-size kd-tree-min-index kd-tree-max-index
65                      kd-tree-is-valid? kd-tree-all-subtrees-are-valid?
66                      )
67                (only parametric-curve 
68                      parametric-curve? linear-curve sample-curve
69                      translate-curve iterate-curve bbox-curve)
70                (prefix bvsp-spline bvsp-spline:)
71               
72                )
73       
74        (define verbose (make-parameter 0))
75       
76        (define (d fstr . args)
77          (let ([port (current-error-port)])
78            (if (positive? (verbose)) 
79                (begin (apply fprintf port fstr args)
80                       (flush-output port) ) )))
81       
82        (define (f64vector-empty? x) (zero? (f64vector-length x)))
83       
84        (import-instance (<KdTree> KdTree3d)
85                         (<Point> Point3d))
86
87
88        ;; convenience procedure to access to results of kd-tree-nearest-neighbor queries
89        (define (kdnn-point x) (cadr x))
90        (define (kdnn-distance x) (caddr x))
91        (define (kdnn-index x) (caar x))
92        (define (kdnn-parent-index x) (car (cadar x)))
93        (define (kdnn-parent-distance x) (cadr (cadar x)))
94
95        (define point->list f64vector->list)
96
97        (define (pointset? x) (every point? x))
98
99        (define-record-type bounds
100          (make-bounds top left bottom right)
101          bounds?
102          (top       bounds-top )
103          (left      bounds-left )
104          (bottom    bounds-bottom )
105          (right     bounds-right )
106          )
107
108
109        (define-record-type originpoint (make-originpoint index coords)
110          originpoint? 
111          (coords originpoint-coords)
112          (index originpoint-index)
113          )
114
115        (define-record-type genpoint (make-genpoint coords parent-index parent-distance)
116          genpoint? 
117          (coords genpoint-coords)
118          (parent-index genpoint-parent-index)
119          (parent-distance genpoint-parent-distance)
120          )
121
122        (define (genpointset? x) (every genpoint? x))
123
124        (define-record-type cell (make-cell ty index origin compartments)
125          cell? 
126          (ty cell-type)
127          (index cell-index)
128          (origin cell-origin)
129          (compartments cell-compartments)
130          )
131
132
133        (define (cell-compartment-ref name cell)
134          (alist-ref name (cell-compartments cell)))
135
136
137        (define (cells-compartments->kd-tree cells compartment-name)
138          (let ((t 
139                 (let recur ((cells cells) (points '()))
140                   (if (null? cells)
141                       (list->kd-tree
142                        points
143                        make-value: (lambda (i v) 
144                                      (list (genpoint-parent-index v)
145                                            (genpoint-parent-distance v)))
146                        make-point: (lambda (v) (genpoint-coords v))
147                        )
148                       (let ((cell (car cells)))
149                         (recur (cdr cells) 
150                                (append (cell-compartment-ref compartment-name cell)
151                                        points)))
152                       ))
153                 ))
154            t))
155
156
157        (define (cells-origins->kd-tree cells)
158          (let ((t 
159                 (let recur ((cells cells) (points '()))
160                   (if (null? cells)
161                       (list->kd-tree
162                        points
163                        make-value: (lambda (i v) 
164                                      (list (genpoint-parent-index v)
165                                            (genpoint-parent-distance v)))
166                        make-point: (lambda (v) (genpoint-coords v))
167                        )
168                       (let ((cell (car cells)))
169                         (recur (cdr cells) 
170                                (cons (make-genpoint (cell-origin cell) (cell-index cell) 0.)
171                                      points)))
172                       ))
173                 ))
174            t))
175
176
177        (define bounds-empty (make-bounds -inf.0 +inf.0 +inf.0 -inf.0))
178
179
180        (define (bounds-translate b dx dy)
181          (make-bounds (+ dy (bounds-top b))
182                       (+ dx (bounds-left b))
183                       (+ dy (bounds-bottom b))
184                       (+ dx (bounds-right b))))
185
186
187        (define (bounds-add b p)
188          (make-bounds (fpmax (coord 1  p) (bounds-top b))
189                       (fpmin (coord 0  p) (bounds-left b))
190                       (fpmin (coord 1  p) (bounds-bottom b))
191                       (fpmax (coord 0  p) (bounds-right b))))
192
193
194        (define-datatype layer-boundary layer-boundary?
195          (Bounds (b bounds?))
196          (BoundsXZ (b bounds?) (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
197          (BoundsYZ (b bounds?) (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
198          )
199
200
201        (define (layer-boundary-bounds b)
202          (cases layer-boundary b
203                 (Bounds (b) b)
204                 (BoundsXZ (b n k x y d d2) b)
205                 (BoundsYZ (b n k x y d d2) b)))
206
207
208        (define (boundary-z-extent-function boundary)
209          (cases layer-boundary boundary
210                 (Bounds (b) 
211                         (lambda (x y) 0.))
212                 (BoundsXZ (b n k x y d d2) 
213                           (lambda (xp yp) 
214                             (let-values (((y0tab y1tab y2tab res)
215                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
216                               (f64vector-ref y0tab 0))))
217                 (BoundsYZ (b n k x y d d2) 
218                           (lambda (xp yp) 
219                             (let-values (((y0tab y1tab y2tab res)
220                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
221                               (f64vector-ref y0tab 0))))
222                 ))
223
224
225        (define (point2d-rejection boundary)
226          (let ((top    (bounds-top boundary))
227                (bottom (bounds-bottom boundary))
228                (left   (bounds-left boundary))
229                (right  (bounds-right boundary)))
230            (lambda (p)
231              (let ((x (coord 0 p)) (y (coord 1 p)))
232                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
233            ))
234
235
236        (define (compute-point3d p zu z-extent-function)
237          (let* ((x (coord 0 p))
238                 (y (coord 1 p))
239                 (z-extent (z-extent-function x y)))
240            (if (zero? zu)
241                (make-point x y zu)
242                (make-point x y (fp* zu z-extent)))
243            ))
244
245
246        (define (list->point-tree.bb points #!key (filter-fn identity))
247         
248          (define (update-bb p x-min y-min z-min x-max y-max z-max)
249            (let ((x (coord 0 p)) (y (coord 1 p)) (z (coord 2 p)))
250              (if (< x (x-min)) (x-min x))
251              (if (< y (y-min)) (y-min y))
252              (if (< z (z-min)) (z-min z))
253              (if (> x (x-max)) (x-max x))
254              (if (> y (y-max)) (y-max y))
255              (if (> z (z-max)) (z-max z))
256              ))
257
258
259          (let* (
260                 (t (list->kd-tree
261                     (filter-fn points)
262                     make-point: (lambda (v) (apply make-point (car v)))
263                     make-value: (lambda (i v) (cdr v))
264                     ))
265                 (n (kd-tree-size t))
266                 (x-min (make-parameter +inf.0))
267                 (y-min (make-parameter +inf.0))
268                 (z-min (make-parameter +inf.0))
269                 (x-max (make-parameter -inf.0))
270                 (y-max (make-parameter -inf.0))
271                 (z-max (make-parameter -inf.0))
272                 (bb (begin (kd-tree-for-each
273                             (lambda (p) (update-bb p x-min y-min z-min
274                                                    x-max y-max z-max))
275                             t)
276                            (list (x-min) (y-min) (z-min) (x-max) (y-max) (z-max))))
277                 )
278
279            (list t bb)
280
281            ))
282
283
284        (define (linear-curve->genpoints.bb parent-index k n origin x-coeffs y-coeffs z-coeffs r-coeffs)
285          (let* ((label  (gensym 'curve))
286                 (c      (translate-curve
287                          (append (point->list origin) (list 0.0))
288                          (linear-curve n k (list x-coeffs y-coeffs z-coeffs r-coeffs) 0. 1. )))
289                 (_      (d "c = ~A~%" c))
290                 (_      (d "bbox-curve c = ~A~%" (bbox-curve c)))
291                 (bb     (match-let (((x1 y1 z1 r1 x2 y2 z2 r2) (bbox-curve c)))
292                                    (list (- x1 r1) (- y1 r1) (- z1 r1)
293                                          (+ x2 r2) (+ y2 r2) (+ z2 r2))))
294                 )
295            (let ((point-data (iterate-curve c n)))
296              (let ((x (coord 0 origin))
297                    (y (coord 1 origin))
298                    (z (coord 2 origin)))
299                (let ((pts 
300                       (map
301                        (lambda (p)
302                          (make-genpoint p parent-index (sqrt (dist2 p origin))))
303                        point-data))
304                       )
305                  (d "linear-curve: pts = ~A~%" pts)
306                  (list pts bb)
307                  ))
308              ))
309          )
310
311
312        (define random-state (random-mtzig:init 19))
313
314
315        (define (random-uniform low high)
316          (let ((rlo (if (< low high) low high))
317                (rhi (if (< low high) high low))) 
318            (let ((delta (+ 1 (- rhi rlo)))
319                  (v (random-mtzig:randu! random-state)))
320              (+ rlo (floor (* delta v)))
321              ))
322          )
323
324
325
326        (define (pointset/uniform-random n x-seed y-seed boundary)
327
328          (let* (
329                 (xybounds  (cases layer-boundary boundary
330                                   (Bounds (b) b)
331                                   (BoundsXZ (b n k x y d d2) b)
332                                   (BoundsYZ (b n k x y d d2) b)))
333                 (x-extent   (- (bounds-right xybounds) (bounds-left xybounds)))
334                 (y-extent   (- (bounds-top xybounds) (bounds-bottom xybounds)))
335                 (z-extent-function (boundary-z-extent-function boundary))
336                 )
337
338            (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
339                  (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
340                  (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
341             
342              (let ((point-rejection1 (point2d-rejection xybounds)))
343               
344                (let recur ((i 0) (pts '()))
345                  (if (< i n)
346                      (let ((x (f64vector-ref x-points i))
347                            (y (f64vector-ref y-points i))
348                            (z (f64vector-ref z-points i)))
349                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
350                          (let ((p3d 
351                                 (and (point-rejection1 p)
352                                      (compute-point3d p z z-extent-function))))
353                            (recur (+ 1 i) (if p3d (cons p3d pts) pts)))))
354                     
355                      pts
356                  ))
357              ))
358          ))
359
360          ;; An line that starts at the given origin and
361          ;; follows the given axis
362          (define (genpointset/axisline origins axis len step pp r)
363           
364            (let* ((k 1)
365                   (n (inexact->exact (floor (/ len step))))
366                   (x-coeffs (list (if (= axis 0) 1. 0.) (coord 0 pp)))
367                   (y-coeffs (list (if (= axis 1) 1. 0.) (coord 1 pp)))
368                   (z-coeffs (list (if (= axis 2) 1. 0.) (coord 2 pp)))
369                   (r-coeffs (list 0. r))
370                   )
371             
372              (d "pointset/axisline: step = ~A~%" step)
373              (d "pointset/axisline: len = ~A~%" len)
374              (d "pointset/axisline: n = ~A~%" n)
375             
376              (fold (lambda (op pts)
377                      (assert (originpoint? op))
378                      (let ((i (originpoint-index op))
379                            (p (originpoint-coords op)))
380                        (d "pointset/axisline: p = ~A~%" p)
381                        (cons
382                         (linear-curve->genpoints.bb i k n p x-coeffs y-coeffs z-coeffs r-coeffs)
383                         pts)))
384                    '()
385                    origins)
386           
387              ))
388
389
390
391          (define (genpointset/angular-perturbation
392                   DendriteLmean AxonLmean nDendrites nAxons 
393                   theta-odd-min theta-odd-max theta-even-min theta-even-max
394                   phi-odd-min phi-odd-max phi-even-min phi-even-max
395                   soma-points )
396           
397            (reverse
398             (car
399              (fold (lambda (p ax)
400                   
401                    (d "AngularPerturbation: p = ~A~%" p)
402
403                    (let ((x (coord 0 p))
404                          (y (coord 1 p))
405                          (z (coord 2 p)))
406                     
407                      (let ((clst (car ax))
408                            (gi   (cadr ax)))
409
410                        (d "AngularPerturbation: gi = ~A~%" gi)
411
412                        (let (
413                              (dendpoints 
414                               (list-tabulate 
415                                nDendrites
416                               
417                                (lambda (i)
418
419                                 
420                                  (let* (
421                                         ;; height of dendrite in z direction
422                                         (theta-degrees (if (odd? i)
423                                                            (random-uniform theta-odd-min theta-odd-max)
424                                                            (random-uniform theta-even-min theta-even-max)
425                                                            ))
426                                         (theta  (* (/ PI 180.) theta-degrees ))
427                                         
428                                         ;; rotation of Dendrite in X-Y plane
429                                         (phi-degrees    (if (odd? i)
430                                                             (random-uniform phi-odd-min phi-odd-max)
431                                                             (random-uniform phi-even-min phi-even-max)))
432                                         (phi (* (/ PI 180.) phi-degrees))
433                                         )
434
435                                    (d "theta = ~A phi = ~A ~%" theta-degrees phi-degrees)
436
437                                    (let (
438                                          (Dend-dX (* DendriteLmean (* (cos theta) (cos phi))))
439                                          (Dend-dY (* DendriteLmean (* (cos theta) (sin phi))))
440                                          (Dend-dZ (* DendriteLmean (sin theta)))
441                                          )
442                                     
443                                      (let ((px (+ Dend-dX x))
444                                            (py (+ Dend-dY y))
445                                            (pz (+ Dend-dZ z)))
446
447                                        (let* ((dpoint (make-point px py pz))
448                                               (soma-distance (sqrt (dist2 p dpoint))))
449
450                                          (d "AngularPerturbation: dpoint = ~A~%" dpoint)
451                                          (d "AngularPerturbation: soma-distance = ~A~%" soma-distance)
452
453                                          (make-genpoint dpoint gi soma-distance)
454                                          ))
455                                      ))
456                                  ))
457                               )
458
459                              (axonpoints
460
461                               (list-tabulate 
462                                nAxons
463                               
464                                (lambda (i)
465                                 
466                                  (let (;; height of Axon in z direction
467                                        (theta  (* (/ PI 180.) (random-uniform -5. -15.))) 
468                                        ;; rotation of Axon in X-Y plane
469                                        (phi    (* (/ PI 180.) (random-uniform 0. 360.))) 
470                                        )
471                                    (let ((Axon-dX (* AxonLmean (* (cos theta) (cos phi))))
472                                          (Axon-dY (* AxonLmean (* (cos theta) (sin phi))))
473                                          (Axon-dZ (* AxonLmean (sin theta))))
474                                     
475                                      (let ((px (+ Axon-dX x))
476                                            (py (+ Axon-dY y))
477                                            (pz (+ Axon-dZ z)))
478                                       
479                                        (let* ((apoint (make-point px py pz))
480                                               (soma-distance (sqrt (dist2 p apoint))))
481
482                                          (d "AngularPerturbation: apoint = ~A~%" apoint)
483                                          (d "AngularPerturbation: soma-distance = ~A~%" soma-distance)
484                                          (make-genpoint apoint gi soma-distance)
485                                          ))
486                                      )))
487                                ))
488                              )
489
490                          (list
491                           (cons
492                            (make-cell 'GoC gi p
493                                       `((Dendrites . ,dendpoints)
494                                         (Axons . ,axonpoints))) clst)
495                           (+ 1 gi))
496                          ))
497                      ))
498                  '(() 0)
499                  soma-points))
500           ))
501       
502
503
504       
505        (define comment-pat (string->irregex "^#.*"))
506
507
508        (define (load-pointset-from-file filename)
509
510          (let ((in (open-input-file filename)))
511           
512            (if (not in) (error 'load-pointset-from-file "file not found" filename))
513
514            (let* ((lines
515                    (filter (lambda (line) (not (irregex-match comment-pat line)))
516                            (read-lines in)))
517
518                   (point-data
519                    (map (lambda (line) (apply make-point (map string->number (string-split line " \t"))))
520                         lines))
521
522                   (point-tree (list->kd-tree point-data))
523                   )
524             
525              (list point-tree)
526             
527              ))
528          )
529
530)
Note: See TracBrowser for help on using the repository browser.