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

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

neurolucida: fix in swc output to print correct types for forks and end points

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