5454
5555// Set version number of the library.
5656#ifndef cimg_version
57- #define cimg_version 377
57+ #define cimg_version 380
5858
5959/*-----------------------------------------------------------
6060 #
@@ -46156,6 +46156,8 @@ namespace cimg_library {
4615646156 /**
4615746157 \param reference Reference image R.
4615846158 \param smoothness Smoothness of estimated displacement field.
46159+ If smoothness is positive, Tikhonov regularization is applied, otherwise TV regularization is applied,
46160+ with specified strength (absolute value of the smoothnes).
4615946161 \param precision Precision required for algorithm convergence.
4616046162 \param nb_scales Number of scales used to estimate the displacement field.
4616146163 \param iteration_max Maximum number of iterations allowed for one scale.
@@ -46186,12 +46188,6 @@ namespace cimg_library {
4618646188 cimg_instance,
4618746189 reference._width,reference._height,reference._depth,reference._spectrum,
4618846190 reference._data);
46189- if (smoothness<0)
46190- throw CImgArgumentException(_cimg_instance
46191- "displacement(): Invalid specified smoothness %g "
46192- "(should be >=0)",
46193- cimg_instance,
46194- smoothness);
4619546191 if (precision<0)
4619646192 throw CImgArgumentException(_cimg_instance
4619746193 "displacement(): Invalid specified precision %g "
@@ -46289,7 +46285,7 @@ namespace cimg_library {
4628946285 }
4629046286
4629146287 // Regularization term.
46292- if (smoothness>0) cimg_forC(U,c) {
46288+ if (smoothness>0) cimg_forC(U,c) { // Tikhonov
4629346289 const double
4629446290 uccc = U(x,y,z,c),
4629546291 upcc = U(_p1x,y,z,c), uncc = U(_n1x,y,z,c),
@@ -46298,6 +46294,30 @@ namespace cimg_library {
4629846294 ux = 0.5*(uncc - upcc), uy = 0.5*(ucnc - ucpc), uz = 0.5*(uccn - uccp);
4629946295 energy+=smoothness*(ux*ux + uy*uy + uz*uz);
4630046296 _V[c]+=smoothness*(upcc + uncc + ucpc + ucnc + uccp + uccn - 6*uccc);
46297+ } else if (smoothness<0) cimg_forC(U,c) { // Total variation
46298+ const double
46299+ ucpp = U(x,_p1y,_p1z,c),
46300+ upcp = U(_p1x,y,_p1z,c), uccp = U(x,y,_p1z,c), uncp = U(_n1x,y,_p1z,c),
46301+ ucnp = U(x,_n1y,_p1z,c),
46302+ uppc = U(_p1x,_p1y,z,c), ucpc = U(x,_p1y,z,c), unpc = U(_n1x,_p1y,z,c),
46303+ upcc = U(_p1x,y,z,c), uccc = U(x,y,z,c), uncc = U(_n1x,y,z,c),
46304+ upnc = U(_p1x,_n1y,z,c), ucnc = U(x,_n1y,z,c), unnc = U(_n1x,_n1y,z,c),
46305+ ucpn = U(x,_p1y,_n1z,c),
46306+ upcn = U(_p1x,y,_n1z,c), uccn = U(x,y,_n1z,c), uncn = U(_n1x,y,_n1z,c),
46307+ ucnn = U(x,_n1y,_n1z,c),
46308+ ux = 0.5*(uncc - upcc), uy = 0.5*(ucnc - ucpc), uz = 0.5*(uccn - uccp),
46309+ gn = 1e-8 + std::sqrt(ux*ux + uy*uy + uz*uz),
46310+ nux = ux/gn, nuy = uy/gn, nuz = uz/gn,
46311+ uxx = uncc + upcc - 2*uccc,
46312+ uxy = 0.25*(uppc + unnc - upnc - unpc),
46313+ uxz = 0.25*(upcp + uncn - upcn - uncp),
46314+ uyy = ucnc + ucpc - 2*uccc,
46315+ uyz = 0.25*(ucpp + ucnn - ucpn - ucnp),
46316+ uzz = uccn + uccp - 2*uccc,
46317+ unn = nux*nux*uxx + nuy*nuy*uyy + nuz*nuz*uzz + 2*(nux*nuy*uxy + nux*nuz*uxz + nuy*nuz*uyz),
46318+ uee = uxx + uyy + uzz - unn;
46319+ energy-=smoothness*gn;
46320+ _V[c]-=smoothness*(uee/gn);
4630146321 }
4630246322
4630346323 // Guide term.
@@ -46338,14 +46358,28 @@ namespace cimg_library {
4633846358 }
4633946359
4634046360 // Regularization term.
46341- if (smoothness>0) cimg_forC(U,c) {
46361+ if (smoothness>0) cimg_forC(U,c) { // Tikhonov
4634246362 const double
4634346363 ucc = U(x,y,c),
4634446364 upc = U(_p1x,y,c), unc = U(_n1x,y,c),
4634546365 ucp = U(x,_p1y,c), ucn = U(x,_n1y,c),
4634646366 ux = 0.5*(unc - upc), uy = 0.5*(ucn - ucp);
4634746367 energy+=smoothness*(ux*ux + uy*uy);
4634846368 _V[c]+=smoothness*(upc + unc + ucp + ucn - 4*ucc);
46369+ } else if (smoothness<0) cimg_forC(U,c) { // Total variation
46370+ const double
46371+ upp = U(_p1x,_p1y,c), ucp = U(x,_p1y,c), unp = U(_n1x,_p1y,c),
46372+ upc = U(_p1x,y,c), ucc = U(x,y,c), unc = U(_n1x,y,c),
46373+ upn = U(_p1x,_n1y,c), ucn = U(x,_n1y,c), unn = U(_n1x,_n1y,c),
46374+ ux = 0.5*(unc - upc), uy = 0.5*(ucn - ucp),
46375+ gn = 1e-8 + std::sqrt(ux*ux + uy*uy),
46376+ nux = ux/gn, nuy = uy/gn,
46377+ uxx = unc + upc - 2*ucc,
46378+ uxy = 0.25*(upp + unn - upn - unp),
46379+ uyy = ucn + ucp - 2*ucc,
46380+ uee = nuy*nuy*uxx + nux*nux*uyy - 2*nux*nuy*uxy;
46381+ energy-=smoothness*gn;
46382+ _V[c]-=smoothness*(uee/gn);
4634946383 }
4635046384
4635146385 // Guide term.
0 commit comments