Imagine++
Schemes.h
1 // ===========================================================================
2 // Imagine++ Libraries
3 // Copyright (C) Imagine
4 // For detailed information: http://imagine.enpc.fr/software
5 // ===========================================================================
6 
7 namespace 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>
25  void neighbourOffsets(const Image<T,dim>& u, const Coords<dim>& p, FVector<size_t,dim>& dp,FVector<size_t,dim>& dm) {
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];
95  FVector<size_t,3> dp,dm;
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>
152  T meanCurvatureMotion(const Image<T,3>& u, const Coords<3>& p) {
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>
238  T meanCurvatureMotion(const Image<T,2>& u, const Coords<2>& p) {
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>
334  FVector<T,dim> normal(const Image<T,dim>& u, const Coords<dim>& p) {
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>
350  FVector<T,dim> gradient(const Image<T,dim> &u, const Coords<dim>& p) {
351  const size_t o = u.offset(p);
352  const T u0 = u[o];
353  FVector<T,dim> g;
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 }
Imagine::Coords
Coordinates.
Definition: Coords.h:16
Imagine::FVector
Vector of fixed size.
Definition: FVector.h:17
Imagine::Image::scalar_type
PixelTraits< T >::scalar_type scalar_type
Scalar type.
Definition: Image.h:32
Imagine::meanCurvatureMotion
T meanCurvatureMotion(const Image< T, 3 > &u, const Coords< 3 > &p)
Level set Mean curvature motion (3D).
Definition: Schemes.h:152
Imagine::MultiArray::size
int size(int i) const
ith size.
Definition: MultiArray.h:225
Imagine::neighbourOffsets
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::Image
Image.
Definition: Image.h:24
Imagine::MultiArray::offset
size_t offset(const Coords< dim > &c) const
Offset.
Definition: MultiArray.h:283
Imagine::normal
FVector< T, dim > normal(const Image< T, dim > &u, const Coords< dim > &p)
Unit normal of iso level.
Definition: Schemes.h:334
Imagine::MultiArray::stride
FArray< size_t, dim > stride() const
Stride.
Definition: MultiArray.h:253
Imagine::neighbourCoords
void neighbourCoords(const Image< T, dim > &u, const Coords< dim > &p, Coords< dim > &pp, Coords< dim > &pm)
Coordinates of neighbours.
Definition: Schemes.h:41
Imagine::gradient
FVector< T, dim > gradient(const Image< T, dim > &u, const Coords< dim > &p)
Gradient.
Definition: Schemes.h:350
Imagine::laplacian
T laplacian(const Image< T, dim > &u, const Coords< dim > &p)
Laplacian.
Definition: Schemes.h:74
Imagine
Imagine++ namespace.
Definition: Array.h:7
Imagine::gaussianCurvature
T gaussianCurvature(const Image< T, 3 > &u, const Coords< 3 > &p)
Gaussian curvature of iso level (3D).
Definition: Schemes.h:278
Imagine::meanCurvature
T meanCurvature(const Image< T, 3 > &u, const Coords< 3 > &p)
Mean curvature (3D).
Definition: Schemes.h:92
Imagine::derivative
T derivative(const Image< T, dim > &u, const Coords< dim > &p, int d)
Derivative.
Definition: Schemes.h:58
Generated by   doxygen 1.8.13