source: project/release/4/neurolucida/tags/1.5/neurolucida.scm @ 27173

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

neurolucida release 1.5

File size: 11.1 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 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(define (mkdir dir)
44  (system* "mkdir -p ~a" (quotewrap dir)))
45
46
47(define *quiet* #f)
48
49
50(define (d fstr . args)
51  (let ([port (if *quiet* (current-error-port) (current-output-port))])
52    (apply fprintf port fstr args)
53    (flush-output port) ) )
54
55
56(define opt-defaults
57  `(
58    (key             . "color")
59    (format          . "swc")
60    ))
61
62
63(define (defopt x)
64  (lookup-def x opt-defaults))
65
66(define (symbol-upcase str)
67  (string->symbol (string-upcase str)))
68
69
70(define opt-grammar
71  `(
72    (data-dir
73     "set download directory (default is a randomly generated name in /tmp)"
74     (single-char #\d)
75     (value (required DIR)))
76
77    (key
78     "extraction key"
79     (single-char #\k)
80     (value       (required "FIELD")
81                  (default  ,(defopt 'key))))
82
83    (format
84     "output format (swc, asc, vcg)"
85     (single-char #\f)
86     (value       (required "FORMAT")
87                  (default  ,(defopt 'format))))
88
89    (help  "Print help"
90            (single-char #\h))
91 
92  ))
93
94
95;; Use args:usage to generate a formatted list of options (from OPTS),
96;; suitable for embedding into help text.
97(define (neurolucida:usage)
98  (print "Usage: " (car (argv)) " [options...] operands ")
99  (newline)
100  (print "Where operands are XML Neurolucida files")
101  (newline)
102  (print "The following options are recognized: ")
103  (newline)
104  (width 35)
105  (print (parameterize ((indent 5)) (usage opt-grammar)))
106  (exit 1))
107
108
109;; Process arguments and collate options and arguments into OPTIONS
110;; alist, and operands (filenames) into OPERANDS.  You can handle
111;; options as they are processed, or afterwards.
112
113(define opts    (getopt-long (command-line-arguments) opt-grammar))
114(define opt     (make-option-dispatch opts opt-grammar))
115
116
117(define data-dir (make-parameter #f))
118
119(define (get-data-dir)
120  (or (opt 'data-dir)
121      (or (data-dir)
122          (let ([dir (create-temporary-directory)])
123            (data-dir dir)
124            dir ) ) ))
125
126
127(define (create-temporary-directory)
128  (let ((dir (or (get-environment-variable "TMPDIR") 
129                 (get-environment-variable "TEMP") 
130                 (get-environment-variable "TMP") 
131                 "/tmp")))
132    (let loop ()
133      (let* ((n (current-milliseconds)) ; current milliseconds
134             (pn (make-pathname dir (string-append "neurolucida-" (number->string n 16)) "tmp")))
135        (cond ((file-exists? pn) (loop))
136              (else (mkdir pn) pn))))))
137
138
139(define neurolucida-xmlns "http://www.mbfbioscience.com/2007/neurolucida")
140
141
142(define (parse-sxml fpath)
143  (with-input-from-file fpath
144    (lambda () (cons '*TOP* (ssax:xml->sxml (current-input-port) `((nl . ,neurolucida-xmlns)))))))
145
146
147(define (extract-element sxml element-name key-attr)
148  (let ((element-fullname (string->symbol (conc "nl:" (->string element-name)))))
149    (let* ((elements ((sxpath `(// nl:mbf ,element-fullname ))  sxml))
150           (keys     ((sxpath `(// ,element-fullname @ ,key-attr)) `(*TOP* ,@elements))))
151      (zip keys elements))))
152
153
154(define (key=? a b)  (string=? (cadar a) (cadar b)))
155
156
157(define (partition-by-key a k key=?)
158  (partition (lambda (x) (equal? (car x) k)) a))
159
160
161(define (make-tree-graph label x)
162
163  (let* ((sxml       `(*TOP* ,@(cdr x)))
164         (g          (make-digraph label #f))
165         (add-node!  (g 'add-node!))
166         (add-edge!  (g 'add-edge!)))
167
168    (add-node! 0 'soma)
169
170    (let recur ((tree     ((sxpath `(nl:tree))  sxml))
171                (parent   0)
172                (node-id  1))
173
174      (let ((points+spines     ((sxpath `((*or* nl:point nl:spine))) tree))
175            (subtrees          (append ((sxpath `(nl:tree)) tree)
176                                       ((sxpath `(nl:branch)) tree))))
177
178        (let ((n (length points+spines)))
179         
180          ;; insert points and spines in the dependency graph
181          (let ((node-id1 (fold (lambda (n i) (add-node! i n) (+ 1 i)) node-id points+spines)))
182           
183            ;; insert edges
184            (if (positive? n)
185                (let ((fin (- node-id1 1)))
186                  (add-edge! (list parent node-id #f))
187                  (let inner-recur ((i node-id))
188                    (if (< i fin)
189                        (let ((j (+ 1 i))) 
190                          (add-edge! (list i j  #f))
191                          (inner-recur j))
192                        ))
193                  ))
194           
195            (if (pair? subtrees)
196                (let ((branch-node (- node-id1 1))) ;; make the last node to be branching node
197                  ((g 'node-info-set!) branch-node (cons 'branch ((g 'node-info) branch-node)))
198                  (fold (lambda (subtree i) (recur subtree branch-node i)) node-id1 subtrees))
199                node-id1
200                ))
201          ))
202      )
203    g))
204
205
206
207;;
208;; SWC format
209;;
210;; n T x y z R P
211;;
212;; n is an integer label that identifies the current point and
213;; increments by one from one line to the next.
214;;
215;; T is an integer representing the type of neuronal segment, such as
216;; soma, axon, apical dendrite, etc. The standard accepted integer
217;; values are given below.
218;;
219;;     * 0 = undefined
220;;     * 1 = soma
221;;     * 2 = axon
222;;     * 3 = dendrite
223;;     * 4 = apical dendrite
224;;     * 5 = fork point
225;;     * 6 = end point
226;;     * 7 = custom
227;;
228;; x, y, z gives the cartesian coordinates of each node.
229;;
230;; R is the radius at that node.
231;;
232;; P indicates the parent (the integer label) of the current point or
233;; -1 to indicate an origin (soma).
234
235
236(define (print-swc x)
237  (define get (compose car alist-ref))
238
239  (let ((sxml `(*TOP* ,@(cdr x))))
240
241    (let* ((offset (make-parameter 0))
242           (soma ((sxpath `(nl:contour))  sxml))
243           (soma-points ((sxpath `(nl:point @))  soma))
244           (soma-npoints (length soma-points))
245           (soma-indices (list-tabulate soma-npoints (lambda (x) (+ 1 x)))))
246
247       (for-each (lambda (p n)
248                  (let ((data (cdr p)))
249                    (let ((T    1) ;; soma
250                          (x  (get 'x data))
251                          (y  (get 'y data))
252                          (z  (get 'z data))
253                          (R  (number->string (/ (string->number (get 'd data)) 2)))
254                          (P  (- n 1))
255                          )
256                    (printf "~A ~A ~A ~A ~A ~A ~A~%" n T x y z R (if (zero? P) -1 P)))))
257                soma-points soma-indices)
258       
259       (offset soma-npoints)
260
261       (let recur ((tree ((sxpath `(nl:tree))  sxml)) 
262                   (parent-index soma-npoints))
263
264         (let ((points     ((sxpath `(nl:point @))  tree))
265               (subtrees   (append ((sxpath `(nl:tree))  tree)
266                                   ((sxpath `(nl:branch))  tree))))
267           
268           (let* ((npoints (length points))
269                  (indices (let ((k (+ 1 (offset)))) (list-tabulate npoints (lambda (x) (+ k x)))))
270                  (parents (cons parent-index indices)))
271             
272             (fold (lambda (pt i par n)
273                         (let ((data (cdr pt)))
274                           (let ((T    (if (fx> n 0) 3 (if (null? subtrees) 6 5))) ;; dendrite, or end point/bifurcation point
275                                 (x  (get 'x data))
276                                 (y  (get 'y data))
277                                 (z  (get 'z data))
278                                 (R  (number->string (/ (string->number (get 'd data)) 2)))
279                                 (P  par)
280                                 )
281                             (printf "~A ~A ~A ~A ~A ~A ~A~%" 
282                                     i T x y z R P)
283                             (fx- n 1)
284                             )))
285                   (fx- npoints 1)
286                   points indices parents)
287             
288             (let ((parent-index1 (+ npoints (offset))))
289
290               (offset (+ (offset) npoints)) 
291
292               (if (pair? subtrees)
293                   (for-each (lambda (subtree) (recur subtree parent-index1) )
294                             subtrees)))
295             )))
296
297    )))
298
299
300(define (print-pts x)
301  (define get (compose car alist-ref))
302
303  (let ((sxml `(*TOP* ,@(cdr x))))
304
305    (let* ((soma ((sxpath `(nl:contour))  sxml))
306           (soma-points ((sxpath `(nl:point @))  soma))
307           (soma-npoints (length soma-points)))
308
309       (for-each (lambda (p)
310                  (let ((data (cdr p)))
311                    (let ((x  (get 'x data))
312                          (y  (get 'y data))
313                          (z  (get 'z data))
314                          (R  (/ (string->number (get 'd data)) 2))
315                          )
316                      (printf "~A ~A ~A~%" x y z))))
317                soma-points)
318
319       (let recur ((tree ((sxpath `(nl:tree))  sxml))  )
320
321         (let ((points     ((sxpath `(nl:point @))  tree))
322               (subtrees   (append ((sxpath `(nl:tree))  tree)
323                                   ((sxpath `(nl:branch))  tree))))
324           
325           (let ((npoints (length points)))
326             
327             (for-each (lambda (pt)
328                         (let ((data (cdr pt)))
329                           (let ((x  (string->number (get 'x data)))
330                                 (y  (string->number (get 'y data)))
331                                 (z  (string->number (get 'z data)))
332                                 (r  (/ (string->number (get 'd data)) 2))
333                                 )
334                             (printf "~A ~A ~A~%" x y z)
335                             )))
336                       points)
337             
338               (if (pair? subtrees)
339                   (for-each (lambda (subtree) (recur subtree) )
340                             subtrees))
341             )))
342
343    )))
344
345
346
347
348
349
350
351(define (main)
352
353  (let ((operands          (opt '@)))
354
355    (if (null? operands) (neurolucida:usage))
356
357    (d "data directory is ~s~%" (get-data-dir))
358
359    (let* ((data-map
360            (concatenate
361             (map (lambda (p)
362                    (let* ((key             (or (opt 'key) (defopt 'key)))
363                           (contents        (parse-sxml p))
364                           (contours        (extract-element contents "contour" key))
365                           (trees           (extract-element contents "tree" key)))
366                      (let recur ((keys (delete-duplicates (map car contours)))
367                                  (elms (append contours trees)) 
368                                  (data-map '()))
369                        (if (null? keys)
370                            (if (null? elms) data-map (error "elements with unknown color" elms))
371                            (let-values (((lst rest)  (partition-by-key elms (car keys) key=?)))
372                              (recur (cdr keys) rest (cons (cons (car keys) (map cadr lst)) data-map)))))
373                        ))
374                  operands)))
375           (ddir (get-data-dir))
376           )
377
378      (let ((format (string->symbol (or (opt 'format) (defopt 'format)))))
379
380        (case format
381          ((pts)
382           (for-each (lambda (f x) 
383                       (with-output-to-file f
384                         (lambda () (print-pts x))))
385                     (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) ".asc")) data-map)
386                     data-map))
387
388          ((swc)
389           (for-each (lambda (f x) 
390                       (with-output-to-file f
391                         (lambda () (print-swc x))))
392                     (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) ".swc")) data-map)
393                     data-map))
394
395          ((vcg)
396           (for-each (lambda (f x) 
397                       (with-output-to-file f
398                         (lambda () 
399                           ((make-format-graph 'vcg) (current-output-port) (make-tree-graph 'neurolucida x)))))
400                     (map (lambda (x) (make-pathname ddir (car (string-split (cadar x) "#")) "vcg")) data-map)
401                     data-map))
402          ))
403      )))
404
405(main)
406
Note: See TracBrowser for help on using the repository browser.