1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
|
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");
|
