A.I.M.S algorithms


meshToVoxelsResampler_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
35 #include <aims/distancemap/stlsort.h> //for Point3dfCompare
37 #include <aims/resampling/motion.h>
38 #include <map>
39 #include <float.h>
40 
41 namespace aims
42 {
43 
44 template<typename O> void MeshToVoxelsResampler<O>::fill_header(
45 PythonHeader &hdr, O &output, const AimsVector<float,3> & offset, float spacing)
46 const
47 {
48  std::vector<float> voxel_size(3, spacing);
49  std::vector<std::string> refs(1, "to_mesh");
50  std::vector<std::vector<float> >trans;
51  Motion mot;
52  mot.setToIdentity();
53  mot.translation() = offset;
54  trans.push_back( mot.toVector() );
55 
56  hdr.setProperty("voxel_size", voxel_size);
57  hdr.setProperty("referentials", refs);
58  hdr.setProperty("transformations", trans);
59  output.setSizeXYZT(spacing, spacing, spacing, 1.);
60 }
61 
67 template<typename O>
69  unsigned int connexity)
70 {
71  const std::vector<AimsVector<float, 3> > &vertices = m.vertex();
72  const std::vector<AimsVector<uint, 3> > &polygons = m.polygon();
73  float xmin, ymin, zmin, xmax, ymax, zmax;
74  unsigned int vsize = vertices.size();
75  unsigned int psize = polygons.size();
76  float r2;
77  float tc;
78 
79  if (connexity == 26)
80  {
81  r2 = 0.75;
82  tc = 0.5;
83  }
84  else if (connexity == 6)
85  {
86  r2 = 0.25;
87  tc = sqrt(2.) / 4.;
88  }
89 
90  xmin = ymin = zmin = FLT_MAX;
91  xmax = ymax = zmax = -FLT_MAX;
92 
93  // bounding box
94  for (unsigned int i = 0; i < vsize; ++i)
95  {
96  const AimsVector<float, 3> &v = vertices[i];
97  if (v[0] < xmin) xmin = v[0];
98  if (v[0] > xmax) xmax = v[0];
99  if (v[1] < ymin) ymin = v[1];
100  if (v[1] > ymax) ymax = v[1];
101  if (v[2] < zmin) zmin = v[2];
102  if (v[2] > zmax) zmax = v[2];
103  }
104  double invspacing = 1./spacing;
105  // borders min/max are not optimals
106  xmin -= 1;
107  ymin -= 1;
108  zmin -= 1;
109  AimsVector<float,3> offset(xmin, ymin, zmin);
110  unsigned int dimx = (unsigned int) round((xmax - xmin) * invspacing) +2;
111  unsigned int dimy = (unsigned int) round((ymax - ymin) * invspacing) +2;
112  unsigned int dimz = (unsigned int) round((zmax - zmin) * invspacing) +2;
113 
114  O output = this->init_data(offset, spacing, dimx, dimy, dimz);
115 
117 
118  for (unsigned int ind = 1; ind < psize + 1; ++ind)
119  {
120  const AimsVector<uint,3> &p = polygons[ind - 1];
121 
122  t[0] = (vertices[p[0]] - offset) * invspacing;
123  t[1] = (vertices[p[1]] - offset) * invspacing;
124  t[2] = (vertices[p[2]] - offset) * invspacing;
125 
127  // vertices (sphere with a radius of sqrt(r2))
128  for (unsigned int j = 0; j < 3; ++j)
129  {
130  t1[0] = ceil(t[j][0] - 0.5);
131  t1[1] = ceil(t[j][1] - 0.5);
132  t1[2] = ceil(t[j][2] - 0.5);
133  for (float dx = -1; dx <= 1; dx += 1)
134  for (float dy = -1; dy <= 1; dy += 1)
135  for (float dz = -1; dz <= 1; dz += 1)
136  {
137  uint x2 = (uint) ceil(t1[0] + dx);
138  uint y2 = (uint) ceil(t1[1] + dy);
139  uint z2 = (uint) ceil(t1[2] + dz);
140  float dix = x2 - t1[0];
141  float diy = y2 - t1[1];
142  float diz = z2 - t1[2];
143  if (dix * dix + diy * diy + diz * diz <= r2)
144  this->set(output, x2, y2, z2, ind);
145  }
146  }
147 
148  // edges (cylinder with a radius of sqrt(r2))
149  float diff[3][3];
150  uint usize[3];
151  AimsVector<float, 3> t3, t4, t5;
152  for (unsigned int j = 0; j < 3; ++j)
153  {
154  const AimsVector<float, 3> &t2 = t[(j + 1) % 3];
155  float dix = (t[j][0] - t2[0] - 0.5);
156  float diy = (t[j][1] - t2[1] - 0.5);
157  float diz = (t[j][2] - t2[2] - 0.5);
158  diff[j][0] = dix;
159  diff[j][1] = diy;
160  diff[j][2] = diz;
161  float size2 = dix * dix + diy * diy + diz * diz;
162  if (!size2) break;
163  //usizej : not optimal value (for speed considerations)
164  uint usizej = (uint) round(sqrt(size2) * 2.);
165  usize[j] = usizej;
166  std::map<Point3df, bool, Point3dfCompare> map;
167  for (uint l = 0; l <= usizej; ++l)
168  {
169  float lambda = ((float) l) / usizej;
170  t3 = t[j] * lambda + t2 * (1 - lambda);
171  t3[0] = ceil(t3[0]);
172  t3[1] = ceil(t3[1]);
173  t3[2] = ceil(t3[2]);
174  for (float dx = -1; dx <= 1; dx += 1)
175  for (float dy = -1; dy <= 1; dy += 1)
176  for (float dz = -1; dz <= 1; dz += 1)
177  {
178  Point3df p;
179  p[0] = t3[0] + dx;
180  p[1] = t3[1] + dy;
181  p[2] = t3[2] + dz;
182  map[p] = true;
183  }
184  }
185  std::map<Point3df,bool>::const_iterator i, e;
186  for (i = map.begin(), e = map.end(); i != e; ++i)
187  {
188  const Point3df &p0 = (*i).first;
189  Point3df p = t[j] - p0;
190  float vx = diy * p[2] - diz * p[1];
191  float vy = diz * p[0] - dix * p[2];
192  float vz = dix * p[1] - diy * p[0];
193  if ((vx * vx + vy * vy + vz * vz) / size2 <= r2)
194  this->set(output, (uint) p0[0],
195  (uint) p0[1],
196  (uint) p0[2], ind);
197  }
198  }
199 
200  // fill triangles (between 2 planes at a distance of dlnk)
201  std::map<Point3df, bool, Point3dfCompare> map;
202  for (uint l = 0; l <= usize[0]; ++l)
203  {
204  float lambda = ((float) l) / usize[0];
205  float msize = usize[1] * lambda;
206  t3 = t[0] * (1 - lambda) + t[1] * lambda;
207  t4 = t[0] * (1 - lambda) + t[2] * lambda;
208  for (uint m = 0; m < msize; ++m)
209  {
210  float mu = ((float) m) / msize;
211  t5 = t3 * (1 - mu) + t4 * mu;
212  t5[0] = ceil(t5[0] - 0.5);
213  t5[1] = ceil(t5[1] - 0.5);
214  t5[2] = ceil(t5[2] - 0.5);
215  for (float dx = -1; dx <= 1; dx += 1)
216  for (float dy = -1; dy <= 1; dy += 1)
217  for (float dz = -1; dz <= 1; dz += 1)
218  {
219  Point3df p;
220  p[0] = t5[0] + dx;
221  p[1] = t5[1] + dy;
222  p[2] = t5[2] + dz;
223  map[p] = true;
224  }
225  }
226  }
227  float nx = -diff[0][1] * diff[1][2] + diff[0][2] * diff[1][1];
228  float ny = -diff[0][2] * diff[1][0] + diff[0][0] * diff[1][2];
229  float nz = -diff[0][0] * diff[1][1] + diff[0][1] * diff[1][0];
230  float norm = sqrt(nx * nx + ny * ny + nz * nz);
231  nx /= norm;
232  ny /= norm;
233  nz /= norm;
234  float nx2 = nx * (nx > 0 ? 1 : -1);
235  float ny2 = ny * (ny > 0 ? 1 : -1);
236  float nz2 = nz * (nz > 0 ? 1 : -1);
237  float dlnk;
238  float nsum = (nx2 + ny2 + nz2);
239  if (connexity == 26)
240  // dlnk = sqrt(3) * cos(alpha) / 2.
241  // = sqrt(3) * <N, K > / (||K|| * 2.)
242  // = <N, [1, 1, 1] > / 2
243  dlnk = tc * nsum;
244  else if (connexity == 6)
245  {
246  // dlkn = cos(beta) / 2 = cos(pi/4. - alpha) / 2
247  // = sqrt(2) * (cos(alpha)+sin(alpha))/4
248  // = sqrt(2) * (cos(alpha)+sqrt(1-cos(alpha)^2))/4
249  float cosalpha = nsum / sqrt(3.);
250  dlnk = tc * (cosalpha + sqrt(1 - cosalpha * cosalpha));
251  }
252  std::map<Point3df,bool>::const_iterator i, e;
253  for (i = map.begin(), e = map.end(); i != e; ++i)
254  {
255  const Point3df &p0 = (*i).first;
256  Point3df p = p0 - t[0];
257  float d = nx * p[0] + ny * p[1] + nz * p[2];
258  if (d < 0) d = -d;
259  if (d <= dlnk)
260  this->set(output, (uint) p0[0],
261  (uint) p0[1],
262  (uint) p0[2], ind);
263  }
264  }
265  return output;
266 }
267 
268 };
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
O doit(const AimsSurfaceTriangle &surface, float spacing=1., unsigned int connexity=26)
surface : input mesh spacing : voxel size (according to mesh metric) connexity : 6 or 26 ...
unsigned int uint
AIMSDATA_API float norm(const Tensor &thing)
std::vector< float > toVector() const