Imagine++
Schemes.h
1// ===========================================================================
2// Imagine++ Libraries
3// Copyright (C) Imagine
4// For detailed information: http://imagine.enpc.fr/software
5// ===========================================================================
6
7namespace Imagine {
10
11#ifndef DOXYGEN_SHOULD_SKIP_THIS
12 const double IMAGE_SCHEME_EPS = 1e-6;
13#endif
14
24 template <typename T,int dim>
26 for( int i = 0 ; i < dim ; i++) {
27 dm[i] = (p[i] > 0) ? u.stride(i) : 0;
28 dp[i] = (p[i] < u.size(i)-1) ? u.stride(i) : 0;
29 }
30 }
40 template <typename T,int dim>
41 void neighbourCoords(const Image<T,dim>& u, const Coords<dim>& p, Coords<dim>& pp, Coords<dim>& pm) {
42 for( int i = 0 ; i < dim ; i++) {
43 pm[i] = (p[i] > 0) ? p[i]-1 : p[i];
44 pp[i] = (p[i] < u.size(i)-1) ? p[i]+1 : p[i];
45 }
46 }
57 template <typename T, int dim>
58 T derivative(const Image<T,dim>& u, const Coords<dim>& p, int d) {
59 Coords<dim> np(p), pp(p);
60 typename Image<T,dim>::scalar_type h(2);
61 if (np[d]>0) np[d]--; else h--;
62 if (pp[d]<u.size(d)-1) pp[d]++; else h--;
63 return ( u(pp) - u(np) ) / h;
64 }
73 template <typename T, int dim>
74 T laplacian(const Image<T,dim>& u, const Coords<dim>& p) {
75 const size_t o = u.offset(p);
76 const T u0 = u[o];
77 T l = -typename Image<T,dim>::scalar_type(2*dim)*u0;
78 for (int i=0;i<dim;i++) {
79 l += ( ( p[i] < u.size(i)-1 ) ? u[ o + u.stride(i) ] : u0 ) + ( ( p[i] > 0 ) ? u[ o - u.stride(i) ] : u0 ) ;
80 }
81 return l;
82 }
91 template <typename T>
92 T meanCurvature(const Image<T,3>& u, const Coords<3>& p) {
93 const size_t o = u.offset(p);
94 const T u0 = u[o];
96 neighbourOffsets(u,p,dp,dm);
97
98 const T upx = u[o+dp[0]];
99 const T umx = u[o-dm[0]];
100 const T upy = u[o+dp[1]];
101 const T umy = u[o-dm[1]];
102 const T upz = u[o+dp[2]];
103 const T umz = u[o-dm[2]];
104
105 const T umxmy = u[o-dm[0]-dm[1]];
106 const T upxmy = u[o+dp[0]-dm[1]];
107 const T umxpy = u[o-dm[0]+dp[1]];
108 const T upxpy = u[o+dp[0]+dp[1]];
109 const T umymz = u[o-dm[1]-dm[2]];
110 const T upymz = u[o+dp[1]-dm[2]];
111 const T umypz = u[o-dm[1]+dp[2]];
112 const T upypz = u[o+dp[1]+dp[2]];
113 const T umzmx = u[o-dm[2]-dm[0]];
114 const T upzmx = u[o+dp[2]-dm[0]];
115 const T umzpx = u[o-dm[2]+dp[0]];
116 const T upzpx = u[o+dp[2]+dp[0]];
117
118 const T ux = ( upx - umx ) / 2;
119 const T uy = ( upy - umy ) / 2;
120 const T uz = ( upz - umz ) / 2;
121
122 const T uxx = upx - 2 * u0 + umx;
123 const T uyy = upy - 2 * u0 + umy;
124 const T uzz = upz - 2 * u0 + umz;
125 const T uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
126 const T uyz = (upypz + umymz - upymz - umypz) / 4;
127 const T uzx = (upzpx + umzmx - upzmx - umzpx) / 4;
128
129 const T ux2 = ux * ux;
130 const T uy2 = uy * uy;
131 const T uz2 = uz * uz;
132 const T grad = ux2 + uy2 + uz2;
133
134 if (grad < T(IMAGE_SCHEME_EPS)) return 0;
135
136 return ( (uyy+uzz) * ux2
137 + (uxx+uzz) * uy2
138 + (uxx+uyy) * uz2
139 - 2*ux*uy*uxy
140 - 2*uz*ux*uzx
141 - 2*uy*uz*uyz ) / grad / std::sqrt(grad) / T(2);
142 }
151 template <typename T>
153 const size_t o = u.offset(p);
154 const T u0 = u[o];
155 FVector<size_t,3> dp,dm;
156 neighbourOffsets(u,p,dp,dm);
157
158 const T upx = u[o+dp[0]];
159 const T umx = u[o-dm[0]];
160 const T upy = u[o+dp[1]];
161 const T umy = u[o-dm[1]];
162 const T upz = u[o+dp[2]];
163 const T umz = u[o-dm[2]];
164
165 const T umxmy = u[o-dm[0]-dm[1]];
166 const T upxmy = u[o+dp[0]-dm[1]];
167 const T umxpy = u[o-dm[0]+dp[1]];
168 const T upxpy = u[o+dp[0]+dp[1]];
169 const T umymz = u[o-dm[1]-dm[2]];
170 const T upymz = u[o+dp[1]-dm[2]];
171 const T umypz = u[o-dm[1]+dp[2]];
172 const T upypz = u[o+dp[1]+dp[2]];
173 const T umzmx = u[o-dm[2]-dm[0]];
174 const T upzmx = u[o+dp[2]-dm[0]];
175 const T umzpx = u[o-dm[2]+dp[0]];
176 const T upzpx = u[o+dp[2]+dp[0]];
177
178 const T ux = ( upx - umx ) / 2;
179 const T uy = ( upy - umy ) / 2;
180 const T uz = ( upz - umz ) / 2;
181
182 const T uxx = upx - 2 * u0 + umx;
183 const T uyy = upy - 2 * u0 + umy;
184 const T uzz = upz - 2 * u0 + umz;
185 const T uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
186 const T uyz = (upypz + umymz - upymz - umypz) / 4;
187 const T uzx = (upzpx + umzmx - upzmx - umzpx) / 4;
188
189 const T ux2 = ux * ux;
190 const T uy2 = uy * uy;
191 const T uz2 = uz * uz;
192 const T grad = ux2 + uy2 + uz2;
193
194 if (grad < T(IMAGE_SCHEME_EPS)) return 0;
195
196 return ( (uyy+uzz) * ux2
197 + (uxx+uzz) * uy2
198 + (uxx+uyy) * uz2
199 - 2*ux*uy*uxy
200 - 2*uz*ux*uzx
201 - 2*uy*uz*uyz ) / grad / T(2);
202 }
204 template <typename T>
205 T meanCurvature(const Image<T,2>& u, const Coords<2>& p) {
206 const size_t o = u.offset(p);
207 const T u0 = u[o];
208 FVector<size_t,2> dp,dm;
209 neighbourOffsets(u,p,dp,dm);
210
211 const T upx = u[o+dp[0]];
212 const T umx = u[o-dm[0]];
213 const T upy = u[o+dp[1]];
214 const T umy = u[o-dm[1]];
215
216 const T umxmy = u[o-dm[0]-dm[1]];
217 const T upxmy = u[o+dp[0]-dm[1]];
218 const T umxpy = u[o-dm[0]+dp[1]];
219 const T upxpy = u[o+dp[0]+dp[1]];
220
221 const T ux = ( upx - umx ) / 2;
222 const T uy = ( upy - umy ) / 2;
223
224 const T uxx = upx - 2 * u0 + umx;
225 const T uyy = upy - 2 * u0 + umy;
226 const T uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
227
228 const T ux2 = ux * ux;
229 const T uy2 = uy * uy;
230 const T grad = ux2 + uy2;
231
232 if (grad < T(IMAGE_SCHEME_EPS)) return 0;
233
234 return (uyy * ux2 - 2 * ux * uy * uxy + uxx * uy2) / grad / std::sqrt(grad);
235 }
237 template <typename T>
239 const size_t o = u.offset(p);
240 const T u0 = u[o];
241 FVector<size_t,2> dp,dm;
242 neighbourOffsets(u,p,dp,dm);
243
244 const T upx = u[o+dp[0]];
245 const T umx = u[o-dm[0]];
246 const T upy = u[o+dp[1]];
247 const T umy = u[o-dm[1]];
248
249 const T umxmy = u[o-dm[0]-dm[1]];
250 const T upxmy = u[o+dp[0]-dm[1]];
251 const T umxpy = u[o-dm[0]+dp[1]];
252 const T upxpy = u[o+dp[0]+dp[1]];
253
254 const T ux = ( upx - umx ) / 2;
255 const T uy = ( upy - umy ) / 2;
256
257 const T uxx = upx - 2 * u0 + umx;
258 const T uyy = upy - 2 * u0 + umy;
259 const T uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
260
261 const T ux2 = ux * ux;
262 const T uy2 = uy * uy;
263 const T grad = ux2 + uy2;
264
265 if (grad < T(IMAGE_SCHEME_EPS)) return 0;
266
267 return (uyy * ux2 - 2 * ux * uy * uxy + uxx * uy2) / grad;
268 }
277 template <typename T>
278 T gaussianCurvature(const Image<T,3>& u, const Coords<3>& p) {
279 const size_t o = u.offset(p);
280 const T u0 = u[o];
281 FVector<size_t,3> dp,dm;
282 neighbourOffsets(u,p,dp,dm);
283
284 const T upx = u[o+dp[0]];
285 const T umx = u[o-dm[0]];
286 const T upy = u[o+dp[1]];
287 const T umy = u[o-dm[1]];
288 const T upz = u[o+dp[2]];
289 const T umz = u[o-dm[2]];
290
291 const T umxmy = u[o-dm[0]-dm[1]];
292 const T upxmy = u[o+dp[0]-dm[1]];
293 const T umxpy = u[o-dm[0]+dp[1]];
294 const T upxpy = u[o+dp[0]+dp[1]];
295 const T umymz = u[o-dm[1]-dm[2]];
296 const T upymz = u[o+dp[1]-dm[2]];
297 const T umypz = u[o-dm[1]+dp[2]];
298 const T upypz = u[o+dp[1]+dp[2]];
299 const T umzmx = u[o-dm[2]-dm[0]];
300 const T upzmx = u[o+dp[2]-dm[0]];
301 const T umzpx = u[o-dm[2]+dp[0]];
302 const T upzpx = u[o+dp[2]+dp[0]];
303
304 const T ux = ( upx - umx ) / 2;
305 const T uy = ( upy - umy ) / 2;
306 const T uz = ( upz - umz ) / 2;
307
308 const T uxx = upx - 2 * u0 + umx;
309 const T uyy = upy - 2 * u0 + umy;
310 const T uzz = upz - 2 * u0 + umz;
311 const T uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
312 const T uyz = (upypz + umymz - upymz - umypz) / 4;
313 const T uzx = (upzpx + umzmx - upzmx - umzpx) / 4;
314
315 const T ux2 = ux * ux;
316 const T uy2 = uy * uy;
317 const T uz2 = uz * uz;
318 const T grad = ux2 + uy2 + uz2;
319
320 if (grad < T(IMAGE_SCHEME_EPS)) return 0;
321
322 return ( ux2 * (uyy*uzz - uyz*uyz) + uy2 * (uxx*uzz - uzx*uzx) + uz2 * (uxx*uyy - uxy*uxy)
323 + 2 * ( ux*uy * (uzx*uyz - uxy*uzz) + uy*uz * (uxy*uzx - uyz*uxx) + ux*uz * (uxy*uyz - uzx*uyy) ) ) / (grad*grad);
324 }
333 template <typename T, int dim>
335 FVector<T,dim> n = gradient(u,p);
336 T nn = norm(n);
337 if (nn < T(IMAGE_SCHEME_EPS)) nn = T(IMAGE_SCHEME_EPS);
338 n /= nn;
339 return n;
340 }
349 template <typename T, int dim>
351 const size_t o = u.offset(p);
352 const T u0 = u[o];
354 for (int i=0;i<dim;i++) {
355 g[i] = ( ( ( p[i] < u.size(i)-1 ) ? u[ o + u.stride(i) ] : u0 ) - ( ( p[i] > 0 ) ? u[ o - u.stride(i) ] : u0 ) ) / T(2);
356 }
357 return g;
358 }
359
361}
Coordinates.
Definition: Coords.h:16
Vector of fixed size.
Definition: FVector.h:17
Image.
Definition: Image.h:24
PixelTraits< T >::scalar_type scalar_type
Scalar type.
Definition: Image.h:32
FArray< size_t, dim > stride() const
Stride.
Definition: MultiArray.h:253
size_t offset(const Coords< dim > &c) const
Offset.
Definition: MultiArray.h:283
int size(int i) const
ith size.
Definition: MultiArray.h:225
void neighbourCoords(const Image< T, dim > &u, const Coords< dim > &p, Coords< dim > &pp, Coords< dim > &pm)
Coordinates of neighbours.
Definition: Schemes.h:41
T meanCurvature(const Image< T, 3 > &u, const Coords< 3 > &p)
Mean curvature (3D).
Definition: Schemes.h:92
T laplacian(const Image< T, dim > &u, const Coords< dim > &p)
Laplacian.
Definition: Schemes.h:74
T meanCurvatureMotion(const Image< T, 3 > &u, const Coords< 3 > &p)
Level set Mean curvature motion (3D).
Definition: Schemes.h:152
FVector< T, dim > gradient(const Image< T, dim > &u, const Coords< dim > &p)
Gradient.
Definition: Schemes.h:350
FVector< T, dim > normal(const Image< T, dim > &u, const Coords< dim > &p)
Unit normal of iso level.
Definition: Schemes.h:334
T derivative(const Image< T, dim > &u, const Coords< dim > &p, int d)
Derivative.
Definition: Schemes.h:58
T gaussianCurvature(const Image< T, 3 > &u, const Coords< 3 > &p)
Gaussian curvature of iso level (3D).
Definition: Schemes.h:278
void neighbourOffsets(const Image< T, dim > &u, const Coords< dim > &p, FVector< size_t, dim > &dp, FVector< size_t, dim > &dm)
Offsets to neighbours.
Definition: Schemes.h:25
Imagine++ namespace.