//Code from F. Hecht - FreeFem++ - LJLL - UPMC //Adapted by Y. Deleuze - LJLL , UPMC - SCCS , NTU - 2013 load "ppm2rnm" load "isoline" string taiwan="taiwan.pgm"; real AreaTaiwan = 35883.; // $Km^2$ real hsize= 50; 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 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. 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