DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkSphereTessellator.hxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Sphere Tessellator
4 
5  Copyright (c) Pew-Thian Yap. All rights reserved.
6  See http://www.unc.edu/~ptyap/ for details.
7 
8  This software is distributed WITHOUT ANY WARRANTY; without even
9  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
10  PURPOSE. See the above copyright notices for more information.
11 
12  =========================================================================*/
13 
14 #ifndef _itkSphereTessellator_hxx
15 #define _itkSphereTessellator_hxx
16 
17 #include <cstdio>
18 #include <cstdlib>
19 #include <cstring>
20 #include <cmath>
21 
22 #include "itkSphereTessellator.h"
23 
24 namespace itk
25 {
26 
27 template <typename TElement>
30 {
31  m_BasicShape = ICOSAHEDRON;
32  m_Order = 1;
33 
34  vertices = NULL;
35  faces = NULL;
36  start = NULL;
37  end = NULL;
38  midpoint = NULL;
39 
40  m_NumberOfVertices = 0;
41  m_NumberOfFaces = 0;
42  m_NumberOfEdges = 0;
43 
44  m_Points.clear();
45  m_Cells.clear();
46 }
47 
48 template <typename TElement>
51 {
52  if (vertices) free (vertices);
53  if (faces) free (faces);
54 }
55 
60 template <typename TElement>
61 void
63 ::Reserve(void)
64 {
65  // Do nothing; vnl_matrix cannot reserve elements
66 }
67 
68 
73 template <typename TElement>
74 void
76 ::Squeeze(void)
77 {
78  // Do nothing; vnl_matrix cannot be squeezed
79 }
80 
81 
85 template <typename TElement>
86 void
89 {
90  m_Points.clear();
91  m_Cells.clear();
92 
93  this->Modified();
94 }
95 
96 template <typename TElement>
97 void
100 {
101  double sqrt3 = 1 / vcl_sqrt(3.0);
102  double tetrahedron_vertices[] = {sqrt3, sqrt3, sqrt3,
103  -sqrt3, -sqrt3, sqrt3,
104  -sqrt3, sqrt3, -sqrt3,
105  sqrt3, -sqrt3, -sqrt3};
106  int tetrahedron_faces[] = {0, 2, 1, 0, 1, 3, 2, 3, 1, 3, 2, 0};
107 
108  n_vertices = 4;
109  n_faces = 4;
110  n_edges = 6;
111  vertices = (double*)malloc(3*n_vertices*sizeof(double));
112  faces = (int*)malloc(3*n_faces*sizeof(int));
113  memcpy ((void*)vertices, (void*)tetrahedron_vertices, 3*n_vertices*sizeof(double));
114  memcpy ((void*)faces, (void*)tetrahedron_faces, 3*n_faces*sizeof(int));
115 }
116 
117 template <typename TElement>
118 void
121 {
122  double octahedron_vertices[] = {0.0, 0.0, -1.0,
123  1.0, 0.0, 0.0,
124  0.0, -1.0, 0.0,
125  -1.0, 0.0, 0.0,
126  0.0, 1.0, 0.0,
127  0.0, 0.0, 1.0};
128  int octahedron_faces[] = {0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1, 5, 2, 1, 5, 3, 2, 5, 4, 3, 5, 1, 4};
129 
130  n_vertices = 6;
131  n_faces = 8;
132  n_edges = 12;
133  vertices = (double*)malloc(3*n_vertices*sizeof(double));
134  faces = (int*)malloc(3*n_faces*sizeof(int));
135  memcpy ((void*)vertices, (void*)octahedron_vertices, 3*n_vertices*sizeof(double));
136  memcpy ((void*)faces, (void*)octahedron_faces, 3*n_faces*sizeof(int));
137 }
138 
139 template <typename TElement>
140 void
143 {
144  double t = (1+vcl_sqrt(5))/2;
145  double tau = t/vcl_sqrt(1+t*t);
146  double one = 1/vcl_sqrt(1+t*t);
147 
148  double icosahedron_vertices[] = {tau, one, 0.0,
149  -tau, one, 0.0,
150  -tau, -one, 0.0,
151  tau, -one, 0.0,
152  one, 0.0 , tau,
153  one, 0.0 , -tau,
154  -one, 0.0 , -tau,
155  -one, 0.0 , tau,
156  0.0 , tau, one,
157  0.0 , -tau, one,
158  0.0 , -tau, -one,
159  0.0 , tau, -one};
160  int icosahedron_faces[] = {4, 8, 7,
161  4, 7, 9,
162  5, 6, 11,
163  5, 10, 6,
164  0, 4, 3,
165  0, 3, 5,
166  2, 7, 1,
167  2, 1, 6,
168  8, 0, 11,
169  8, 11, 1,
170  9, 10, 3,
171  9, 2, 10,
172  8, 4, 0,
173  11, 0, 5,
174  4, 9, 3,
175  5, 3, 10,
176  7, 8, 1,
177  6, 1, 11,
178  7, 2, 9,
179  6, 10, 2};
180 
181  n_vertices = 12;
182  n_faces = 20;
183  n_edges = 30;
184  vertices = (double*)malloc(3*n_vertices*sizeof(double));
185  faces = (int*)malloc(3*n_faces*sizeof(int));
186  memcpy ((void*)vertices, (void*)icosahedron_vertices, 3*n_vertices*sizeof(double));
187  memcpy ((void*)faces, (void*)icosahedron_faces, 3*n_faces*sizeof(int));
188 }
189 
190 template <typename TElement>
191 int
193 ::SearchMidpoint (int index_start, int index_end)
194 {
195  unsigned int i;
196  for (i=0; i<edge_walk; i++)
197  if ((start[i] == index_start && end[i] == index_end) ||
198  (start[i] == index_end && end[i] == index_start))
199  {
200  int res = midpoint[i];
201 
202  /* update the arrays */
203  start[i] = start[edge_walk-1];
204  end[i] = end[edge_walk-1];
205  midpoint[i] = midpoint[edge_walk-1];
206  edge_walk--;
207 
208  return res;
209  }
210 
211  /* vertex not in the list, so we add it */
212  start[edge_walk] = index_start;
213  end[edge_walk] = index_end;
214  midpoint[edge_walk] = n_vertices;
215 
216  /* create new vertex */
217  vertices[3*n_vertices] = (vertices[3*index_start] + vertices[3*index_end]) / 2.0;
218  vertices[3*n_vertices+1] = (vertices[3*index_start+1] + vertices[3*index_end+1]) / 2.0;
219  vertices[3*n_vertices+2] = (vertices[3*index_start+2] + vertices[3*index_end+2]) / 2.0;
220 
221  /* normalize the new vertex */
222  double length = vcl_sqrt (vertices[3*n_vertices] * vertices[3*n_vertices] +
223  vertices[3*n_vertices+1] * vertices[3*n_vertices+1] +
224  vertices[3*n_vertices+2] * vertices[3*n_vertices+2]);
225  length = 1/length;
226  vertices[3*n_vertices] *= length;
227  vertices[3*n_vertices+1] *= length;
228  vertices[3*n_vertices+2] *= length;
229 
230  n_vertices++;
231  edge_walk++;
232  return midpoint[edge_walk-1];
233 }
234 
235 template <typename TElement>
236 void
238 ::Subdivide (void)
239 {
240  unsigned int n_vertices_new = n_vertices+2*n_edges;
241  unsigned int n_faces_new = 4*n_faces;
242  unsigned int i;
243 
244  edge_walk = 0;
245  n_edges = 2*n_vertices + 3*n_faces;
246  start = (int*)malloc(n_edges*sizeof (int));
247  end = (int*)malloc(n_edges*sizeof (int));
248  midpoint = (int*)malloc(n_edges*sizeof (int));
249 
250  int *faces_old = (int*)malloc (3*n_faces*sizeof(int));
251  faces_old = (int*)memcpy((void*)faces_old, (void*)faces, 3*n_faces*sizeof(int));
252  vertices = (double*)realloc ((void*)vertices, 3*n_vertices_new*sizeof(double));
253  faces = (int*)realloc ((void*)faces, 3*n_faces_new*sizeof(int));
254  n_faces_new = 0;
255 
256  for (i=0; i<n_faces; i++)
257  {
258  int a = faces_old[3*i];
259  int b = faces_old[3*i+1];
260  int c = faces_old[3*i+2];
261 
262  int ab_midpoint = SearchMidpoint (b, a);
263  int bc_midpoint = SearchMidpoint (c, b);
264  int ca_midpoint = SearchMidpoint (a, c);
265 
266  faces[3*n_faces_new] = a;
267  faces[3*n_faces_new+1] = ab_midpoint;
268  faces[3*n_faces_new+2] = ca_midpoint;
269  n_faces_new++;
270  faces[3*n_faces_new] = ca_midpoint;
271  faces[3*n_faces_new+1] = ab_midpoint;
272  faces[3*n_faces_new+2] = bc_midpoint;
273  n_faces_new++;
274  faces[3*n_faces_new] = ca_midpoint;
275  faces[3*n_faces_new+1] = bc_midpoint;
276  faces[3*n_faces_new+2] = c;
277  n_faces_new++;
278  faces[3*n_faces_new] = ab_midpoint;
279  faces[3*n_faces_new+1] = b;
280  faces[3*n_faces_new+2] = bc_midpoint;
281  n_faces_new++;
282  }
283  n_faces = n_faces_new;
284  free (start);
285  free (end);
286  free (midpoint);
287  free (faces_old);
288 }
289 
293 template <typename TElement>
294 void
297 {
298  switch (m_BasicShape)
299  {
300  case TETRAHEDRON:
301  InitTetrahedron ();
302  break;
303  case OCTAHEDRON:
304  InitOctahedron ();
305  break;
306  case ICOSAHEDRON:
307  InitIcosahedron ();
308  break;
309  }
310 
311  for (unsigned int k=1; k < m_Order; k++)
312  Subdivide ();
313 
314  OutputSphere ();
315 }
316 
317 template <typename TElement>
318 void
321 {
322 // std::cout fprintf (ptr, "OFF\n%d %d %d\n", n_vertices, n_faces, n_edges);
323  m_NumberOfVertices = n_vertices;
324  m_NumberOfFaces = n_faces;
325  m_NumberOfEdges = n_edges;
326 
327  m_Points.set_size(n_vertices, 3);
328  m_Cells.set_size(n_faces, 3);
329 
330  for (unsigned int k=0; k<n_vertices; k++)
331  {
332  m_Points(k,0) = vertices[3*k];
333  m_Points(k,1) = vertices[3*k+1];
334  m_Points(k,2) = vertices[3*k+2];
335  }
336  for (unsigned int k=0; k<n_faces; k++)
337  {
338  m_Cells(k,0) = faces[3*k];
339  m_Cells(k,1) = faces[3*k+1];
340  m_Cells(k,2) = faces[3*k+2];
341  }
342 }
343 
344 template <typename TElement>
345 void
347 ::PrintSelf(std::ostream& os, Indent indent) const
348 {
349  Superclass::PrintSelf(os,indent);
350  os << indent << "tessellation order = " << m_Order << std::endl;
351 }
352 
353 } // end namespace itk
354 
355 #endif
void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
int SearchMidpoint(int index_start, int index_end)