00001
00009 #ifndef SCHEMES_H
00010 #define SCHEMES_H
00011
00012 #include <Globals.h>
00013 #include <Image.H>
00014
00015 using namespace Images;
00016
00017 namespace levelset {
00018
00023 template<unsigned DIM, typename Pixel = float>
00024 class Schemes
00025 {
00026 public:
00027 typedef BaseImage<DIM, Pixel> ImageType;
00028 typedef typename ImageType::Index Index;
00029
00032
00034 static void find_neighbours(const ImageType & phi, const Index & ind, int inf_coords[], int sup_coords[]);
00036 static Pixel laplacian(const ImageType & phi, const Index & ind);
00038 static Pixel mean_curvature(const ImageType & phi, const Index & ind);
00040 static void normal(const ImageType & phi, const Index & ind, Pixel norm[]);
00042 static void gradient(const ImageType & phi, const Index & ind, Pixel grad[]);
00044 static float sq_norm_gradient(const ImageType & phi, const Index & ind);
00046 static float norm_gradient(const ImageType & phi, const Index & ind);
00048 static Pixel other_mean_curvature(const ImageType & phi, const Index & ind);
00049
00051 };
00052
00053
00055
00056
00057
00058 template<unsigned DIM, typename Pixel>
00059 void Schemes<DIM, Pixel>::find_neighbours(const ImageType & phi, const Index & ind, int inf_coords[], int sup_coords[])
00060 {
00061 for (int i=0 ; i<DIM ; i++){
00062 inf_coords[i] = (ind(i+1)>0) ? ind(i+1) - 1 : ind(i+1);
00063 sup_coords[i] = (ind(i+1)<phi.extent(i)-1) ? ind(i+1) + 1 : ind(i+1);
00064 }
00065 }
00066
00067 template <unsigned DIM, typename Pixel>
00068 Pixel Schemes<DIM, Pixel>::laplacian(const ImageType & phi, const Index & ind)
00069 {
00070 Pixel result = DIM * phi(ind);
00071 result *= -2;
00072
00073 int inf_coords[DIM];
00074 int sup_coords[DIM];
00075 Schemes::find_neighbours(phi, ind, inf_coords, sup_coords);
00076
00077 Index sup;
00078 Index inf;
00079 for (int i=0 ; i<DIM ; i++){
00080 sup = ind;
00081 inf = ind;
00082 sup(i+1) = sup_coords[i];
00083 inf(i+1) = inf_coords[i];
00084 result += phi(sup) + phi(inf);
00085 }
00086
00087 return result;
00088 }
00089
00090 template <unsigned DIM, typename Pixel>
00091 Pixel Schemes<DIM, Pixel>::mean_curvature(const ImageType & phi, const Index & ind)
00092 {
00093 int inf_coords[DIM];
00094 int sup_coords[DIM];
00095 Schemes::find_neighbours(phi, ind, inf_coords, sup_coords);
00096
00097 const blitz::Array<Pixel,DIM> & array = (* dynamic_cast<const blitz::Array<Pixel,DIM> *>(&phi));
00098
00099 if (DIM == 2){
00100
00101 int x = ind(1);
00102 int y = ind(2);
00103 int mx = inf_coords[0];
00104 int px = sup_coords[0];
00105 int my = inf_coords[1];
00106 int py = sup_coords[1];
00107
00108 const Pixel u0 = array(x,y);
00109
00110 const Pixel upx = array(px,y);
00111 const Pixel umx = array(mx,y);
00112 const Pixel upy = array(x,py);
00113 const Pixel umy = array(x,my);
00114
00115 const Pixel umxmy = array(mx,my);
00116 const Pixel upxmy = array(px,my);
00117 const Pixel umxpy = array(mx,py);
00118 const Pixel upxpy = array(px,py);
00119
00120
00121 const Pixel ux = ( upx - umx ) / 2;
00122 const Pixel uy = ( upy - umy ) / 2;
00123
00124 const Pixel uxx = upx - 2 * u0 + umx;
00125 const Pixel uyy = upy - 2 * u0 + umy;
00126 const Pixel uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
00127
00128 const Pixel ux2 = ux * ux;
00129 const Pixel uy2 = uy * uy;
00130 const Pixel grad = ux2 + uy2;
00131
00132 if (grad < Globals::EPS) return 0;
00133
00134 return (uyy * ux2 - 2 * ux * uy * uxy + uxx * uy2) / grad / sqrt(grad);
00135 }
00136
00137 else if (DIM == 3){
00138
00139 int x = ind(1);
00140 int y = ind(2);
00141 int z = ind(3);
00142 int mx = inf_coords[0];
00143 int px = sup_coords[0];
00144 int my = inf_coords[1];
00145 int py = sup_coords[1];
00146 int mz = inf_coords[2];
00147 int pz = sup_coords[2];
00148
00149 const Pixel u0 = array(x,y,z);
00150 const Pixel upx = array(px,y,z);
00151 const Pixel umx = array(mx,y,z);
00152 const Pixel upy = array(x,py,z);
00153 const Pixel umy = array(x,my,z);
00154 const Pixel upz = array(x,y,pz);
00155 const Pixel umz = array(x,y,mz);
00156
00157 const Pixel umxmy = array(mx,my,z);
00158 const Pixel upxmy = array(px,my,z);
00159 const Pixel umxpy = array(mx,py,z);
00160 const Pixel upxpy = array(px,py,z);
00161 const Pixel umymz = array(x,my,mz);
00162 const Pixel upymz = array(x,py,mz);
00163 const Pixel umypz = array(x,my,pz);
00164 const Pixel upypz = array(x,py,pz);
00165 const Pixel umzmx = array(mx,y,mz);
00166 const Pixel upzmx = array(mx,y,pz);
00167 const Pixel umzpx = array(px,y,mz);
00168 const Pixel upzpx = array(px,y,pz);
00169
00170 const Pixel ux = ( upx - umx ) / 2;
00171 const Pixel uy = ( upy - umy ) / 2;
00172 const Pixel uz = ( upz - umz ) / 2;
00173
00174 const Pixel uxx = upx - 2 * u0 + umx;
00175 const Pixel uyy = upy - 2 * u0 + umy;
00176 const Pixel uzz = upz - 2 * u0 + umz;
00177 const Pixel uxy = (upxpy + umxmy - upxmy - umxpy) / 4;
00178 const Pixel uyz = (upypz + umymz - upymz - umypz) / 4;
00179 const Pixel uzx = (upzpx + umzmx - upzmx - umzpx) / 4;
00180
00181 const Pixel ux2 = ux * ux;
00182 const Pixel uy2 = uy * uy;
00183 const Pixel uz2 = uz * uz;
00184 const Pixel grad = ux2 + uy2 + uz2;
00185
00186 if (grad < Globals::EPS) return 0;
00187
00188 return ( (uyy+uzz) * ux2
00189 + (uxx+uzz) * uy2
00190 + (uxx+uyy) * uz2
00191 - 2*ux*uy*uxy
00192 - 2*uz*ux*uzx
00193 - 2*uy*uz*uyz ) / grad / sqrt(grad) / 2;
00194 }
00195
00196 else std::cerr << "Numerical operation not implemented in this dimension." << std::endl;
00197
00198 return 0;
00199 }
00200
00201 template <unsigned DIM, typename Pixel>
00202 void Schemes<DIM, Pixel>::normal(const ImageType & phi, const Index & ind, Pixel norm[])
00203 {
00204 Schemes::gradient(phi, ind, norm);
00205
00206 float n = 0;
00207
00208 for (int i=0 ; i<DIM ; i++)
00209 n += Maths::sqr(norm[i]);
00210
00211 n = sqrt( n );
00212 if (n < Globals::EPS) n = (Pixel) Globals::EPS;
00213
00214 for (int i=0 ; i<DIM ; i++)
00215 norm[i] /= n;
00216 }
00217
00218 template <unsigned DIM, typename Pixel>
00219 void Schemes<DIM, Pixel>::gradient(const ImageType & phi, const Index & ind, Pixel grad[])
00220 {
00221 int inf_coords[DIM];
00222 int sup_coords[DIM];
00223 Schemes::find_neighbours(phi, ind, inf_coords, sup_coords);
00224
00225 Index inf;
00226 Index sup;
00227 for (int i=0 ; i<DIM ; i++){
00228 inf = ind;
00229 sup = ind;
00230 inf(i+1) = inf_coords[i];
00231 sup(i+1) = sup_coords[i];
00232 grad[i] = (phi(sup) - phi(inf)) / 2.0;
00233 }
00234 }
00235
00236 template <unsigned DIM, typename Pixel>
00237 float Schemes<DIM, Pixel>::sq_norm_gradient(const ImageType & phi, const Index & ind)
00238 {
00239 float result = 0;
00240 Pixel grad[DIM];
00241 Schemes::gradient(phi, ind, grad);
00242
00243 for (int i=0 ; i<DIM ; i++)
00244 result += Maths::sqr(grad[i]);
00245
00246 return result;
00247 }
00248
00249 template <unsigned DIM, typename Pixel>
00250 float Schemes<DIM, Pixel>::norm_gradient(const ImageType & phi, const Index & ind)
00251 {
00252 return sqrt( Schemes::sq_norm_gradient(phi, ind));
00253 }
00254
00255 template <unsigned DIM, typename Pixel>
00256 Pixel Schemes<DIM, Pixel>::other_mean_curvature(const ImageType & phi, const Index & ind)
00257 {
00258
00259 int inf_coords[DIM];
00260 int sup_coords[DIM];
00261 Schemes::find_neighbours(phi, ind, inf_coords, sup_coords);
00262
00263 const blitz::Array<Pixel,DIM> & array = (* dynamic_cast<const blitz::Array<Pixel,DIM> *>(&phi));
00264
00265 if (DIM == 2){
00266
00267 int x = ind(1);
00268 int y = ind(2);
00269 int mx = inf_coords[0];
00270 int px = sup_coords[0];
00271 int my = inf_coords[1];
00272 int py = sup_coords[1];
00273 int ppx = (px<phi.extent(0)-1) ? px+1:px;
00274 int mmx = (mx>0) ? mx-1:mx;
00275 int ppy = (py<phi.extent(1)-1) ? py+1:py;
00276 int mmy = (my>0) ? my-1:my;
00277
00278 const Pixel u0 = array(x,y);
00279 const Pixel upx = array(px,y);
00280 const Pixel umx = array(mx,y);
00281 const Pixel upy = array(x,py);
00282 const Pixel umy = array(x,my);
00283 const Pixel uppx = array(ppx,y);
00284 const Pixel ummx = array(mmx,y);
00285 const Pixel uppy = array(x,ppy);
00286 const Pixel ummy = array(x,mmy);
00287 const Pixel umxmy = array(mx,my);
00288 const Pixel upxmy = array(px,my);
00289 const Pixel umxpy = array(mx,py);
00290 const Pixel upxpy = array(px,py);
00291
00292 const Pixel upx_x = ( uppx - u0 ) / 2;
00293 const Pixel upx_y = ( upxpy - upxmy ) / 2;
00294 const Pixel umx_x = ( u0 - ummx ) / 2;
00295 const Pixel umx_y = ( umxpy - umxmy ) / 2;
00296 const Pixel upy_x = ( upxpy - umxpy ) / 2;
00297 const Pixel upy_y = ( uppy - u0 ) / 2;
00298 const Pixel umy_x = ( upxmy - umxmy ) / 2;
00299 const Pixel umy_y = ( u0 - ummy ) / 2;
00300
00301 const Pixel gradpx = upx_x*upx_x + upx_y*upx_y;
00302 const Pixel gradmx = umx_x*umx_x + umx_y*umx_y;
00303 const Pixel gradpy = upy_x*upy_x + upy_y*upy_y;
00304 const Pixel gradmy = umy_x*umy_x + umy_y*umy_y;
00305
00306 const Pixel normalpx_x = upx_x / sqrt(gradpx + 1e-10);
00307 const Pixel normalmx_x = umx_x / sqrt(gradmx + 1e-10);
00308 const Pixel normalpy_y = upy_y / sqrt(gradpy + 1e-10);
00309 const Pixel normalmy_y = umy_y / sqrt(gradmy + 1e-10);
00310
00311 const Pixel normal_xx = ( normalpx_x - normalmx_x ) / 2;
00312 const Pixel normal_yy = ( normalpy_y - normalmy_y ) / 2;
00313
00314 return normal_xx + normal_yy;
00315 }
00316
00317 else std::cerr << "Numerical operation not yet implemented in this dimension." << std::endl;
00318
00319 return 0;
00320 }
00321 }
00322
00323
00324
00325
00326 #endif // SCHEMES_H