# Redistancing of a perturbed distance field

Test case originally written by Alexandre Limare.

This case is extracted from Russo et al., 2000. We initialize a perturbed distance field, where the zero levelset is an ellipse of the form \displaystyle \phi (x,y,0) = f(x,y) \times g(x,y) where the perturbation is \displaystyle f(x,y) = \epsilon + (x - x_0)^2 +(y - y_0)^2 and the ellipse is \displaystyle g(x,y) = \left( \sqrt{\frac{x^2}{A^2}+\frac{y^2}{B^2}} -R \right) with A=2, B=1, R = 1 , x_0 = 3.5, y_0 = 2..

We want to recover a perfect distance field, i.e. remove the initial perturbation.

We show here the initial and final level-set for the same isovalues.

f(x) = a + b*x
f1(x) = a1 + b1*x
f2(x) = a2 + b2*x
fit f(x) 'log' u (log($1)):(log($2)) via a,b
fit f1(x) 'log' u (log($1)):(log($3)) via a1,b1
fit f2(x) 'log' u (log($1)):(log($4)) via a2,b2
ftitle(a,b) = sprintf("%.2f-%4.2f x", a, -b)
set logscale xy
set xrange [32:512]
set xtics 32,2,512
set format y "%.1e"
plot 'log' u 1:2 t 'avg', exp(f(log(x))) t ftitle(a,b), \
'log' u 1:3 t 'rms', exp(f1(log(x))) t ftitle(a1,b1), \
'log' u 1:4 t 'max', exp(f2(log(x))) t ftitle(a2,b2)
f(x) = a + b*x
f1(x) = a1 + b1*x
f2(x) = a2 + b2*x
unset logscale
unset xrange
fit f(x) 'log' u (log($1)):(log($5)) via a,b
fit f1(x) 'log' u (log($1)):(log($6)) via a1,b1
fit f2(x) 'log' u (log($1)):(log($7)) via a2,b2
ftitle(a,b) = sprintf("%.2f-%4.2f x", a, -b)
set logscale xy
set xrange [32:512]
set xtics 32,2,512
set format y "%.1e"
plot 'log' u 1:5 t 'avg', exp(f(log(x))) t ftitle(a,b), \
'log' u 1:6 t 'rms', exp(f1(log(x))) t ftitle(a1,b1), \
'log' u 1:7 t 'max', exp(f2(log(x))) t ftitle(a2,b2)

Here we study the value of the level-set function on a set of points where it is theoretically 0, we show that we have also a 3^{rd} convergence.

## References

 [russo2000] Giovanni Russo and Peter Smereka. A remark on computing distance functions. Journal of Computational Physics, 163(1):51–67, 2000.
#include "distance_point_ellipse.h"
#include "redistance.h"
#include "view.h"

double perturb (double x, double y, double eps, coord center)
{
return eps + sq(x - center.x) + sq(y - center.y);
}

void draw_isolines (scalar s, double smin, double smax, int niso, int w)
{
vertex scalar vdist[];
foreach_vertex()
vdist[] = (s[] + s[-1] + s[0,-1] + s[-1,-1])/4.;
for (double sval = smin ; sval <= smax; sval += (smax - smin)/niso)
isoline ("vdist", sval, lw = w);
}

int main()
{
origin (-5., -5.);
L0 = 10;

for (int MAXLEVEL = 6; MAXLEVEL < 9; MAXLEVEL++){
init_grid (1 << MAXLEVEL);

scalar dist[];
double A = 4., B = 2.;
coord  center_perturb = {3.5, 2.};
foreach(){
double a, b;
dist[] = DistancePointEllipse (A, B, x, y, &a, &b)*
perturb (x, y, 0.1, center_perturb);
}

if (MAXLEVEL == 8) {
squares ("dist", map = cool_warm, min = -2, max = 2);
draw_isolines (dist, -2, 2, 20, 1);
save ("start.png");
}

int nbit = redistance (dist, imax = 1 << (MAXLEVEL + 1));

if (MAXLEVEL == 8) {
squares ("dist", map = cool_warm, min = -2, max = 2);
draw_isolines (dist, -2, 2, 20, 1);
save ("final.png");
}

scalar err[], errt[], errt1[];
foreach() {
double a, b;
err[] = dist[] - DistancePointEllipse (A, B, x, y, &a, &b);
errt[] = dist[] > - 0.8 ? err[] : nodata;
errt1[] = fabs (dist[]) < 1.2*L0/(1 << grid->maxdepth) ? err[] : nodata;
}

norm n = normf (errt), n2 = normf (errt1);
fprintf (stderr, "%d %g %g %g %g %g %g %d\n",
1 << MAXLEVEL, n.avg, n.rms, n.max,
n2.avg, n2.rms, n2.max, nbit);

if (MAXLEVEL == 8) {
squares ("log(fabs(err) + 1e-16)", min = log(n.max) - 6, max = log(n.max));
save ("err.png");
}
}
}