load "ppm2rnm" load "isoline"
string leman="lg.pgm";
real AreaLac =  580.03; // $Km^2$
real hsize= 5;
real[int,int] Curves(3,1);
int[int] be(1);
int nc;// nb of curve
{
  real[int,int] ff1(leman); // read  image and set to an rect. array
  int nx = ff1.n, ny=ff1.m; // grey value between 0 to 1 (dark)
  // build a Cartesian mesh such that the origne is qt the right place.
  mesh Th=square(nx-1,ny-1,[(nx-1)*(x),(ny-1)*(1-y)]);
   // warning  the numbering is of the vertices (x,y) is
   // given by $  i = x/nx + nx* y/ny $
  fespace Vh(Th,P1);
  Vh f1; f1[]=ff1; //  transforme array in finite element function.
  plot(f1,wait=1);
  nc=isoline(Th,f1,iso=0.25,close=1,Curves,beginend=be,smoothing=.1,ratio=0.5);
  verbosity=1;
}
// the longuest isoline
int ic0=be(0), ic1=be(1)-1;
plot([Curves(0,ic0:ic1),Curves(1,ic0:ic1)], wait=1);
int NC= Curves(2,ic1)/hsize;
real xl = Curves(0,ic0:ic1).max-5;
real yl = Curves(1,ic0:ic1).min+5;
border G(t=0,1) {  P=Curve(Curves,ic0,ic1,t);  label= 1 + (x>xl)*2 + (y<yl);}

plot(G(-NC),wait=1);
mesh Th=buildmesh(G(-NC));
plot(Th,wait=1);
real scale = sqrt(AreaLac/Th.area);
Th=movemesh(Th,[x*scale,y*scale]); // resize the  mesh to have the correct scale
cout << " Th.area = " << Th.area << " Km^2 " << " == " << AreaLac <<  "   Km^2 " << endl ;
plot(Th,wait=1,ps="o/leman.eps");
