source: project/release/4/neurolucida/tags/1.7/neurolucida.scm @ 27325

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

neurolucida version 1.7

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