Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <string>
34 : #include <fstream>
35 : #include <sstream>
36 : #include <vector>
37 : #include <queue>
38 : #include <map>
39 : #include "file_io_swc.h"
40 : #include "common/error.h"
41 : #include "common/util/string_util.h"
42 : #include "lib_grid/global_attachments.h"
43 : #include "lib_grid/algorithms/subset_util.h" // EraseEmptySubsets
44 : #include "lib_grid/algorithms/element_side_util.h" // GetOpposingSide
45 : #include <boost/lexical_cast.hpp>
46 :
47 : using namespace std;
48 :
49 : namespace ug{
50 :
51 :
52 0 : bool FileReaderSWC::load_file(const char* fileName)
53 : {
54 : m_vPts.clear();
55 :
56 0 : std::ifstream inFile(fileName);
57 0 : if (!inFile)
58 : {
59 0 : UG_LOGN("SWC input file '" << fileName << "' could not be opened for reading.");
60 : return false;
61 : }
62 :
63 : // read line by line
64 : std::string line;
65 : size_t lineCnt = 0;
66 0 : size_t curInd = 0;
67 : std::map<int, size_t> indexMap;
68 0 : while (std::getline(inFile, line))
69 : {
70 0 : ++lineCnt;
71 :
72 : // trim whitespace
73 0 : line = TrimString(line);
74 :
75 : // ignore anything from possible '#' onwards
76 : size_t nChar = line.size();
77 0 : for (size_t i = 0; i < nChar; ++i)
78 : {
79 0 : if (line.at(i) == '#')
80 : {
81 0 : line = line.substr(0, i);
82 0 : break;
83 : }
84 : }
85 :
86 : // empty lines can be ignored
87 0 : if (line.empty()) continue;
88 :
89 : // split the line into tokens
90 0 : std::istringstream buf(line);
91 0 : std::istream_iterator<std::string> beg(buf), end;
92 0 : std::vector<std::string> strs(beg, end);
93 :
94 : // assert number of tokens is correct
95 0 : if (strs.size() != 7)
96 : {
97 0 : UG_LOGN("Error reading SWC file '" << fileName << "': Line " << lineCnt
98 : << " does not contain exactly 7 values.");
99 : return false;
100 : }
101 :
102 :
103 : // collect data for SWC point
104 0 : m_vPts.resize(m_vPts.size() + 1);
105 : SWCPoint& pt = m_vPts.back();
106 :
107 : // get index from file and map to our index
108 0 : indexMap[boost::lexical_cast<int>(strs[0])] = curInd;
109 :
110 : // type
111 0 : int type = boost::lexical_cast<int>(strs[1]);
112 0 : switch (type)
113 : {
114 0 : case 0: pt.type = swc_types::SWC_UNDF; break;
115 0 : case 1: pt.type = swc_types::SWC_SOMA; break;
116 0 : case 2: pt.type = swc_types::SWC_AXON; break;
117 0 : case 3: pt.type = swc_types::SWC_DEND; break;
118 0 : case 4: pt.type = swc_types::SWC_APIC; break;
119 0 : case 5: pt.type = swc_types::SWC_FORK; break;
120 0 : case 6: pt.type = swc_types::SWC_END; break;
121 0 : default: pt.type = swc_types::SWC_CUSTOM;
122 : }
123 :
124 : // coordinates
125 0 : pt.coords.x() = boost::lexical_cast<number>(strs[2]);
126 0 : pt.coords.y() = boost::lexical_cast<number>(strs[3]);
127 0 : pt.coords.z() = boost::lexical_cast<number>(strs[4]);
128 :
129 : // radius
130 0 : pt.radius = boost::lexical_cast<number>(strs[5]);
131 :
132 : // connections
133 0 : int conn = boost::lexical_cast<int>(strs[6]);
134 0 : if (conn >= 0)
135 : {
136 : std::map<int, size_t>::const_iterator it = indexMap.find(conn);
137 0 : if (it == indexMap.end())
138 : {
139 0 : UG_LOGN("Error reading SWC file '" << fileName << "': Line " << lineCnt
140 : << " refers to unknown parent index " << conn << ".");
141 0 : return false;
142 : }
143 :
144 0 : size_t parentID = indexMap[conn];
145 0 : pt.conns.push_back(parentID);
146 0 : m_vPts[parentID].conns.push_back(curInd);
147 : }
148 :
149 : // increase current point index
150 0 : ++curInd;
151 0 : }
152 :
153 : return true;
154 0 : }
155 :
156 :
157 :
158 0 : bool FileReaderSWC::create_grid(Grid& g, ISubsetHandler* pSH, number scale_length)
159 : {
160 0 : if (!g.has_vertex_attachment(aPosition))
161 : g.attach_to_vertices(aPosition);
162 : Grid::VertexAttachmentAccessor<APosition> aaPos(g, aPosition);
163 :
164 : // handle diameter attachment
165 0 : if (!GlobalAttachments::is_declared("diameter"))
166 : {
167 0 : try {GlobalAttachments::declare_attachment<ANumber>("diameter");}
168 0 : catch (ug::UGError& err)
169 : {
170 : UG_LOGN(err.get_msg());
171 : return false;
172 0 : }
173 : }
174 :
175 0 : ANumber aDiam = GlobalAttachments::attachment<ANumber>("diameter");
176 0 : if (!g.has_vertex_attachment(aDiam))
177 : g.attach_to_vertices(aDiam, true);
178 :
179 : Grid::AttachmentAccessor<Vertex, ANumber> aaDiam(g, aDiam);
180 :
181 : // create grid
182 : const size_t nP = m_vPts.size();
183 0 : std::vector<Vertex*> vrts(nP, NULL);
184 0 : for (size_t i = 0; i < nP; ++i)
185 : {
186 : const SWCPoint& pt = m_vPts[i];
187 :
188 : // create vertex and save
189 0 : Vertex* v = vrts[i] = *g.create<RegularVertex>();
190 : VecScale(aaPos[v], pt.coords, scale_length);
191 0 : pSH->assign_subset(v, pt.type - 1);
192 0 : aaDiam[v] = 2 * pt.radius * scale_length;
193 :
194 : // create edge connections to already created vertices
195 0 : for (size_t j = 0; j < pt.conns.size(); ++j)
196 : {
197 0 : if (pt.conns[j] < i)
198 : {
199 0 : Edge* e = *g.create<RegularEdge>(EdgeDescriptor(vrts[pt.conns[j]], v));
200 0 : pSH->assign_subset(e, pt.type - 1);
201 : }
202 : }
203 : }
204 :
205 : // final subset management
206 0 : AssignSubsetColors(*pSH);
207 0 : pSH->set_subset_name("soma", 0);
208 0 : pSH->set_subset_name("axon", 1);
209 0 : pSH->set_subset_name("dend", 2);
210 0 : pSH->set_subset_name("apic", 3);
211 0 : pSH->set_subset_name("fork", 4);
212 0 : pSH->set_subset_name("end", 5);
213 0 : pSH->set_subset_name("custom", 6);
214 0 : EraseEmptySubsets(*pSH);
215 :
216 : return true;
217 0 : }
218 :
219 :
220 :
221 0 : const std::vector<swc_types::SWCPoint>& FileReaderSWC::swc_points() const
222 : {
223 0 : return m_vPts;
224 : }
225 :
226 :
227 :
228 0 : std::vector<swc_types::SWCPoint>& FileReaderSWC::swc_points()
229 : {
230 0 : return m_vPts;
231 : }
232 :
233 :
234 :
235 0 : bool LoadGridFromSWC(Grid& g, ISubsetHandler* pSH, const char* fileName, AVector3& aPos)
236 : {
237 : // now read the file
238 : FileReaderSWC reader;
239 0 : return reader.load_file(fileName) && reader.create_grid(g, pSH, 1.0);
240 : }
241 :
242 :
243 :
244 0 : bool ExportGridToSWC(Grid& g, ISubsetHandler* pSH, const char* fileName, AVector3& aPos)
245 : {
246 : // get access to positions
247 0 : if (!g.has_vertex_attachment(aPosition))
248 : {
249 : UG_LOGN("Position attachment not attached to grid.");
250 0 : return false;
251 : }
252 : Grid::VertexAttachmentAccessor<APosition> aaPos(g, aPosition);
253 :
254 : // get access to diameter attachment
255 0 : if (!GlobalAttachments::is_declared("diameter"))
256 : {
257 0 : try {GlobalAttachments::declare_attachment<ANumber>("diameter");}
258 0 : catch (ug::UGError& err)
259 : {
260 : UG_LOGN(err.get_msg());
261 : return false;
262 0 : }
263 : }
264 0 : ANumber aDiam = GlobalAttachments::attachment<ANumber>("diameter");
265 0 : if (!g.has_vertex_attachment(aDiam))
266 : {
267 : UG_LOGN("WARNING: No diameter attachment attached to grid. "
268 : "Will use 1.0 as default diameter.");
269 0 : g.attach_to_vertices_dv(aDiam, 1.0);
270 : }
271 : Grid::AttachmentAccessor<Vertex, ANumber> aaDiam(g, aDiam);
272 :
273 : // analyze subset names to find out corresponding swc-types
274 0 : size_t nss = pSH->num_subsets();
275 0 : std::vector<size_t> vType(nss);
276 : bool soma_subset_present = false;
277 0 : for (size_t i = 0; i < nss; ++i)
278 : {
279 0 : std::string name(pSH->get_subset_name(i));
280 : std::transform(name.begin(), name.end(), name.begin(), ::toupper);
281 0 : if (name.find("SOMA") != std::string::npos)
282 : {
283 : soma_subset_present = true;
284 0 : vType[i] = 1;
285 : }
286 0 : else if (name.find("AXON") != std::string::npos)
287 0 : vType[i] = 2;
288 0 : else if (name.find("APIC") != std::string::npos)
289 0 : vType[i] = 4;
290 0 : else if (name.find("DEND") != std::string::npos)
291 0 : vType[i] = 3;
292 0 : else if (name.find("END") != std::string::npos)
293 0 : vType[i] = 6;
294 0 : else if (name.find("FORK") != std::string::npos)
295 0 : vType[i] = 5;
296 0 : else if (name.find("CUSTOM") != std::string::npos)
297 0 : vType[i] = 7;
298 0 : else vType[i] = 0;
299 : }
300 :
301 0 : if (!soma_subset_present)
302 : UG_LOGN("Warning: No somatic subset could be identified.")
303 :
304 0 : if (g.begin<Vertex>() == g.end<Vertex>())
305 : {
306 : UG_LOGN("Warning: No vertices contained in grid.")
307 : return true;
308 : }
309 :
310 : // find soma vertex (if identifiable)
311 0 : Vertex* start = *g.begin<Vertex>();
312 0 : if (soma_subset_present)
313 : {
314 0 : g.begin_marking();
315 : std::queue<Vertex*> q; // corresponds to breadth-first
316 : q.push(start);
317 0 : while (!q.empty())
318 : {
319 0 : Vertex* v = q.front();
320 0 : if (vType[pSH->get_subset_index(v)] == 1) break;
321 : g.mark(v);
322 : q.pop();
323 :
324 : // push neighboring elems to queue
325 : Grid::traits<Edge>::secure_container edges;
326 : g.associated_elements(edges, v);
327 :
328 : size_t sz = edges.size();
329 0 : for (size_t e = 0; e < sz; ++e)
330 : {
331 0 : Vertex* otherEnd = GetOpposingSide(g, edges[e], v);
332 0 : if (!g.is_marked(otherEnd))
333 : q.push(otherEnd);
334 : }
335 : }
336 0 : g.end_marking();
337 :
338 0 : if (q.empty())
339 : UG_LOGN("Warning: No soma vertex could be found in the requested neuron.")
340 : else
341 0 : start = q.front();
342 : }
343 :
344 : // write the neuron to file
345 0 : std::ofstream outFile(fileName, std::ios::out);
346 0 : UG_COND_THROW(!outFile.is_open(), "Could not open output file '" << fileName << "'.");
347 :
348 : outFile << "# This file has been generated by UG4." << std::endl;
349 :
350 : std::stack<std::pair<Vertex*, int> > stack; // corresponds to depth-first
351 0 : stack.push(std::make_pair(start, -1));
352 :
353 0 : g.begin_marking();
354 : int ind = 0; // by convention, swc starts with index 1
355 : bool all_types_identified = true;
356 0 : while (!stack.empty())
357 : {
358 : // get all infos regarding vertex
359 : std::pair<Vertex*, int>& info = stack.top();
360 0 : Vertex* v = info.first;
361 0 : int conn = info.second;
362 : stack.pop();
363 :
364 : // mark curr vrt
365 : g.mark(v);
366 :
367 0 : size_t type = vType[pSH->get_subset_index(v)];
368 0 : if (!type) all_types_identified = false;
369 :
370 : const vector3& coord = aaPos[v];
371 :
372 0 : number radius = 0.5*aaDiam[v];
373 :
374 : // write line to file
375 0 : outFile << ++ind << " " << type << " "
376 0 : << coord[0] << " " << coord[1] << " " << coord[2] << " "
377 0 : << radius << " " << conn << std::endl;
378 :
379 : // push neighboring elems to queue
380 : Grid::traits<Edge>::secure_container edges;
381 : g.associated_elements(edges, v);
382 :
383 : size_t sz = edges.size();
384 0 : for (size_t e = 0; e < sz; ++e)
385 : {
386 0 : Vertex* otherEnd = GetOpposingSide(g, edges[e], v);
387 0 : if (!g.is_marked(otherEnd))
388 0 : stack.push(std::make_pair(otherEnd, ind));
389 : }
390 : }
391 0 : g.end_marking();
392 :
393 0 : if (!all_types_identified)
394 : UG_LOGN("WARNING: Some vertex type(s) - soma, dendrite, axon, etc. -\n"
395 : "could not be identified using the subset names.\n"
396 : << "To ensure correct types in the resulting swc file, the ugx subset names\n"
397 : "need to contain one of the strings \"SOMA\", \"AXON\", \"DEND\", \"APIC\",\n"
398 : "upper/lower case can be ignored.");
399 :
400 0 : outFile.close();
401 :
402 : return true;
403 0 : }
404 :
405 :
406 :
407 : }// end of namespace
|
