pgm2msh

// Original version : Chesapeake-mesh.edp by F.Hecht

//  pgm2msh.edp

//   Setouchi.pgm -> Th-Setouchi.msh

load "ppm2rnm"

real hmin=10; // 最小メッシュサイズ [狭い領域に合わせて設定。元データのサイズの 1/50 程度]

string FigureName="Setouchi"; // 今回は 600x300程度

real[int,int] ff1(FigureName+".pgm"); // 画像データを2次元配列に読み込み[原点は左上] 

int nx // 画像の2値化 : f2 = 1 (海), 0 (陸)

real f2threshold = 0.2;

f2= f1 < f2threshold; // f2[] = 1 (f1[] < f2threshold) , 0 (f1[]>=f2threshold)

cout << "plot f2 = f1 < " << f2threshold <<endl;

plot(f1,cmm="元画像",value=1,wait=1,dim=2,fill=1);

plot(f2,cmm="2値化画像",value=1,wait=1,dim=2,fill=1);

cout << "hmin=" << hmin << endl;

Th=adaptmesh(Th,f2,nbvx=1000000,hmin=hmin); // 最小メッシュがhminになるようにメッシュ最適化

cout << "Step 1 ( adaptmesh(hmin=" << hmin << ") )" << ", hmin=" << hmin << endl;

plot(Th,cmm="step 1",wait=1);

Th=adaptmesh(Th,f2,nbvx=1000000,hmin=hmin); // もう一度最適化して，メッシュをならす

cout << "Step 2 ( adaptmesh(hmin=" << hmin <<") )" << ", hmin=" << hmin << endl;

plot(Th,cmm="step 2",wait=1);

real vtrunmesh Th("Th-Setouchi.msh");c=0.2;

f2=f2; // 新しいメッシュに f2 を補間 [必要]

Th=trunc(Th,f2>vtrunc); // 海辺の細かいメッシュを除去 (f2[]>vtrunc)

cout << "Truncate (f2 > "+vtrunc+")" << endl;

plot(Th,cmm="[final] truncate f2>"+vtrunc, wait=1);

savemesh(Th,"Th-"+FigureName+".msh"); // メッシュデータ出力 = ff1.n, ny=ff1.m; 

mesh Th=square(nx-1,ny-1,[(nx-1)*(x),(ny-1)*(1-y)]); // 原点(0,0)を左下に変更

cout << "nx=" << nx << ", ny=" << ny <<endl;

// warning  the numbering is of the vertices (x,y) is 

// given by $  i = x/nx + nx* y/ny $

fespace Vh(Th,P1);

Vh f1,f2;

f1[]=ff1; // 画像データ配列を有限要素にコピー

// 画像の2値化 : f2 = 1 (海), 0 (陸)

real f2threshold = 0.2;

f2= f1 < f2threshold; // f2[] = 1 (f1[] < f2threshold) , 0 (f1[]>=f2threshold)

cout << "plot f2 = f1 < " << f2threshold <<endl;

plot(f1,cmm="元画像",value=1,wait=1,dim=2,fill=1);

plot(f2,cmm="2値化画像",value=1,wait=1,dim=2,fill=1);

cout << "hmin=" << hmin << endl;

Th=adaptmesh(Th,f2,nbvx=1000000,hmin=hmin); // 最小メッシュがhminになるようにメッシュ最適化

cout << "Step 1 ( adaptmesh(hmin=" << hmin << ") )" << ", hmin=" << hmin << endl;

plot(Th,cmm="step 1",wait=1);

Th=adaptmesh(Th,f2,nbvx=1000000,hmin=hmin); // もう一度最適化して，メッシュをならす

cout << "Step 2 ( adaptmesh(hmin=" << hmin <<") )" << ", hmin=" << hmin << endl;

plot(Th,cmm="step 2",wait=1);

real vtrunc=0.2;

f2=f2; // 新しいメッシュに f2 を補間 [必要]

Th=trunc(Th,f2>vtrunc); // 海辺の細かいメッシュを除去 (f2[]>vtrunc)

cout << "Truncate (f2 > "+vtrunc+")" << endl;

plot(Th,cmm="[final] truncate f2>"+vtrunc, wait=1);

savemesh(Th,"Th-"+FigureName+".msh"); // メッシュデータ出力