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

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

cell-geometry: proper infrastructure scripts and directory layout

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