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

Last change on this file since 27232 was 27232, checked in by Ivan Raikov, 9 years ago

neurolucida: corrected command line help string

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