//Code from F. Hecht - FreeFem++ - LJLL - UPMC //Adapted by Y. Deleuze - LJLL , UPMC - SCCS , NTU - 2014 //yannick.deleuze @ ljll.math.upmc.fr load "msh3" load "medit" load "ppm2rnm" load "isoline" string taiwan="Taiwan2.pgm"; string taiwanelevation="Taiwanelevation.pgm"; real AreaTaiwan = 35883.; // $Km^2$ real hsize= 15; real[int,int] Curves(3,1); int[int] be(1); int nc;// nb of curve real[int,int] ff1(taiwan); // read image and set to an rect. array real[int,int] ffe(taiwanelevation); // read image and set to an rect. array int nx = ff1.n, ny=ff1.m; // grey value beetween 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. Vh fe; fe=0; fe[]=ffe; // transforme array in finite element function. // plot(fe,wait=1,value=1); fe=1-fe; // plot(fe,wait=1,value=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