C and MATLAB codes: Phase-field method for 3D volume reconstruction

%%%%%The following MATLAB code is used to preprocess the scattered points of a teapot.

clear;clc;close all; %.stl files store the points positions of a teapot %.stl files should be downloaded from other sites fe11=stlread(‘hofmeek11.stl’);fe12=stlread(‘hofmeek12.stl’); fe13=stlread(‘hofmeek13.stl’);fe14=stlread(‘hofmeek14.stl’); fe15=stlread(‘hofmeek15.stl’);fe16=stlread(‘hofmeek16.stl’); fe21=stlread(‘hofmeek21.stl’);fe22=stlread(‘hofmeek22.stl’); fe23=stlread(‘hofmeek23.stl’);fe24=stlread(‘hofmeek24.stl’); fe25=stlread(‘hofmeek25.stl’);fe26=stlread(‘hofmeek26.stl’); points=[fe11.vertices;fe12.vertices;fe13.vertices;fe14.vertices;fe15.vertices;fe16.vertices; fe21.vertices;fe22.vertices;fe23.vertices;fe24.vertices;fe25.vertices;fe26.vertices]; figure(1);centerx=sum(points(:,1))/size(points,1); centery=sum(points(:,2))/size(points,1);centerz=sum(points(:,3))/size(points,1); points(:,1)=points(:,1)-centerx;points(:,2)=points(:,2)-centery;points(:,3)=points(:,3)-centerz; maxx=max(max(abs(points)));points=points/(maxx+2.5);points=(points+1)/2;points2=points; clear points;points=points2(1:2:end,:);clear points2;nx=50;theta=linspace(0,1,nx); for i=1:10 x=0.5+(0.10-i0.01)cos(theta2pi); y=0.5+(0.10-i0.01)sin(theta2pi); points2=cat(1,points,([x;y;0.0*x+0.5])’);points=points2; end plot3(points(:,1),points(:,2),points(:,3),’k.’);axis image fid = fopen(‘fun_data.m’, ‘wt’); %store the new points positions for i=1:size(points,1) fprintf(fid, ‘%10.7f\n%10.7f\n%10.7f\n’, points(i,1),points(i,2),points(i,3)); end fclose(fid);fid = fopen(‘nt_points.m’, ‘wt’); %store the number of points fprintf(fid, ‘%10.7f\n’, size(points,1));fclose(fid);axis image

/** The following C code is used to calculate phase-file function. ***/

#include #include #include #include #include #include "util.h" #include "complex.h" #define NR_END 1 int nx, ny, nz, it,count, n_level, c_relax; float gam,lam, dt, h, h2, xleft, xright, yleft, yright, zleft, zright,pi, ***muu, ***mF, ***mFd, ***omFd, ***ct, ***sc, ***ooc, Sb, ***sorc, MM; void print_vol(float ***aa,int yu, int yd, int xl, int xr,int zl,int zr, int num); int main() {extern int nx, ny, nz, it,count, n_level,c_relax; extern float gam, dt, h, h2, xleft, xright, yleft, yright, zleft, zright,pi, ***muu, ***mF, ***mFd, ***omFd, ***ct, ***sc, ***ooc, Sb, ***sorc, MM; int i, j, k, max_it,ns; float ***oc, ***nc, ***gf, cahn, alpha; float elapsed,start,T,end; FILE *my; start = clock(); pi = 4.0*atan(1.0); nx = gnx; ny = gny;nz = gnz; xleft=0.0, xright= 2.0; yleft=0.0, yright= xright*ny/(1.0*nx); zleft=0.0, zright= xright*nz/(1.0*nx); c_relax = 5; count = 1; n_level=(int)(log2(ny)+0.001); h = (xright-xleft) / (1.0*nx); h2 = h*h; gam = 12.0*h/(2.0*sqrt(2)*atanh(0.9)); cahn = pow(gam, 2); dt=2.5e-5; T=64*dt; max_it = (int)(T/dt+0.5); ns =(int)( max_it/8+0.5); MM = 1.5; Sb = 3.0*(3.0*MM*MM-1.0); printf("nx = %d , ny = %d , nz = %d\n", nx, ny, nz); printf("dt = %f\n",dt); printf("max_it = %d\n", max_it); printf("ns = %d\n", ns); printf("gam = %f\n", gam); printf("h = %f\n", h); printf("n_level = %d\n", n_level); my = fopen("data/remarks.m","w"); fprintf(my, "\n\n"); fprintf(my, "nx = %d\n", nx); fprintf(my, "ny = %d\n", ny); fprintf(my, "nz = %d\n", nz); fprintf(my, "dt = %f\n", dt); fprintf(my, "ns = %d\n", ns); fprintf(my, "gam = %f\n", gam); fprintf(my, "h = %f\n", h); fprintf(my, "n_level = %d\n\n", n_level); fclose(my); mF = cube(0, nx+1, 0, ny+1, 0, nz+1); mFd = cube(0, nx+1, 0, ny+1, 0, nz+1); omFd = cube(0, nx+1, 0, ny+1, 0, nz+1);muu = cube(0, nx+1, 0, ny+1, 0, nz+1); sorc = cube(0, nx+1, 0, ny+1, 0, nz+1); ooc = cube(0, nx+1, 0, ny+1, 0, nz+1); gf = cube(0, nx+1, 0, ny+1, 0, nz+1); oc = cube(0, nx+1, 0, ny+1, 0, nz+1); nc = cube(0, nx+1, 0, ny+1, 0, nz+1);ct = cube(1, nx, 1, ny, 1, nz); sc = cube(1, nx, 1, ny, 1, nz); initialization(oc,gf); cube_copy(nc, oc, 1, nx, 1, ny, 1, nz); cube_copy(ooc, oc, 1, nx, 1, ny, 1, nz); augmen_phi(oc, nx,ny,nz); print_data3(oc,count); count++; augmen_phi(gf, nx,ny,nz); for (it=1; it<= max_it; it++){ ijkloop{ if ( fabs(oc[i][j][k]) < MM ) mFd[i][j][k] = pow(oc[i][j][k],3) - oc[i][j][k]; else mFd[i][j][k] = (3.0*MM*MM-1.0)*oc[i][j][k] - 2.0*pow(MM,3); if ( fabs(ooc[i][j][k]) < MM ) omFd[i][j][k] = pow(ooc[i][j][k],3) - ooc[i][j][k]; else omFd[i][j][k] = (3.0*MM*MM-1.0)*ooc[i][j][k] - 2.0*pow(MM,3);} PF3d(oc,nc,gf);augmen_phi(nc, nx,ny,nz);augmen_phi(oc, nx,ny,nz); cube_copy(ooc, oc, 1, nx, 1, ny, 1, nz);cube_copy(oc, nc, 1, nx, 1, ny, 1, nz); printf("iteration= %d \n", it); if ((it % ns==0)) {print_data3(oc,count);count++;printf("\n counts= %d \n",count); end = clock();elapsed = ((float) (end - start)) / CLOCKS_PER_SEC; printf("Time elapsed %f\n",elapsed); my = fopen("data/remarks.m","a"); fprintf(my, "\n counts = %d CPU Time= %f\n", count-1, elapsed); fclose(my); }} end = clock(); elapsed = ((float) (end - start)) / CLOCKS_PER_SEC; my = fopen("data/remarks.m","a"); fclose(my);return 0;} void initialization(float ***phi,float ***gf) { extern float h;float gam2; double valuexyz; ILE *fp;int i, j, k,ik,npo,ia,ja,ka,msize; float x, y, z,d,length,max_mun=-1.0;gam2 =3.0*h/(2.0*sqrt(2)*atanh(0.9));length=2.0*h; fp = fopen(".../nt_points.m","r"); /*input scattered points*/ fscanf(fp, "%lf", &valuexyz); npo = (int)(valuexyz+0.01);fclose(fp); fp = fopen("D:/make_initial/fun_data.m","r"); ijkloop{ phi[i][j][k]=1.0;} msize=(int)(gam2/h*6+0.1); for (ik=1;ik<=npo;ik++){ fscanf(fp, "%lf", &valuexyz);x=valuexyz;fscanf(fp, "%lf", &valuexyz); y=valuexyz;fscanf(fp, "%lf", &valuexyz);z=valuexyz; ia = (int)(x/h+0.499999)-1;ja = (int)(y/h+0.499999)-1;ka = (int)(z/h+0.499999)-1; for (i=ia-msize+1; i<=ia+msize; i++) { for (j=ja-msize+1; j<=ja+msize; j++) { for (k=ka-msize+1; k<=ka+msize; k++) { d=sqrt(pow((i-0.5)*h-x,2)+pow((j-0.5)*h-y,2)+pow((k-0.5)*h-z,2)); if (d<phi[i][j][k]) {phi[i][j][k]=d;} } } } } fclose(fp); ijkloop{ phi[i][j][k] =-tanh((phi[i][j][k]-sqrt(2)*gam2-2.0*h)/(sqrt(2)*gam2)); gf[i][j][k]= 1.0-phi[i][j][k]*phi[i][j][k] + 1.0e-10; } } void PF3d(float ***c_old, float ***c_new,float ***gf) { extern int nx, ny, nz;extern float gam, h, Sb, dt, ***ooc, ***mFd, ***omFd, ***sorc, ***intc; int i,j,k, max_it_CH = 100, it_mg = 1;float tol = 1.0e-6, resid = 1.0; float ***sor1;sor1 = cube(1, nx, 1, ny, 1, nz); cube_copy(sor1, c_new, 1, nx, 1, ny, 1, nz);augmen_phi(c_old, nx,ny,nz); ijkloop{ sorc[i][j][k] = c_old[i][j][k]/dt - gf[i][j][k]*( (1.5/pow(gam,2))*( mFd[i][j][k]

• Sbc_old[i][j][k] ) - (0.5/pow(gam,2))( mFd[i][j][k] - Sb*ooc[i][j][k] )
• 0.5(c_old[i+1][j][k]+c_old[i-1][j][k]+c_old[i][j+1][k]+c_old[i][j-1][k] +c_old[i][j][k+1]+c_old[i][j][k-1]-6.0c_old[i][j][k])/(h*h)
• 0.5Sbc_old[i][j][k]/pow(gam,2) ); } while (it_mg <= max_it_CH && resid > tol) { vcycle(c_new, sorc, gf,nx, ny, nz, 1);resid = error3(sor1, c_new, nx, ny, nz); cube_copy(sor1, c_new, 1, nx, 1, ny, 1, nz);it_mg++; } free_cube(sor1, 1, nx, 1, ny, 1, nz); return; } void vcycle(float **uf_new, float **su, float **gf, int nxf, int nyf, int nzf, int ilevel) { relax(uf_new, su,gf, ilevel, nxf, nyf, nzf); } void relax(float **c, float **sc,float **gf, int ilevel, int nxt, int nyt, int nzt) { extern float dt, Sb, gam, xleft,xright; extern int c_relax; int i, j, k,ik; float ht2, x_fac, y_fac, z_fac, a[1], f[1],xs,xl,ys,yl,zs,zl; ht2 = pow((xright - xleft)/(float)nxt,2); xs=1.0;xl=1.0;ys=1.0;yl=1.0;zs=1.0;zl=1.0; for (ik=0;ik<c_relax;ik++) ijkloopt { f[0] = 0.0; if (i == 1){ x_fac =xl; f[0] += c[i+1][j][k]xl;} else if(i == nxt){x_fac =xs;f[0] +=c[i-1][j][k]xs; } else{x_fac =xs+xl;f[0] += c[i-1][j][k]xs+c[i+1][j][k]xl;} if (j == 1){y_fac =yl;f[0] += c[i][j+1][k]yl;} else if(j == nyt){y_fac =ys;f[0] +=c[i][j-1][k]ys;} else{ y_fac =ys+yl;f[0] += c[i][j-1][k]ys+c[i][j+1][k]yl; } if (k == 1){ z_fac =zl; f[0] += c[i][j][k+1]zl; } else if(k == nzt){z_fac =zs;f[0] +=c[i][j][k-1]zs;} else{z_fac =zs+zl;f[0] += c[i][j][k-1]zs+c[i][j][k+1]zl;} a[0] = 1.0/dt + 0.5Sbgf[i][j][k]/pow(gam,2) + 0.5(x_fac + y_fac + z_fac)/ht2gf[i][j][k]; f[0] = 0.5f[0]/ht2gf[i][j][k]+sc[i][j][k]; c[i][j][k] = f[0]/a[0];} } void laplace_ch_gf(float **a, float **lap_a,float **gf, int nxt, int nyt, int nzt) { extern float xright;int i, j, k; float dxt2, dadx_L, dadx_R, dady_B, dady_F, dadz_D, dadz_U; dxt2 = pow(xright/(float)nxt,2); ijkloopt {if (i > 1) dadx_L = a[i][j][k] - a[i-1][j][k]; else dadx_L = 0.0; if (i < nxt) dadx_R = a[i+1][j][k] - a[i][j][k]; else dadx_R = 0.0; if (j > 1) dady_B = a[i][j][k] - a[i][j-1][k]; else dady_B = 0.0; if (j < nyt) dady_F = a[i][j+1][k] - a[i][j][k]; else dady_F = 0.0; if (k > 1) dadz_D = a[i][j][k] - a[i][j][k-1]; else dadz_D = 0.0; if (k < nzt) dadz_U = a[i][j][k+1] - a[i][j][k]; else dadz_U = 0.0; lap_a[i][j][k] =gf[i][j][k](dadx_R-dadx_L + dady_F- dady_B +dadz_U- dadz_D) / dxt2;}} float error3(float **c_old, float **c_new, int nxt, int nyt, int nzt) { extern int nx, ny, nz;int i, j, k;float value; value = 0.0; ijkloop{value += pow(c_old[i][j][k]-c_new[i][j][k],2);} value = sqrt(value/(1.0nxnynz));return value;} void print_data3( float **c, int count) { char buffer[20];FILE fc;int i,j,k; sprintf(buffer,”data/phi%d.m”,count);fc = fopen(buffer,”w”); ijkloop {fprintf(fc, “ %10.7f\n”, c[i][j][k]);}fclose(fc);return;} float error(float **c_old, float **c_new, int nxt, int nyt, int nzt) { extern int nx, ny, nz;int i, j, k;float value;value = 0.0; ijkloop{value += pow(c_old[i][j][k]-c_new[i][j][k],2);} value = sqrt(value/(1.0nxnynz));return value;} float **cube(int xl, int xr, int yl, int yr, int zl, int zr) { int i,j,nrow=xr-xl+1, ncol=yr-yl+1, ndep=zr-zl+1; float **t;t = (float **) malloc(((nrow+1)sizeof(float *))); t += 1;t -= xl; t[xl]=(float **) malloc(((nrowncol+1)sizeof(float *))); t[xl] += 1;t[xl] -= yl; t[xl][yl] = (float *) malloc(((nrowncolndep+1)sizeof(float ))); t[xl][yl] += 1;t[xl][yl] -= zl; for(j=yl+1; j<=yr; j++) t[xl][j]=t[xl][j-1]+ndep; for(i=xl+1; i<=xr; i++){ t[i] = t[i-1]+ncol;t[i][yl] = t[i-1][yl]+ncolndep; for(j=yl+1; j<=yr; j++)t[i][j] = t[i][j-1]+ndep;}return t;} void free_cube(float **t, int xl, int xr, int yl, int yr, int zl, int zr) { free((char ) (t[xl][yl]+zl-1));free((char *) (t[xl]+yl-1));free((char *) (t+xl-1));} void cube_copy(float **a, float **b,int xl, int xr, int yl, int yr, int zl, int zr) { int i, j, k;for (i=xl; i<=xr; i++)for (j=yl; j<=yr; j++)for (k=zl; k<=zr; k++){ a[i][j][k]=b[i][j][k];}return;} void augmen_phi(float **c,int nxt,int nyt,int nzt) { int i, j, k;for (j=1; j<=nyt; j++) for (k=1; k<=nzt; k++) { c[0][j][k] = c[1][j][k];c[nxt+1][j][k] = c[nxt][j][k]; } for (i=0; i<=nxt+1; i++)for (k=1; k<=nzt; k++) { c[i][0][k] = c[i][1][k];c[i][nyt+1][k] = c[i][nyt][k];} for (i=0; i<=nxt+1; i++) for (j=0; j<=nyt+1; j++) { c[i][j][0] = c[i][j][1];c[i][j][nzt+1] = c[i][j][nzt]; }}

/** The following util.h file is associated with C code. ***/

#define gnx 1282 #define gny 128 #define gnz 128 #define iloop for(i=1;i<=gnx;i++) #define jloop for(j=1;j<=gny;j++) #define kloop for(k=1;k<=gnz;k++) #define ijkloop iloop jloop kloop #define kjiloop kloop jloop iloop #define iloopp for(i=0;i<=gnx;i++) #define jloopp for(j=0;j<=gny;j++) #define kloopp for(k=0;k<=gnz;k++) #define ijkloopp iloopp jloopp kloopp #define kjiloopp kloopp jloopp iloopp #define iloopt for(i=1;i<=nxt;i++) #define jloopt for(j=1;j<=nyt;j++) #define kloopt for(k=1;k<=nzt;k++) #define ijkloopt iloopt jloopt kloopt #define kjiloopt kloopt jloopt iloopt float **cube(int xl, int xr, int yl, int yr, int zl, int zr); void free_cube(float **t, int xl, int xr, int yl, int yr, int zl, int zr); void print_data(float **phi); void augmen_phi(float ***c,int nxt,int nyt,int nzt);

/** The following complex.h file is associated with C code. **/

void initialization(float **phi,float **gf); void PF3d(float **c_old, float **c_new, float **gf); void vcycle(float **uf_new, float **su,float **gf, int nxf, int nyf, int nzf, int ilevel); void relax(float **c, float **sc,float **gf, int ilevel, int nxt, int nyt, int nzt); float error(float **c_old, float **c_new, int nxt, int nyt, int nzt); void laplace_ch_gf(float **a, float **lap_a,float **gf, int nxt, int nyt, int nzt); float error3(float **c_old, float **c_new, int nxt, int nyt, int nzt); void print_data3( float ***c, int count);

%The following MATLAB code is used to show the 3D profile of a teapot.

clear; nx=256; ny=128; nz=128; h=2/nx; x = linspace(0,2,nx); y = linspace(0,1,ny); z = linspace(0,1,nz); [xx,yy,zz] = meshgrid(y,x,z); for ip= 1 ss=sprintf(‘…/phi%d.m’,ip); u=load(ss); S(1:nx,1:ny,1:nz)=0; for i=1:nx for k=1:nz for j=1:ny aa=u(nynz(i-1)+nz*(j-1)+k); S(i,j,k)=aa; end end end figure(ip); p=patch(isosurface(xx,yy,zz,S,0.0)); set(p,’FaceColor’,[0.9 0.9 0.9],’EdgeColor’,’none’); daspect([1 1 1]); axis([0 1 0 2 0 1]); view(-132,25); camlight;lighting phong; axis image; end