source: project/release/4/neurolucida/trunk/neurolucida.scm @ 30623

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

neurolucida: removed usage of system*

File size: 24.8 KB
Line 
1
2;;
3;; Neurolucida XML file manipulation.
4;;
5;; Copyright 2011-2012 Ivan Raikov and the Okinawa Institute of Science and
6;; 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(require-extension files data-structures srfi-1 srfi-13  getopt-long sxml-transforms sxpath sxpath-lolevel ssax)
23(require-extension typeclass kd-tree iset digraph format-graph  )
24
25(import-instance (<KdTree> KdTree3d)
26                 (<Point> Point3d))
27
28
29(define lookup-def 
30  (lambda (k lst . rest)
31    (let-optionals rest ((default #f))
32      (alist-ref k lst eq? default))))
33
34
35(define (quotewrapped? str)
36  (and (string? str) (string-prefix? "\"" str) (string-suffix? "\"" str) ))
37
38
39(define (quotewrap str)
40  (cond ((quotewrapped? str) str)
41        ((string-any char-whitespace? str)
42         (string-append "\"" str "\""))
43        (else str)))
44
45
46; Returns attr-list node for a given obj
47;   or #f if it is absent
48(define (sxml:attr-alist obj)
49  (if (and (not (null? (cdr obj)))
50            (pair? (cadr obj)) 
51            (eq? '@ (caadr obj)))
52         (cdr (cadr obj))
53         #f))
54
55;; obtain  child named n of a node
56(define (sxml:kidn name node)
57  ((select-first-kid (lambda (x)  (eq? (car x) name))) node))
58
59;; obtain all children of a node named n
60(define (sxml:kidsn name node)
61  ((select-kids (lambda (x) (eq? (car x) name))) node))
62
63(define (mkdir dir)
64  (system (sprintf "mkdir -p ~a" (quotewrap dir))))
65
66
67(define *quiet* #f)
68
69
70(define (d fstr . args)
71  (let ([port (if *quiet* (current-error-port) (current-output-port))])
72    (apply fprintf port fstr args)
73    (flush-output port) ) )
74
75
76(define opt-defaults
77  `(
78    (key             . "color")
79    (format          . "swc")
80    (forest-grid     . "d:10,p:5,n:1")
81    ))
82
83
84(define (defopt x)
85  (lookup-def x opt-defaults))
86
87(define (symbol-upcase str)
88  (string->symbol (string-upcase str)))
89
90
91(define opt-grammar
92  `(
93    (forest-grid
94     "print grid coordinates covering the entire forest (resolution is given in number of points)"
95     (value       (optional "P:N[,D:N,N:1]")
96                  (transformer ,(compose (lambda (xs) 
97                                           (map (lambda (x) 
98                                                  (let ((kv (string-split x ":")))
99                                                    (cons (string->symbol (string-downcase (car kv)))
100                                                          (string->number (cadr kv))))) xs ))
101                                         (lambda (x) (string-split x ","))))
102                  (default  ,(defopt 'forest-grid))))
103
104    (permute-coords
105     "permute coordinates of node points in NG format according to the given list of 1-based indices"
106     (value (required INDICES)))
107
108    (data-dir
109     "set download directory (default is a randomly generated name in /tmp)"
110     (single-char #\d)
111     (value (required DIR)))
112
113    (key
114     "extraction key"
115     (single-char #\k)
116     (value       (required "FIELD")
117                  (default  ,(defopt 'key))))
118
119    (format
120     "output format (swc, pts, vcg, ng)"
121     (single-char #\f)
122     (value       (required "FORMAT")
123                  (default  ,(defopt 'format))))
124
125    (version
126     "print the current version of neurolucida and exit")
127
128    (help  "Print help"
129            (single-char #\h))
130 
131  ))
132
133
134;; Use args:usage to generate a formatted list of options (from OPTS),
135;; suitable for embedding into help text.
136(define (neurolucida:usage)
137  (print "Usage: " (car (argv)) " [options...] operands ")
138  (newline)
139  (print "Where operands are XML Neurolucida files")
140  (newline)
141  (print "The following options are recognized: ")
142  (newline)
143  (width 35)
144  (print (parameterize ((indent 5)) (usage opt-grammar)))
145  (exit 1))
146
147
148;; Process arguments and collate options and arguments into OPTIONS
149;; alist, and operands (filenames) into OPERANDS.  You can handle
150;; options as they are processed, or afterwards.
151
152(define opts    (getopt-long (command-line-arguments) opt-grammar))
153(define opt     (make-option-dispatch opts opt-grammar))
154
155
156(define data-dir (make-parameter #f))
157
158(define (get-data-dir)
159  (or (opt 'data-dir)
160      (or (data-dir)
161          (let ([dir (create-temporary-directory)])
162            (data-dir dir)
163            dir ) ) ))
164
165
166(define (create-temporary-directory)
167  (let ((dir (or (get-environment-variable "TMPDIR") 
168                 (get-environment-variable "TEMP") 
169                 (get-environment-variable "TMP") 
170                 "/tmp")))
171    (let loop ()
172      (let* ((n (current-milliseconds)) ; current milliseconds
173             (pn (make-pathname dir (string-append "neurolucida-" (number->string n 16)) "tmp")))
174        (cond ((file-exists? pn) (loop))
175              (else (mkdir pn) pn))))))
176
177
178(define neurolucida-xmlns "http://www.mbfbioscience.com/2007/neurolucida")
179
180
181(define (parse-sxml fpath)
182  (with-input-from-file fpath
183    (lambda () (cons '*TOP* (ssax:xml->sxml (current-input-port) `((nl . ,neurolucida-xmlns)))))))
184
185
186(define (extract-element sxml element-name key-attr)
187  (let ((element-fullname (string->symbol (conc "nl:" (->string element-name)))))
188    (let* ((elements ((sxpath `(// nl:mbf ,element-fullname ))  sxml))
189           (keys     ((sxpath `(// ,element-fullname @ ,key-attr)) `(*TOP* ,@elements))))
190      (zip keys elements))))
191
192
193(define (key=? a b)  (string=? (cadar a) (cadar b)))
194
195
196(define (partition-by-key a k key=?)
197  (partition (lambda (x) (equal? (car x) k)) a))
198
199
200(define (get-node-origin x)
201  (let ((x (if (equal? 'nl:point (car x)) x (sxml:kidn 'nl:point x))))
202    (let ((attrs (sxml:attr-alist x)))
203      (list (string->number (car (lookup-def 'x attrs)))
204            (string->number (car (lookup-def 'y attrs)))
205            (string->number (car (lookup-def 'z attrs))))
206      )))
207
208
209(define (get-node-radius x)
210  (let ((x (if (equal? 'nl:point (car x)) x (sxml:kidn 'nl:point x))))
211    (let ((attrs (sxml:attr-alist x)))
212      (/ (string->number (car (lookup-def 'd attrs))) 2.))))
213
214
215(define (make-tree-graph label x)
216
217  (let* ((sxml           `(*TOP* ,@(cdr x)))
218         (g              (make-digraph label #f))
219         (node-info      (g 'node-info))
220         (node-info-set! (g 'node-info-set!))
221         (add-node!      (g 'add-node!))
222         (add-edge!      (g 'add-edge!))
223         )
224
225    (add-node! 0 'soma)
226
227    (let ((initial-tree ((sxpath `(nl:tree))  sxml)))
228
229      (let recur ((tree      initial-tree)
230                   (parent   0)
231                   (node-id  1))
232
233      (let ((points+spines     ((sxpath `((*or* nl:point nl:spine))) tree))
234            (subtrees          (append ((sxpath `(nl:tree)) tree)
235                                       ((sxpath `(nl:branch)) tree))))
236
237        ;; insert points in the dependency graph
238        (let ((node-id1
239               (car
240               (fold (lambda (n i.lp) 
241                       (case (car n)
242                         ((nl:point)
243                          (let ((i (car i.lp))
244                                (lp (cdr i.lp)))
245
246                            (if (not (equal? lp n))
247                                (begin
248                                  (add-node! i n) (cons (+ 1 i) n))
249                                i.lp)))
250                         ((nl:spine)
251                          (let ((i (car i.lp)) (lp (cdr i.lp)))
252                            (node-info-set! (- i 1) (append (node-info (- i 1)) (list n))) 
253                            i.lp))
254                         (else (error 'make-tree-graph "unknown node type" n))))
255                     (cons node-id  #f)
256                     points+spines))))
257         
258          ;; insert edges
259          (if (pair? points+spines)
260              (let ((fin (- node-id1 1)))
261                (add-edge! (list parent node-id #f))
262                (let inner-recur ((i node-id))
263                  (if (< i fin)
264                      (let ((j (+ 1 i))) 
265                        (add-edge! (list i j  #f))
266                        (inner-recur j))
267                      ))
268                ))
269         
270          (if (pair? subtrees)
271              (let ((branch-node (- node-id1 1)))
272                (node-info-set! branch-node `(nl:branch ,(node-info branch-node)))
273                (fold (lambda (subtree i) 
274                        (recur subtree branch-node (+ i 1))) node-id1 subtrees))
275              node-id1)
276          ))
277      ))
278    g))
279
280
281
282;;
283;; NG (Neurolucida Graph) format
284;;
285;; ((Branch ... (Node ...) (Branch ...) ) ... )
286;;
287;; Where ... is one or more attributes:
288;;
289;;   origin x y z
290;;   radius r
291;;   length l
292;;   spine-density-linear
293;;   branch-order-Soma
294;;   branch-order-Terminal
295;;
296
297(define (make-ng key g #!key (permute-coords #f))
298
299  (define (sum lst) (fold + 0. lst))
300
301
302  (define (pdist2 a b)
303    (let ((diff2 (lambda (i) (let ((v (- (list-ref a i) (list-ref b i)))) (* v v)))))
304      (sum (list-tabulate 3 diff2))))
305 
306
307  ;; DFS that increases depth only when branches are encountered
308  (define (branch-dfs-fold g fn fb roots x y 
309                           #!key
310                           (edge-target cadr)
311                           (next-edges (g 'out-edges)))
312   
313
314    (define (traverse visited n x y)
315      (if (bit-vector-ref visited n)
316          (values visited x y)
317          (let ((visited (bit-vector-set! visited n #t))
318                (x1 (fn n x y)))
319            (traverse-edges visited (next-edges n) x1 y))
320          ))
321
322   
323    (define (traverse-edges visited elst x y) =
324      (if (null? elst)
325          (values visited x y)
326          (let ((n (edge-target (car elst)))
327                (es (cdr elst)))
328            (if (bit-vector-ref visited n)
329                (traverse-edges visited es x y)
330                (let ((visited (bit-vector-set! visited n #t))
331                      (x (fn n x y)))
332                  (let ((enext (next-edges n))
333                        (out-edges ((g 'out-edges) n)))
334                    (let ((y (if (> (length out-edges) 1)
335                                 (fb n out-edges x y) y)))
336                      (let-values (((visited x y) 
337                                    (traverse-edges visited enext x y)))
338                        (traverse-edges visited es x y)
339                        ))
340                    ))
341                ))
342          ))
343
344   
345    (define (traverse-roots visited ns x y)
346      (if (null? ns) (values visited x y)
347          (let-values (((visited x y) (traverse visited (car ns) x y)))
348            (traverse-roots visited (cdr ns) x y))))
349
350   
351    (traverse-roots (make-bit-vector ((g 'capacity))) roots x y)
352    )
353
354
355   
356  ;;
357  ;; Compute node order from the soma and from the terminals
358  ;;
359 
360  (define (compute-node-orders g)
361     
362    (let ((roots ((g 'roots)))
363          (terminals  ((g 'terminals))))
364
365      (let-values (
366                   ((visited node-branch-order final-branch-depth) ;; branch order from soma
367                    (branch-dfs-fold g 
368                                     (lambda (n x y) (s32vector-set! x n y) x) 
369                                     (lambda (n es x y) (+ 1 y))
370                                     roots
371                                     (make-s32vector ((g 'capacity)) -1) ;; branch order per node
372                                     0 ;; branch order
373                                     ))
374                   ((visited node-inverse-branch-order final-inverse-branch-depth) ;; branch order from terminals
375                    (branch-dfs-fold g 
376                                     (lambda (n x y) (s32vector-set! x n y) x) 
377                                     (lambda (n es x y) (+ 1 y))
378                                     terminals
379                                     (make-s32vector ((g 'capacity)) -1) ;; branch order per node
380                                     0 ;; branch order
381                                     edge-target: car
382                                     next-edges: (g 'in-edges)
383                                     ))
384                   )
385
386        (let (
387              ;; branch order from soma
388              (node-depth node-branch-order)
389              ;; branch order from terminals
390              (inverse-node-depth node-inverse-branch-order)
391              )
392       
393          `((node-depth . ,node-depth)
394            (inverse-node-depth . ,inverse-node-depth)
395            ))
396          )))
397
398  (define (compute-node-tree roots out-edges node-info 
399                             node-orders node-depths inverse-node-depths)
400
401    (let recur ((node (car roots)) (ax '()))
402             
403      (let* ((info   (node-info node))
404             (spines (sxml:kidsn 'nl:spine info))
405             )
406       
407       
408        (if (equal? info 'soma)
409           
410            (fold recur '() (map cadr (out-edges node)))
411           
412            (case (car info)
413             
414              ((nl:point)
415               (let* ((oes (out-edges node))
416                      (next-node (and (pair? oes) (cadr (car oes))))
417                      (cylinder-origin (get-node-origin info))
418                      (cylinder-radius (get-node-radius info))
419                      (cylinder-length (and next-node
420                                            (let ((node1-origin (get-node-origin (node-info next-node))))
421                                              (sqrt (pdist2 cylinder-origin node1-origin))
422                                              ))
423                                       )
424                      )
425                 
426                 (assert (<= (length oes) 1))
427                 (if cylinder-length
428                     (if (not (positive? cylinder-length))
429                         (error 'neurolucida "zero cylinder length" info (node-info next-node))))
430                 
431                 (let ((ax1 (cons
432                             (if next-node
433                                 `(Node (origin . ,(or (and permute-coords (permute-coords cylinder-origin) )
434                                                       cylinder-origin))
435                                        (radius . ,cylinder-radius)
436                                        (length . ,cylinder-length)
437                                        (branch-order-Soma . ,(s32vector-ref node-depths node))
438                                        (branch-order-Terminal . ,(s32vector-ref inverse-node-depths node))
439                                        ,@(if (null? spines) '() `((spine-density-linear . ,(/ (length spines) cylinder-length))))
440                                        )
441                                 `(Terminal (origin . ,(or (and permute-coords (permute-coords cylinder-origin) )
442                                                           cylinder-origin))
443                                            (radius . ,cylinder-radius)
444                                            (branch-order-Soma . ,(s32vector-ref node-depths node))
445                                            (branch-order-Terminal . ,(s32vector-ref inverse-node-depths node))
446                                            ))
447                             ax)))
448                   
449                   (if next-node (recur next-node ax1) (reverse ax1))
450                   
451                   )
452                 ))
453             
454              ((nl:branch)
455               (let* ((branch-origin (get-node-origin info))
456                      (branch-radius (get-node-radius info))
457                      )
458                 
459                 (cons
460                  `(Branch (origin . ,(or (and permute-coords (permute-coords branch-origin) )
461                                          branch-origin))
462                           (radius . ,branch-radius)
463                           (branch-order-Soma . ,(s32vector-ref node-depths node))
464                           (branch-order-Terminal . ,(s32vector-ref inverse-node-depths node))
465                           (children . ,(map (lambda (x) (recur x '())) (map cadr (out-edges node))))
466                           )
467                  ax)
468                 ))
469             
470              (else (recur (cdr nodes) ax)))
471           
472            ))
473      ))
474 
475  (let* ((roots ((g 'roots)))
476         (out-edges (g 'out-edges))
477         (node-info (g 'node-info))
478         (node-orders (compute-node-orders g))
479         (node-depths (lookup-def 'node-depth node-orders))
480         (inverse-node-depths (lookup-def 'inverse-node-depth node-orders))
481         )
482
483
484    (let ((tree
485          `(((key . ,key)) . 
486            ,(compute-node-tree roots out-edges node-info
487                                node-orders node-depths inverse-node-depths))))
488
489      tree
490      )))
491     
492 
493
494
495
496;;
497;; SWC format
498;;
499;; n T x y z R P
500;;
501;; n is an integer label that identifies the current point and
502;; increments by one from one line to the next.
503;;
504;; T is an integer representing the type of neuronal segment, such as
505;; soma, axon, apical dendrite, etc. The standard accepted integer
506;; values are given below.
507;;
508;;     * 0 = undefined
509;;     * 1 = soma
510;;     * 2 = axon
511;;     * 3 = dendrite
512;;     * 4 = apical dendrite
513;;     * 5 = fork point
514;;     * 6 = end point
515;;     * 7 = custom
516;;
517;; x, y, z gives the cartesian coordinates of each node.
518;;
519;; R is the radius at that node.
520;;
521;; P indicates the parent (the integer label) of the current point or
522;; -1 to indicate an origin (soma).
523
524
525(define (print-swc x)
526  (define get (compose car alist-ref))
527
528  (let ((sxml `(*TOP* ,@(cdr x))))
529
530    (let* ((offset (make-parameter 0))
531           (soma ((sxpath `(nl:contour))  sxml))
532           (soma-points ((sxpath `(nl:point @))  soma))
533           (soma-npoints (length soma-points))
534           (soma-indices (list-tabulate soma-npoints (lambda (x) (+ 1 x)))))
535
536       (for-each (lambda (p n)
537                  (let ((data (cdr p)))
538                    (let ((T    1) ;; soma
539                          (x  (get 'x data))
540                          (y  (get 'y data))
541                          (z  (get 'z data))
542                          (R  (number->string (/ (string->number (get 'd data)) 2)))
543                          (P  (- n 1))
544                          )
545                    (printf "~A ~A ~A ~A ~A ~A ~A~%" n T x y z R (if (zero? P) -1 P)))))
546                soma-points soma-indices)
547       
548       (offset soma-npoints)
549
550       (let recur ((tree ((sxpath `(nl:tree))  sxml)) 
551                   (parent-index soma-npoints))
552
553         (let ((points     ((sxpath `(nl:point @))  tree))
554               (subtrees   (append ((sxpath `(nl:tree))  tree)
555                                   ((sxpath `(nl:branch))  tree))))
556           
557           (let* ((npoints (length points))
558                  (indices (let ((k (+ 1 (offset)))) (list-tabulate npoints (lambda (x) (+ k x)))))
559                  (parents (cons parent-index indices)))
560             
561             (fold (lambda (pt i par n)
562                         (let ((data (cdr pt)))
563                           (let ((T    (if (fx> n 0) 3 (if (null? subtrees) 6 5))) ;; dendrite, or end point/bifurcation point
564                                 (x  (get 'x data))
565                                 (y  (get 'y data))
566                                 (z  (get 'z data))
567                                 (R  (number->string (/ (string->number (get 'd data)) 2)))
568                                 (P  par)
569                                 )
570                             (printf "~A ~A ~A ~A ~A ~A ~A~%" 
571                                     i T x y z R P)
572                             (fx- n 1)
573                             )))
574                   (fx- npoints 1)
575                   points indices parents)
576             
577             (let ((parent-index1 (+ npoints (offset))))
578
579               (offset (+ (offset) npoints)) 
580
581               (if (pair? subtrees)
582                   (for-each (lambda (subtree) (recur subtree parent-index1) )
583                             subtrees)))
584             )))
585
586    )))
587
588
589(define (print-pts x)
590  (define get (compose car alist-ref))
591
592  (let ((sxml `(*TOP* ,@(cdr x))))
593
594    (let* ((soma ((sxpath `(nl:contour))  sxml))
595           (soma-points ((sxpath `(nl:point @))  soma))
596           (soma-npoints (length soma-points)))
597
598       (for-each (lambda (p)
599                  (let ((data (cdr p)))
600                    (let ((x  (get 'x data))
601                          (y  (get 'y data))
602                          (z  (get 'z data))
603                          (R  (/ (string->number (get 'd data)) 2))
604                          )
605                      (printf "~A ~A ~A~%" x y z))))
606                soma-points)
607
608       (let recur ((tree ((sxpath `(nl:tree))  sxml))  )
609
610         (let ((points     ((sxpath `(nl:point @))  tree))
611               (subtrees   (append ((sxpath `(nl:tree))  tree)
612                                   ((sxpath `(nl:branch))  tree))))
613           
614           (let ((npoints (length points)))
615             
616             (for-each (lambda (pt)
617                         (let ((data (cdr pt)))
618                           (let ((x  (string->number (get 'x data)))
619                                 (y  (string->number (get 'y data)))
620                                 (z  (string->number (get 'z data)))
621                                 (r  (/ (string->number (get 'd data)) 2))
622                                 )
623                             (printf "~A ~A ~A~%" x y z)
624                             )))
625                       points)
626             
627               (if (pair? subtrees)
628                   (for-each (lambda (subtree) (recur subtree) )
629                             subtrees))
630             )))
631
632    )))
633
634
635(define (print-grid data-map options)
636
637  (define get (compose car alist-ref))
638
639  (let ((nbndpts (alist-ref 'p options))
640        (b-radius (alist-ref 'd options))
641        (min-neighbors (alist-ref 'n options))
642        (bb.all-points
643         (fold 
644          (lambda (x ax)
645
646            (let ((sxml `(*TOP* ,@(cdr x))))
647             
648              (let recur ((trees ((sxpath `(nl:tree))  sxml))
649                          (exs (car ax))
650                          (all-points (cdr ax)))
651               
652                (if (null? trees) (cons exs all-points)
653                   
654                    (let ((points     ((sxpath `(nl:point @))  (car trees)))
655                          (subtrees   (append ((sxpath `(nl:tree))  (car trees))
656                                              ((sxpath `(nl:branch))  (car trees)))))
657                     
658                      (let (
659                            (exs.all-points
660                             (fold (lambda (pt ax)
661
662                                     (let ((data (cdr pt)))
663
664                                       (let ((x  (string->number (get 'x data)))
665                                             (y  (string->number (get 'y data)))
666                                             (z  (string->number (get 'z data))))
667                                         
668                                         (let* ((exs (car ax))
669                                                (all-points (cdr ax))
670                                                (x-min (list-ref exs 0))
671                                                (x-max (list-ref exs 1))
672                                                (y-min (list-ref exs 2))
673                                                (y-max (list-ref exs 3))
674                                                (z-min (list-ref exs 4))
675                                                (z-max (list-ref exs 5)))
676                                           
677                                           (cons
678                                            (list (min x x-min) (max x x-max)
679                                                  (min y y-min) (max y y-max)
680                                                  (min z z-min) (max z z-max))
681                                            (cons (make-point x y z) all-points))
682                                           ))
683                                       ))
684                                   (cons exs all-points) points))
685                            )
686                       
687                        (let ((exs.all-points
688                               (fold (lambda (subtree ax) (recur subtree (car ax) (cdr ax )))
689                                     exs.all-points subtrees)))
690                          (recur (cdr trees) (car exs.all-points) (cdr exs.all-points)))
691                        ))
692                    ))
693              ))
694          (cons (list +inf.0 -inf.0
695                      +inf.0 -inf.0
696                      +inf.0 -inf.0) 
697                '())
698          data-map))
699        )
700   
701      (if (null? (cdr bb.all-points))
702          (error "empty point set in forest"))
703
704      (let* (
705             (bb (car bb.all-points))
706             (x-min (list-ref bb 0))
707             (x-max (list-ref bb 1))
708             (y-min (list-ref bb 2))
709             (y-max (list-ref bb 3))
710             (z-min (list-ref bb 4))
711             (z-max (list-ref bb 5))
712             (point-tree (list->kd-tree (cdr bb.all-points)))
713             )
714       
715        (let* ((x-extent (- x-max x-min))
716               (y-extent (- y-max y-min))
717               (z-extent (- z-max z-min))
718               (x-step (ceiling (abs (/ x-extent nbndpts))))
719               (y-step (ceiling (abs (/ y-extent nbndpts))))
720               (z-step (ceiling (abs (/ z-extent nbndpts))))
721               )
722         
723          (printf "# x-min: ~A x-max: ~A y-min: ~A y-max ~A z-min: ~A z-max ~A~%" 
724                  x-min x-max y-min y-max z-min z-max)
725         
726         
727          (let ( (bndpts 
728                  (concatenate
729                   (concatenate
730                    (list-tabulate nbndpts 
731                         (lambda (i) 
732                           (list-tabulate nbndpts 
733                              (lambda (j) 
734                                (list-tabulate nbndpts 
735                                   (lambda (k) 
736                                     (make-point (+ x-min (* i x-step)) 
737                                                 (+ y-min (* j y-step)) 
738                                                 (+ z-min (* k z-step)) ))))))
739                         ))
740                   ))
741                  )
742
743
744            (if (and b-radius min-neighbors)
745                (for-each
746                 (lambda (bp)
747                   (let ((pts (kd-tree-near-neighbors point-tree b-radius bp)))
748                     (if (and (not (null? pts))
749                              (<= min-neighbors (length pts)))
750                         (printf "~A ~A ~A~%" (coord 0 bp) (coord 1 bp) (coord 2 bp)))
751                     ))
752                 bndpts)
753                (for-each
754                 (lambda (bp)
755                   (printf "~A ~A ~A~%" (coord 0 bp) (coord 1 bp) (coord 2 bp)))
756                 bndpts)
757                )
758           
759            )))
760      ))
761
762
763
764(define (main)
765
766  (let ((operands          (opt '@)))
767
768    (if (opt 'version) 
769        (begin
770          (printf "neurolucida version ~?" "~A~%" 
771                  (lookup-def 'version (extension-information 'neurolucida)))
772          (exit 0)
773          ))
774
775    (if (null? operands) (neurolucida:usage))
776
777    (d "data directory is ~s~%" (get-data-dir))
778
779    (let* ((format (string->symbol (or (opt 'format) (defopt 'format))))
780           (data-map
781            (concatenate
782             (map (lambda (p)
783                    (let* ((key             (or (opt 'key) (defopt 'key)))
784                           (contents        (parse-sxml p))
785                           (contours        (extract-element contents "contour" key))
786                           (trees           (extract-element contents "tree" key)))
787                      (let recur ((keys (delete-duplicates (map car contours)))
788                                  (elms (append contours trees)) 
789                                  (data-map '()))
790
791                        (if (null? keys)
792                            (if (null? elms) data-map (if (equal? format 'ng) data-map (error "elements with unknown color" elms)))
793                            (let-values (((lst rest)  (partition-by-key elms (car keys) key=?)))
794                              (recur (cdr keys) rest (cons (cons (car keys) (map cadr lst)) data-map)))))
795                        ))
796                  operands)))
797           (ddir (get-data-dir))
798           )
799
800      (if (opt 'forest-grid)
801          (with-output-to-file (make-pathname ddir "forest.grid")
802            (lambda () (print-grid data-map (opt 'forest-grid))))
803
804          (case format
805           
806            ((pts)
807             (for-each (lambda (f x) 
808                         (with-output-to-file f
809                           (lambda () (print-pts x))))
810                       (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) ".pts")) data-map)
811                       data-map))
812           
813            ((swc)
814             (for-each (lambda (f x) 
815                         (with-output-to-file f
816                           (lambda () (print-swc x))))
817                       (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) ".swc")) data-map)
818                       data-map))
819           
820            ((vcg)
821             (for-each (lambda (f x) 
822                         (with-output-to-file f
823                           (lambda () 
824                             ((make-format-graph 'vcg) (current-output-port) (make-tree-graph 'neurolucida x)))))
825                       (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) "vcg")) data-map)
826                       data-map))
827           
828            ((ng)
829             (let* ((permute-coords (opt 'permute-coords))
830                    (permute-coords 
831                     (and permute-coords
832                          (let ((permute-coords (map string->number (string-split permute-coords ","))))
833                            (if (not (and (list? permute-coords)
834                                          (= (length permute-coords) 3)
835                                          (every (lambda (x) (and (integer? x) (<= 1 x) (>= 3 x))) permute-coords)))
836                                (error "invalid point coordinate permutation indices" permute-coords))
837                            (let ((permute-coords (map (lambda (i) (- i 1)) permute-coords)))
838                              (lambda (p) (map (lambda (i) (list-ref p i)) permute-coords))
839                              ))
840                          ))
841                    )
842               
843               (for-each (lambda (f k x) 
844                           (with-output-to-file f
845                             (lambda () 
846                               (pp (make-ng k (make-tree-graph 'neurolucida x) 
847                                            permute-coords: permute-coords)
848                                   (current-output-port)))))
849                         (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) "ng")) data-map)
850                         (map (lambda (x) (car (string-split (cadar x) "#"))) data-map)
851                         data-map)
852               ))
853          ))
854      )))
855
856(main)
857
Note: See TracBrowser for help on using the repository browser.