aimsalgo 6.0.0
Neuroimaging image processing
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
34#include <aims/resampling/meshToVoxelsResampler.h>
35#include <aims/distancemap/stlsort.h> //for Point3dfCompare
36#include <cartobase/object/object.h>
37#include <aims/resampling/motion.h>
38#include <map>
39#include <float.h>
40
41namespace aims
42{
43
44template<typename O> void MeshToVoxelsResampler<O>::fill_header(
45carto::PropertySet &hdr, O &output, const AimsVector<float,3> & offset, float spacing)
46const
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.setTranslation( 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 setVoxelSize( output, spacing, spacing, spacing, 1. );
60}
61
67template<typename O>
68O MeshToVoxelsResampler<O>::doit(const AimsSurfaceTriangle &m, float spacing,
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 = 0.0f;
77 float tc = 0.0f;
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
116 AimsVector<float, 3> t[3];
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
126 AimsVector<float, 3> t1;
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 = 0.0f;
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};
aims::MeshToVoxelsResampler::doit
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
Definition meshToVoxelsResampler_d.h:68
carto::PropertySet::setProperty
void setProperty(const std::string &, const T &)
soma::AffineTransformation3dBase::setTranslation
void setTranslation(Point3df trans)
soma::AffineTransformationBase::toVector
std::vector< float > toVector() const
Motion
aims::AffineTransformation3d Motion
AimsSurfaceTriangle
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
norm
AIMSDATA_API float norm(const Tensor &thing)
uint
unsigned int uint
Point3df
AimsVector< float, 3 > Point3df