source: project/release/4/picnic/trunk/picnic-utils.scm @ 31146

Last change on this file since 31146 was 31146, checked in by Ivan Raikov, 6 years ago

picnic: ignoring axon points in swc loader by default

File size: 47.8 KB
Line 
1;;
2;; Neural Parametric Curve Connectivity spatial and geometric utility procedures.
3;;
4;; Copyright 2012-2014 Ivan Raikov and the Okinawa Institute of Science and Technology
5;;
6;; This program is free software: you can redistribute it and/or
7;; modify it under the terms of the GNU General Public License as
8;; published by the Free Software Foundation, either version 3 of the
9;; License, or (at your option) any later version.
10;;
11;; This program is distributed in the hope that it will be useful, but
12;; WITHOUT ANY WARRANTY; without even the implied warranty of
13;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14;; General Public License for more details.
15;;
16;; A full copy of the GPL license can be found at
17;; <http://www.gnu.org/licenses/>.
18;;
19
20(module picnic-utils
21
22        *
23
24
25        (import scheme chicken)
26
27       
28        (require-extension srfi-69 datatype matchable vector-lib regex
29                           mpi mathh typeclass kd-tree random-mtzig
30                           lalr-driver digraph graph-dfs)
31
32
33        (require-library srfi-1 srfi-4 srfi-13 irregex files posix data-structures
34                         bvsp-spline parametric-curve)
35
36
37        (import 
38                (only srfi-1 
39                      fold fold-right filter-map filter every zip list-tabulate delete-duplicates partition 
40                      first second third take concatenate)
41                (only srfi-4 
42                      s32vector s32vector-length s32vector-ref s32vector-set! make-s32vector
43                      f64vector f64vector? f64vector-ref f64vector-set! f64vector-length f64vector->list list->f64vector make-f64vector)
44                (only srfi-13 string= string< string-null? string-prefix? string-trim-both)
45                (only irregex string->irregex irregex-match)
46                (only files make-pathname)
47                (only posix glob find-files)
48                (only extras read-lines pp fprintf )
49                (only ports with-output-to-port )
50                (only data-structures ->string alist-ref compose identity string-split sort atom?)
51                (only lolevel extend-procedure procedure-data extended-procedure?)
52                (prefix bvsp-spline bvsp-spline:)
53                (prefix parametric-curve pcurve:)
54                )
55
56        (define picnic-verbose (make-parameter 0))
57
58
59        (define (d fstr . args)
60          (let ([port (current-error-port)])
61            (if (positive? (picnic-verbose)) 
62                (begin (apply fprintf port fstr args)
63                       (flush-output port) ) )))
64
65
66        (include "mathh-constants")
67
68        (define (sign x) (if (negative? x) -1.0 1.0))
69
70        (define (f64vector-empty? x) (zero? (f64vector-length x)))
71
72        (define (random-init seed) (random-mtzig:init seed))
73
74        (define (random-uniform low high st)
75          (let ((rlo (if (< low high) low high))
76                (rhi (if (< low high) high low))) 
77            (let ((delta (+ 1 (- rhi rlo)))
78                  (v (random-mtzig:randu! st)))
79              (+ rlo (floor (* delta v)))
80              ))
81          )
82
83
84        (define (random-normal mean sdev st)
85          (+ mean (* sdev (random-mtzig:randn! st))))
86
87        (import-instance (<KdTree> KdTree3d)
88                         (<Point> Point3d))
89       
90        ;; convenience procedure to access to results of kd-tree-nearest-neighbor queries
91        (define (kdnn-point x) (cadr x))
92        (define (kdnn-distance x) (caddr x))
93        (define (kdnn-index x) (caar x))
94        (define (kdnn-parent-index x) (car (cadar x)))
95        (define (kdnn-parent-distance x) (cadr (cadar x)))
96
97
98        (define point->list f64vector->list)
99       
100
101       
102        (define-record-type pointset (make-pointset prefix id points boundary)
103          pointset? 
104          (prefix pointset-prefix)
105          (id pointset-id)
106          (points pointset-points)
107          (boundary pointset-boundary)
108          )
109
110        (define-record-type genpoint (make-genpoint coords parent-index parent-distance segment)
111          genpoint? 
112          (coords genpoint-coords)
113          (parent-index genpoint-parent-index)
114          (parent-distance genpoint-parent-distance)
115          (segment genpoint-segment)
116          )
117
118        (define-record-type swcpoint (make-swcpoint id type coords radius pre)
119          swcpoint? 
120          (id swcpoint-id)
121          (type swcpoint-type)
122          (coords swcpoint-coords)
123          (radius swcpoint-radius)
124          (pre swcpoint-pre)
125          )
126
127       
128        (define-record-type cell (make-cell ty index origin sections)
129          cell? 
130          (ty cell-type)
131          (index cell-index)
132          (origin cell-origin)
133          (sections cell-sections)
134          )
135       
136       
137        (define (cell-section-ref name cell)
138          (alist-ref name (cell-sections cell)))
139
140       
141        (define (write-pointset name pts)
142            (call-with-output-file (sprintf "~A.pointset.dat" name)
143              (lambda (out)
144                (for-each (match-lambda
145                           ((gid p)
146                            (fprintf out "~A ~A ~A~%" 
147                                     (coord 0 p)
148                                     (coord 1 p)
149                                     (coord 2 p))))
150                          pts))
151              ))
152
153       
154        (define (write-layout name pts #!optional rank)
155            (call-with-output-file (if rank 
156                                       (sprintf "~A.~A.layout.dat" name rank)
157                                       (sprintf "~A.layout.dat" name))
158              (lambda (out)
159                (for-each (match-lambda
160                           ((gid p)
161                            (fprintf out "~A ~A ~A ~A~%" 
162                                     gid
163                                     (coord 0 p)
164                                     (coord 1 p)
165                                     (coord 2 p))))
166                          pts))
167              ))
168
169        (define (write-sections forest-name section-name layout sections #!optional rank)
170          (call-with-output-file (if rank
171                                     (sprintf "~A.~A.~A.section.dat" forest-name section-name rank)
172                                     (sprintf "~A.~A.section.dat" forest-name section-name))
173            (lambda (out)
174              (for-each
175               (match-lambda*
176                (((gid p) section)
177                 (fprintf out "~A " gid)
178                 (for-each
179                  (lambda (neurites)
180                    (for-each
181                     (lambda (gd)
182                       (let ((p (genpoint-coords gd)))
183                         (fprintf out "~A ~A ~A "
184                                  (coord 0 p)
185                                  (coord 1 p)
186                                  (coord 2 p))))
187                     (cdr neurites)))
188                  (cdr section))
189                 (fprintf out "~%")
190                 ))
191               layout 
192               sections))))
193
194
195        (define (cells-sections->kd-tree cells section-name #!key
196                                         (make-value
197                                          (lambda (i v) 
198                                            (list (genpoint-parent-index v)
199                                                  (genpoint-parent-distance v))))
200                                         (make-point
201                                          (lambda (v) (genpoint-coords v))))
202          (let ((t 
203                 (let recur ((cells cells) (points '()))
204                   (if (null? cells)
205                       (list->kd-tree
206                        points
207                        make-value: make-value
208                        make-point: make-point)
209                       (let ((cell (car cells)))
210                         (recur (cdr cells) 
211                                (let inner ((sections (append (cell-section-ref section-name cell)))
212                                            (points points))
213                                  (if (null? sections) points
214                                      (inner
215                                       (cdr sections)
216                                       (append (cdr (car sections)) points))
217                                      ))
218                                ))
219                       ))
220                 ))
221            t))
222
223        (define (sections->kd-tree cells #!key
224                                   (make-value
225                                    (lambda (i v) 
226                                      (list (genpoint-parent-index v)
227                                            (genpoint-parent-distance v))))
228                                   (make-point
229                                    (lambda (v) (genpoint-coords v))))
230          (let ((t 
231                 (let recur ((cells cells) (points '()))
232                   (if (null? cells)
233                       (list->kd-tree
234                        points
235                        make-value: make-value
236                        make-point: make-point)
237                       (let ((sections (car cells)))
238                         (let inner ((sections sections) (points points))
239                           (if (null? sections)
240                               (recur (cdr cells) points)
241                               (let ((section (car sections)))
242                                 (inner (cdr sections) 
243                                        (append (cdr (cadr section)) points))
244                                 ))
245                           ))
246                       ))
247                 ))
248            t))
249
250
251        (define (cells-origins->kd-tree cells)
252          (let ((t 
253                 (let recur ((cells cells) (points '()))
254                   (if (null? cells)
255                       (list->kd-tree
256                        points
257                        make-value: (lambda (i v) 
258                                      (list (genpoint-parent-index v)
259                                            (genpoint-parent-distance v)))
260                        make-point: (lambda (v) (genpoint-coords v))
261                        )
262                       (let ((cell (car cells)))
263                         (recur (cdr cells) 
264                                (cons (make-genpoint (cell-origin cell) (cell-index cell) 0. 0)
265                                      points)))
266                       ))
267                 ))
268            t))
269
270        (define (pointCoord axis p)
271          (coord (inexact->exact axis) p))
272
273        ;; Samples a parametric curve at regular intervals in the range xmin..xmax inclusive.
274        (define (sample-uniform)
275          (lambda (n)
276            (let* (
277                   (xmin  0.0)
278                   (xmax  1.0)
279                   (delta (- xmax xmin))
280                   (dx    (if (zero? delta) 0
281                              (if (< n 2)
282                                (error 'sample-uniform "number of iterations must be >= 2")
283                                (/ (- xmax xmin) (- n 1)))))
284                   )
285              (list-tabulate n (lambda (i) (+ xmin (* dx i))))
286              )))
287
288
289        ;; Samples a parametric curve according to a polynomial
290        ;; c0 + c1 x + c2 x^2 in the range xmin..xmax inclusive.
291        (define (sample-polynomial c0 c1 c2)
292          (lambda (n)
293            (let* (
294                   (f (lambda (x) (+ c0 (* x c1) (* x x c2))))
295                   (dx (/ 1.0 n))
296                   )
297              (let recur ((i 0.0) (fi 0.0) (lst '()))
298                (let ((fi1 (+ fi (f i))))
299                  (if (<= fi1 1.0)
300                      (recur (+ i dx) fi1 (cons fi1 lst))
301                      (reverse lst))))
302              )))
303
304
305        (define (compose-curves c1 c2)
306          (let ((c (pcurve:compose-curve (list + + +) c1 c2)))
307;            (print "compose-curves: c1 = ") (pp (pcurve:iterate-curve c1 10))
308;            (print "compose-curves: c2 = ") (pp (pcurve:iterate-curve c2 10))
309;            (print "compose-curves: c = ") (pp (pcurve:iterate-curve c 10))
310            c))
311             
312
313        (define (make-simple-curve fx fy fz n) 
314          (let ((c (pcurve:simple-curve n 1 (list fx fy fz) 0.0 1.0)))
315            c
316            ))
317       
318
319        (define (make-harmonic axis amp period phase n)
320          (let* ((freq (/ (* 2 PI) (/ 1.0 period)))
321                 (pf   (lambda (t) (* amp (cos (+ (* freq t) phase)))))
322                 (zf   (lambda (t) 0.0))
323                 (c (pcurve:simple-curve 
324                    (inexact->exact n) 1
325                    (cond ((= axis 0)
326                           (list pf zf zf))
327                          ((= axis 1)
328                           (list zf pf zf))
329                          ((= axis 2)
330                           (list pf zf zf))
331                          (else
332                           (error 'make-harmonic "invalid axis" axis))
333                          )
334                    0.0 1.0)))
335            c
336            ))
337       
338
339        (define (make-line-segment p dx dy dz) 
340          (let ((c (pcurve:line-segment 3 (list dx dy dz))))
341            (pcurve:translate-curve (point->list p) c)
342            ))
343       
344       
345        ;; auxiliary function to create segment points
346        (define (make-segment si np sp xyz)
347          (list-tabulate 
348           np
349           (lambda (i) (let* ((pi (+ sp i))
350                              (p (make-point 
351                                  (f64vector-ref (car xyz) pi)
352                                  (f64vector-ref (cadr xyz) pi)
353                                  (f64vector-ref (caddr xyz) pi))))
354                         (list si p)
355                         ))
356           ))
357
358
359         
360
361        ;;
362        ;; Creates a process of the given segments and number of points per
363        ;; segment from the given curve.
364        ;;
365        (define (make-segmented-process c f ns np)
366          (let ((xyz ((pcurve:sample-curve* c) (f (* ns np)))))
367            (let ((np1 (floor (/ (f64vector-length (car xyz)) ns))))
368              (let recur ((si ns) (ax '()))
369                (if (positive? si)
370                    (let ((sp (* (- si 1) np1)))
371                      (recur (- si 1) (append (make-segment si np1 sp xyz) ax)))
372                    ax)
373                ))
374            ))
375
376        ;;
377        ;; Non-segmented process
378        ;;
379        (define (make-process c f np)
380          (let ((xyz ((pcurve:sample-curve* c) (f np))))
381            (list-tabulate 
382             (f64vector-length (car xyz))
383             (lambda (i) 
384               (make-point 
385                (f64vector-ref (car xyz) i)
386                (f64vector-ref (cadr xyz) i)
387                (f64vector-ref (caddr xyz) i)))
388             ))
389          )
390       
391       
392        ;;
393        ;; Creates a named section containing points from the given segmented processes.
394        ;;
395        (define (make-segmented-section gid p label ps)
396          `(,label . 
397                   ,(fold (lambda (i+proc ax)
398                            (let ((seci (car i+proc)) 
399                                  (proc (cadr i+proc)))
400                              (cons
401                               `(,seci . 
402                                       ,(map (lambda (sp)
403                                               (let ((segi (car sp))
404                                                     (dpoint (cadr sp)))
405                                                 (let ((soma-distance (sqrt (dist2 p dpoint))))
406                                                   (make-genpoint dpoint gid soma-distance segi))
407                                                 ))
408                                             proc))
409                               ax)
410                              ))
411                          '() ps)
412                   ))
413       
414        ;;
415        ;; Non-segmented named section
416        ;;
417        (define (make-section gid p label ps)
418          `(,label . 
419                   ,(fold (lambda (i+proc ax)
420                            (let* ((seci (car i+proc)) 
421                                   (proc (cadr i+proc))
422                                   (pts (map (lambda (dpoint)
423                                               (let ((soma-distance (sqrt (dist2 p dpoint))))
424                                                 (make-genpoint dpoint gid soma-distance #f)))
425                                             proc)))
426                              (d "make-section: label = ~A gid = ~A p = ~A pts = ~A~%" 
427                                 label gid p proc)
428                              (cons `(,seci . ,pts) ax)
429                              ))
430                          '() ps)
431                   ))
432       
433       
434       
435        (define (make-gen-kd-tree data #!key (threshold 0.0))
436         
437          (define (update-bb p x-min y-min z-min x-max y-max z-max)
438            (let ((x (coord 0 p)) (y (coord 1 p)) (z (coord 2 p)))
439              (if (< x (x-min)) (x-min x))
440              (if (< y (y-min)) (y-min y))
441              (if (< z (z-min)) (z-min z))
442              (if (> x (x-max)) (x-max x))
443              (if (> y (y-max)) (y-max y))
444              (if (> z (z-max)) (z-max z))
445              ))
446
447
448          (let* (
449                 (t (list->kd-tree
450                     (filter (lambda (x) (<= threshold (cdr x))) data)
451                     make-value: (lambda (i v) (cdr v))
452                     make-point: (lambda (v) (apply make-point (car v)))
453                     offset: 0
454                     ))
455                 (n (kd-tree-size t))
456                 (x-min (make-parameter +inf.0))
457                 (y-min (make-parameter +inf.0))
458                 (z-min (make-parameter +inf.0))
459                 (x-max (make-parameter -inf.0))
460                 (y-max (make-parameter -inf.0))
461                 (z-max (make-parameter -inf.0))
462                 (bb (begin (kd-tree-for-each
463                             (lambda (p) (update-bb p x-min y-min z-min
464                                                    x-max y-max z-max))
465                             t)
466                            (list (x-min) (y-min) (z-min) (x-max) (y-max) (z-max))))
467                 )
468
469            (cons bb t)
470
471            ))
472
473
474
475
476
477        (define (genpoint-projection prefix my-comm myrank size cells fibers zone cell-start fiber-start)
478
479          (d "rank ~A: prefix = ~A zone = ~A~%" myrank prefix zone)
480
481          (fold (lambda (cell ax)
482
483                  (let* ((i (+ cell-start (car cell)))
484                         (root (modulo i size))
485                         (sections (cadr cell)))
486                   
487                    (fold 
488                     
489                     (lambda (sec ax)
490                       
491                       (let ((seci (car sec))
492                             (gxs  (cdr sec)))
493                         
494                         (let ((query-data
495                                (fold 
496                                 (lambda (gd ax)
497                                   (d "rank ~A: querying point ~A (~A)~%" 
498                                      myrank (genpoint-coords gd) 
499                                      (genpoint-parent-index gd))
500                                   (fold
501                                    (lambda (x ax) 
502                                      (let (
503                                            (source (car x))
504                                            (target i)
505                                            (distance (cadr x))
506                                            (segi (genpoint-segment gd))
507                                            )
508                                        (if (not segi)
509                                            (error 'genpoint-projection "point does not have segment index" gd))
510                                        (append (list source target distance segi seci) ax)
511                                        ))
512                                    ax
513                                    (delete-duplicates
514                                     (map (lambda (x) 
515                                            (d "rank ~A: query result = ~A (~A) (~A) ~%" 
516                                               myrank (kdnn-point x) (kdnn-distance x) (kdnn-parent-index x))
517                                            (list (+ fiber-start (kdnn-parent-index x))
518                                                  (+ (kdnn-distance x) (genpoint-parent-distance gd)  (kdnn-parent-distance x))
519                                                  ))
520                                          (kd-tree-near-neighbors* fibers zone (genpoint-coords gd)))
521                                     (lambda (u v)
522                                       (= (car u) (car v)))
523                                     )
524                                    ))
525                                 '() gxs)))
526                           (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
527
528                                  (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
529                             
530                             (append res1 ax))
531                           
532                           ))
533                       )
534                     ax sections)
535                    ))
536                '() (filter cadr cells)
537                ))
538
539
540       
541
542        (define (point-projection prefix my-comm myrank size pts fibers zone point-start nn-filter)
543          (fold (lambda (px ax)
544
545                  (d "~A: rank ~A: px = ~A zone=~A~%"  prefix myrank px zone)
546
547                  (let* ((i (+ point-start (car px)))
548                         (root (modulo i size))
549                         (dd (d "~A: rank ~A: querying point ~A (root ~A)~%" prefix myrank px root))
550                         (query-data
551                          (let ((pd (cadr px))) 
552                            (fold
553                             (lambda (x ax) 
554                               (let ((source (car x))
555                                     (target i)
556                                     (distance (cadr x)))
557                                 (if (and (> distance  0.) (not (= source target)))
558                                     (append (list source target distance) ax)
559                                     ax)
560                                 ))
561                             '()
562                             (delete-duplicates
563                              (map (lambda (x) 
564                                     (let ((res (list (car (cadar x)) 
565                                                      (+ (caddr x) (cadr (cadar x))))))
566                                       (d "~A: x = ~A res = ~A~%" prefix x res)
567                                       res))
568                                   (nn-filter pd (kd-tree-near-neighbors* fibers zone pd))
569                                   )
570                              (lambda (u v) (= (car u) (car v)))
571                              )
572                             ))
573                          ))
574
575
576                    (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
577                           (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
578                      (append res1 ax))
579                    ))
580
581                '() pts))
582       
583
584
585
586
587        (define-datatype layer-boundary layer-boundary?
588          (Bounds (top number?) (left number?) (bottom number?) (right number?))
589          (BoundsXZ (top number?) (left number?) (bottom number?) (right number?)
590                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
591          (BoundsYZ (top number?) (left number?) (bottom number?) (right number?)
592                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
593          )
594
595
596
597        (define (boundary-z-extent-function boundary)
598          (cases layer-boundary boundary
599                 (Bounds (top left bottom right) 
600                         (lambda (x y) 0.))
601                 (BoundsXZ (top left bottom right n k x y d d2) 
602                           (lambda (xp yp) 
603                             (let-values (((y0tab y1tab y2tab res)
604                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
605                               (f64vector-ref y0tab 0))))
606                 (BoundsYZ (top left bottom right n k x y d d2) 
607                           (lambda (xp yp) 
608                             (let-values (((y0tab y1tab y2tab res)
609                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
610                               (f64vector-ref y0tab 0))))
611                 ))
612
613
614        (define (point2d-rejection top left bottom right)
615            (lambda (p)
616              (let ((x (coord 0 p)) (y (coord 1 p)))
617                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
618            )
619
620
621        (define (compute-point3d p zu z-extent-function)
622          (let* ((x (coord 0 p))
623                 (y (coord 1 p))
624                 (z-extent (z-extent-function x y)))
625            (if (zero? zu)
626                (make-point x y zu)
627                (make-point x y (fp* zu z-extent)))
628            ))
629
630
631        (define (cluster-point-rejection p cluster-pts mean-distance)
632          (let ((D (* 4 mean-distance mean-distance))
633                (nn (kd-tree-nearest-neighbor cluster-pts p)))
634            (and (< (dist2 p nn) D) p)))
635
636
637
638        (define (XZAxis n k x-points z-points boundary)
639         
640          (if (not (>= n 3)) 
641              (error 'generate-boundary "XZAxis requires at least 3 interpolation points (n >= 3)"))
642         
643          (cases layer-boundary boundary
644                 (Bounds (top left bottom right) 
645
646                         (let-values (((d d2 constr errc diagn)
647                                       (bvsp-spline:compute n k x-points z-points)))
648                           
649                           (if (not (zero? errc)) 
650                               (error 'generate-boundary "error in constructing spline from boundary points" errc))
651                           
652                           (BoundsXZ top left bottom right n k x-points z-points d d2)))
653                 
654                 (else (error 'generate-boundary "boundary argument to XZAxis is already a pseudo-3D boundary")))
655          )
656
657
658        (define (Grid x-spacing y-spacing z-spacing boundary)
659
660          (match-let (((top left bottom right)
661                       (cases layer-boundary boundary
662                                   (Bounds (top left bottom right) 
663                                           (list top left bottom right))
664                                   (BoundsXZ (top left bottom right n k x y d d2) 
665                                             (list top left bottom right))
666                                   (BoundsYZ (top left bottom right n k x y d d2) 
667                                             (list top left bottom right))
668                                   )))
669
670          (let* (
671                 (x-extent   (- right left))
672                 (y-extent   (- top bottom))
673                 (z-extent-function
674                  (boundary-z-extent-function boundary))
675                 (compute-grid-points3d
676                  (lambda (p z-spacing z-extent-function)
677                    (let* ((x (coord 0 p))
678                           (y (coord 1 p))
679                           (z-extent (z-extent-function x y)))
680                      (let recur ((points '()) (z (/ z-spacing 2.)))
681                        (if (> z z-extent)
682                            points
683                            (recur (cons (make-point x y z) points) (+ z z-spacing)))
684                        ))
685                    ))
686                 )
687           
688            (d "Grid: x-spacing = ~A~%" x-spacing)
689            (d "Grid: y-spacing = ~A~%" y-spacing)
690            (d "Grid: z-spacing = ~A~%" z-spacing)
691           
692            (d "Grid: x-extent = ~A~%" x-extent)
693            (d "Grid: y-extent = ~A~%" y-extent)
694           
695            (let ((x-points (let recur ((points '()) (x (/ x-spacing 2.)))
696                              (if (> x x-extent)
697                                  (list->f64vector (reverse points))
698                                  (recur (cons x points) (+ x x-spacing)))))
699                  (y-points (let recur ((points '()) (y (/ y-spacing 2.)))
700                              (if (> y y-extent)
701                                  (list->f64vector (reverse points))
702                                  (recur (cons y points) (+ y y-spacing)))))
703                  )
704             
705              (let ((nx (f64vector-length x-points))
706                    (ny (f64vector-length y-points))
707                    )
708               
709                (let recur ((i 0) (j 0) (ax '()))
710                  (if (< i nx)
711                      (let ((x (f64vector-ref x-points i)))
712                        (if (< j ny)
713                            (let* ((y   (f64vector-ref y-points j))
714                                   (p   (make-point x y))
715                                   (p3ds (if (zero? z-spacing)
716                                             (list (make-point x y 0.0))
717                                             (compute-grid-points3d p z-spacing z-extent-function)))
718                                   )
719                              (recur i (+ 1 j) (if p3ds (append p3ds ax) ax)))
720                            (recur (+ 1 i) 0 ax)))
721                      (let* ((t (list->kd-tree ax))
722                             (n (kd-tree-size t)))
723                        (list t boundary)
724                        ))
725                  )))
726            ))
727          )
728
729        (define (UniformRandomPointProcess n x-seed y-seed boundary)
730
731          (match-let (((top left bottom right)
732                       (cases layer-boundary boundary
733                                   (Bounds (top left bottom right) 
734                                           (list top left bottom right))
735                                   (BoundsXZ (top left bottom right n k x y d d2) 
736                                             (list top left bottom right))
737                                   (BoundsYZ (top left bottom right n k x y d d2) 
738                                             (list top left bottom right))
739                                   )))
740
741          (let* (
742                 (x-extent   (- right left))
743                 (y-extent   (- top bottom))
744                 (z-extent-function (boundary-z-extent-function boundary))
745                 )
746
747            (let ((x-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init x-seed)))
748                  (y-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init y-seed)))
749                  (z-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init (current-milliseconds)))))
750             
751              (let ((point-rejection1 (point2d-rejection top left bottom right)))
752               
753                (let recur ((i 0) (ax '()))
754                  (if (< i n)
755                      (let ((x (f64vector-ref x-points i))
756                            (y (f64vector-ref y-points i))
757                            (z (f64vector-ref z-points i)))
758                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
759                          (let ((p3d 
760                                 (and (point-rejection1 p)
761                                      (compute-point3d p z z-extent-function))))
762                            (recur (+ 1 i) (if p3d (cons p3d ax) ax)))))
763                      (let* ((t (list->kd-tree ax))
764                             (n (kd-tree-size t)))
765
766                        (list t boundary))))
767                ))
768            ))
769          )
770
771
772        (define (ClusteredRandomPointProcess cluster-pts n mean-distance x-seed y-seed boundary)
773
774          (match-let (((top left bottom right)
775                       (cases layer-boundary boundary
776                                   (Bounds (top left bottom right) 
777                                           (list top left bottom right))
778                                   (BoundsXZ (top left bottom right n k x y d d2) 
779                                             (list top left bottom right))
780                                   (BoundsYZ (top left bottom right n k x y d d2) 
781                                             (list top left bottom right))
782                                   )))
783
784
785          (let* (
786                 (x-extent   (- right left))
787                 (y-extent   (- top bottom))
788                 (z-extent-function (boundary-z-extent-function boundary))
789                 )
790
791            (let recur ((pts '()) (x-seed x-seed) (y-seed y-seed))
792
793              (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
794                    (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
795                    (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
796               
797                (let ((point-rejection1 (point2d-rejection top left bottom right)))
798                 
799                  (let inner-recur ((j 0) (ax pts))
800                    (if (< j n)
801                        (let ((x (f64vector-ref x-points j))
802                              (y (f64vector-ref y-points j))
803                              (z (f64vector-ref z-points j)))
804                          (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
805                            (let ((p3d 
806                                   (and (point-rejection1 p)
807                                        (compute-point3d p z z-extent-function))))
808                              (let ((pp (cluster-point-rejection p3d cluster-pts mean-distance)))
809                                (inner-recur (+ 1 j)  (if pp (cons pp ax) ax))))))
810
811                        (if (< (length ax) n)
812                            (recur ax (+ 1 x-seed) (+ 1 y-seed))
813                            (let* ((t (list->kd-tree (take ax n)))
814                                   (n (kd-tree-size t)))
815                             
816                              (list t boundary))))
817                    ))
818                ))
819            ))
820          )
821
822
823
824       
825        (define comment-pat (string->irregex "^#.*"))
826
827
828        (define (load-points-from-file filename)
829
830          (let ((in (open-input-file filename)))
831           
832            (if (not in) (error 'load-points-from-file "file not found" filename))
833
834            (let* ((lines
835                    (filter (lambda (line) (not (irregex-match comment-pat line)))
836                            (read-lines in)))
837
838                   (point-data
839                    (filter-map
840                     (lambda (line) 
841                       (let ((lst (map string->number (string-split line " \t"))))
842                         (and (not (null? lst)) (apply make-point lst))))
843                     lines))
844
845                   (point-tree (list->kd-tree point-data))
846                   )
847             
848              (list point-tree)
849             
850              ))
851          )
852
853
854        (define (make-tree-graph lst label)
855         
856          (let* (
857                 (g              (make-digraph label #f))
858                 (node-info      (g 'node-info))
859                 (node-info-set! (g 'node-info-set!))
860                 (add-node!      (g 'add-node!))
861                 (add-edge!      (g 'add-edge!))
862                 )
863
864            ;; insert nodes
865            (let recur ((lst lst))
866
867              (if (not (null? lst))
868
869                  (let ((point (car lst)))
870
871                    (let ((node-id (swcpoint-id point)))
872
873                      (add-node! node-id point)
874
875                      (recur (cdr lst))))))
876
877            ;; insert edges
878            (let recur ((lst lst))
879             
880              (if (not (null? lst))
881                 
882                  (let ((point (car lst)))
883                   
884
885                    (let ((node-id (swcpoint-id point))
886                          (pre-id  (swcpoint-pre point)))
887
888                      (if (> pre-id 0)
889
890                          (let* ((pre-point   (node-info pre-id))
891                                 (pre-coords  (and pre-point (swcpoint-coords pre-point)))
892                                 (node-coords (swcpoint-coords point))
893                                 (distance    (sqrt (dist2 node-coords pre-coords))))
894                           
895                            (add-edge! (list pre-id node-id distance))))
896                       
897                      (recur (cdr lst))
898                      ))
899                  ))
900            g 
901            ))
902
903
904        (define (tree-graph-distances+segments g nseg)
905
906
907          (define n        ((g 'capacity)))
908          (define distv    (make-f64vector (+ 1 n) -1.0))
909          (define rdistv   (make-f64vector (+ 1 n) -1.0))
910          (define segv     (make-s32vector (+ 1 n) -1))
911
912
913          ;; determine distances from root
914          (define (traverse-dist es)
915            (if (null? es) distv
916                (match-let (((i j dist) (car es)))
917                  (if (>= (f64vector-ref distv j) 0.0)
918                      (traverse-dist (cdr es))
919                      (let ((idist (f64vector-ref distv i)))
920                        (f64vector-set! distv j (+ idist dist))
921                        (let ((distv1 (traverse-dist ((g 'out-edges) j))))
922                          (traverse-dist es)))
923                      ))
924                ))
925         
926         
927          ;; determine distances from end (reverse distance)
928          (define (traverse-rdist es)
929            (if (null? es) rdistv
930                (match-let (((i j dist) (car es)))
931                  (if (>= (f64vector-ref rdistv i) 0.0)
932                      (traverse-rdist (cdr es))
933                      (let ((jdist (f64vector-ref distv j)))
934                        (f64vector-set! rdistv i (+ jdist dist))
935                        (let ((rdistv1 (traverse-rdist ((g 'in-edges) i))))
936                          (traverse-rdist es)))
937                      ))
938                ))
939
940
941          (define (compute-segv distv rdistv)
942            (let recur ((n n))
943              (if (>= n 1)
944                  (let* ((dist  (f64vector-ref distv n))
945                         (rdist (f64vector-ref rdistv n))
946                         (len   (and (positive? dist) (positive? rdist) (+ dist rdist)))
947                         (delta (and len (round (/ len nseg))))
948                         (seg   (and delta (round (/ dist delta)))))
949                    (if seg (s32vector-set! segv n (exact->inexact seg)))
950                    (recur (- n 1))
951                    ))
952              ))
953         
954          (let ((in-edges (g 'in-edges)) 
955                (out-edges (g 'out-edges)) 
956                (terminals ((g 'terminals)))
957                (roots ((g 'roots))))
958            (for-each (lambda (x) (f64vector-set! distv x 0.0)) roots)
959            (for-each (lambda (x) (s32vector-set! segv x 0)) roots)
960            (for-each (lambda (x) (f64vector-set! rdistv x 0.0)) terminals)
961            (traverse-dist (concatenate (map (lambda (x) (out-edges x)) roots)))
962            (traverse-rdist (concatenate (map (lambda (x) (in-edges x)) terminals)))
963            (compute-segv distv rdistv)
964            (list distv segv)
965          ))
966
967
968        (define (tree-graph->section-points cell-index cell-origin type g gdistv gsegv)
969         
970          (let* ((node-info (g 'node-info))
971                 (succ      (g 'succ))
972                 (offset    (let ((cell-loc (point->list cell-origin))
973                                  (root-loc (point->list (swcpoint-coords (node-info 1)))))
974                              (map - cell-loc root-loc))))
975
976            (d "tree-graph->section-points: offset = ~A~%" offset)
977
978
979            (let recur ((n 1) (lst '()))
980
981              (let* (
982                     (point (node-info n))
983                     (point-type (swcpoint-type point))
984                     (point-pre (swcpoint-pre point))
985                     (proceed? (or (= point-type type)
986                                   (case (swcpoint-type point)
987                                     ((0 1 5 6) #t)
988                                     (else #f))))
989                     )
990
991                (d "tree-graph->section-points: n = ~A point-type = ~A proceed? = ~A~%" 
992                   n point-type proceed?)
993
994                (d "tree-graph->section-points: succ n = ~A~%" (succ n))
995                 
996                (if proceed?
997
998                    (let (
999                          (point1 (list
1000                                   (s32vector-ref gsegv n)
1001                                   (apply make-point (map + offset (point->list (swcpoint-coords point))))))
1002                          )
1003
1004                      (fold (lambda (x ax) (recur x ax))
1005                            (cons point1 lst)
1006                            (succ n)))
1007
1008                    lst)
1009
1010                ))
1011            ))
1012
1013 
1014        (define (load-swc filename label type nseg)
1015         
1016          (let ((in (open-input-file filename)))
1017           
1018            (if (not in) (error 'load-swc "file not found" filename))
1019           
1020            (let* (
1021                   (lines
1022                    (let ((lines (read-lines in)))
1023                      (close-input-port in)
1024                      (filter (lambda (line) (not (irregex-match comment-pat line))) 
1025                              lines)))
1026
1027                   (swc-data
1028                    (filter-map
1029                     (lambda (line) 
1030                       (let ((lst (map string->number (string-split line " \t"))))
1031                         (and (not (null? lst)) 
1032                              (match-let (((id my-type x y z radius pre) lst))
1033                                         (make-swcpoint id my-type (make-point x y z)
1034                                                        radius pre)))
1035                         ))
1036                     lines))
1037
1038                   (swc-graph (make-tree-graph swc-data label))
1039
1040                   (dist+segs  (tree-graph-distances+segments swc-graph nseg))
1041
1042                   )
1043
1044              (cons type (cons swc-graph dist+segs)))
1045          ))
1046
1047
1048        (define (load-swcdir path label type nseg)
1049         
1050          (let ((pat ".*.swc"))
1051
1052            (let ((flst (find-files path
1053                                    test: (regexp pat)
1054                                    action: cons
1055                                    seed: (list) 
1056                                    limit: 0)))
1057
1058              (d "load-swcdir: flst = ~A~%" flst)
1059
1060              (map (lambda (fn) (load-swc fn label type nseg)) (sort flst string<?))
1061              ))
1062          )
1063         
1064
1065        (define (segment-projection label source-tree target-sections zone my-comm myrank size)
1066
1067          (MPI:barrier my-comm)
1068         
1069          (let ((my-results
1070                 (genpoint-projection label my-comm myrank size target-sections source-tree zone 0 0)))
1071
1072            (MPI:barrier my-comm)
1073
1074            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
1075              (lambda (out-sources)
1076                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
1077                  (lambda (out-targets)
1078                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
1079                      (lambda (out-distances)
1080                        (call-with-output-file (sprintf "~Asegments~A.dat"  label (if (> size 1) myrank ""))
1081                          (lambda (out-segments)
1082                            (for-each
1083                             (lambda (my-data)
1084                               (let* ((my-entry-len 5)
1085                                      (my-data-len (/ (f64vector-length my-data) my-entry-len)))
1086                                 (d "rank ~A: length my-data = ~A~%" myrank my-data-len)
1087                                 (let recur ((m 0))
1088                                   (if (< m my-data-len)
1089                                       (let* (
1090                                              (my-entry-offset (* m my-entry-len))
1091                                              (source (inexact->exact (f64vector-ref my-data my-entry-offset)))
1092                                              (target (inexact->exact (f64vector-ref my-data (+ 1 my-entry-offset))))
1093                                              (distance (f64vector-ref my-data (+ 2 my-entry-offset)))
1094                                              (segment (f64vector-ref my-data (+ 3 my-entry-offset)))
1095                                              (section (f64vector-ref my-data (+ 4 my-entry-offset)))
1096                                              )
1097                                         (fprintf out-sources   "~A~%" source)
1098                                         (fprintf out-targets   "~A~%" target)
1099                                         (fprintf out-distances "~A~%" distance)
1100                                         (fprintf out-segments  "~A ~A~%" segment section)
1101                                         (recur (+ 1 m)))))
1102                                 ))
1103                             my-results))
1104                          ))
1105                      ))
1106                  ))
1107              ))
1108          )
1109       
1110
1111        (define (projection label source-tree target zone my-comm myrank size) 
1112
1113          (MPI:barrier my-comm)
1114         
1115          (let ((my-results (point-projection label my-comm myrank size target source-tree zone 0 (lambda (x nn) nn))))
1116
1117          (MPI:barrier my-comm)
1118           
1119            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
1120              (lambda (out-sources)
1121                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
1122                  (lambda (out-targets)
1123                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
1124                      (lambda (out-distances)
1125                        (for-each
1126                         (lambda (my-data)
1127                           (let* ((my-entry-len 3)
1128                                  (my-data-len (/ (f64vector-length my-data) my-entry-len)))
1129                             (d "~A: rank ~A: length my-data = ~A~%" label myrank my-data-len)
1130                             (let recur ((m 0))
1131                               (if (< m my-data-len)
1132                                   (let ((source (inexact->exact (f64vector-ref my-data (* m my-entry-len))))
1133                                         (target (inexact->exact (f64vector-ref my-data (+ 1 (* m my-entry-len)))))
1134                                         (distance (f64vector-ref my-data (+ 2 (* m my-entry-len)))))
1135                                     (fprintf out-sources "~A~%" source)
1136                                     (fprintf out-targets "~A~%" target)
1137                                     (fprintf out-distances "~A~%" distance)
1138                                     (recur (+ 1 m)))))
1139                             ))
1140                         my-results)
1141                        ))
1142                    ))
1143                ))
1144            ))
1145       
1146
1147        (include "calc-parser.scm")
1148
1149       
1150        (define (load-config-file filename)
1151          (let ((in (open-input-file filename)))
1152            (if (not in) (error 'load-config-file "file not found" filename))
1153            (init-bindings)
1154            (let* ((lines (reverse (filter (lambda (x) (not (string-null? x))) (read-lines in))))
1155                   (properties (filter-map
1156                                (lambda (line) 
1157                                  (and (not (string-prefix? "//" line))
1158                                       (let ((lst (string-split line "=")))
1159                                         (and (> (length lst) 1)
1160                                              (let ((key (string->symbol (string-trim-both (car lst) #\space)))
1161                                                    (val (calc-eval-string (cadr lst))))
1162                                                (add-binding key val)
1163                                                (cons key val))))))
1164                                lines)))
1165              properties
1166              ))
1167          )
1168
1169        (define (make-output-fname dirname sysname suffix . rest) 
1170          (let-optionals 
1171           rest ((x #t))
1172           (and x
1173                (let ((dirname (if (string? x) x dirname)))
1174                  (let ((fname (sprintf "~A~A" sysname suffix)))
1175                    (or (and dirname (make-pathname dirname fname)) fname)))
1176                )))
1177
1178       
1179)
1180
1181       
Note: See TracBrowser for help on using the repository browser.