diff options
Diffstat (limited to 'Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h')
-rw-r--r-- | Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h | 2127 |
1 files changed, 2127 insertions, 0 deletions
diff --git a/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h new file mode 100644 index 00000000000..27a5434fe41 --- /dev/null +++ b/Master/texmf-dist/scripts/ketcindy/ketlibC/ketcommon.h @@ -0,0 +1,2127 @@ +//20180308 norm added + +int length1(double gdata[]){ + return floor(gdata[0]+0.5); +} + +int length2(double gdata[][2]){ + return floor(gdata[0][0]+0.5); +} + +int length3(double gdata[][3]){ + return floor(gdata[0][0]+0.5); +} + +int length4(double gdata[][4]){ + return floor(gdata[0][0]+0.5); +} + +int length5(double gdata[][5]){ + return floor(gdata[0][0]+0.5); +} + +int length6(double gdata[][6]){ + return floor(gdata[0][0]+0.5); +} + +int length2i(int gdata[][2]){ + return gdata[0][0]; +} + +void copyd(int from, int upto, double q[], double p[]){ + int i; + for(i=from; i<=upto; i++){ + p[i]=q[i]; + } +} + +void copyi(int from, int upto, int q[], int p[]){ + int i; + for(i=from; i<=upto; i++){ + p[i]=q[i]; + } +} + +int push2(double x, double y, double ptL[][2]){ + int nall=length2(ptL); + nall++; + ptL[nall][0]=x; ptL[nall][1]=y; + ptL[0][0]=nall; + return nall; +} + +int push3(double x, double y, double z, double ptL[][3]){ + int nall=length3(ptL); + nall++; + ptL[nall][0]=x; ptL[nall][1]=y; ptL[nall][2]=z; + ptL[0][0]=nall; + return nall; +} + +int push4(double x, double y, double z, double w, double ptL[][4]){ + int nall=length4(ptL); + nall++; + ptL[nall][0]=x; ptL[nall][1]=y; + ptL[nall][2]=z; ptL[nall][3]=w; + ptL[0][0]=nall; + return nall; +} + +int push5(double x, double y, double z, double w, double u, double ptL[][5]){ + int nall=length5(ptL); + nall++; + ptL[nall][0]=x; ptL[nall][1]=y; + ptL[nall][2]=z; ptL[nall][3]=w; ptL[nall][4]=u; + ptL[0][0]=nall; + return nall; +} + +int push6(double x, double y, double z, double w, double u, double v,double ptL[][6]){ + int nall=ptL[0][0]; + nall++; + ptL[nall][0]=x; ptL[nall][1]=y;ptL[nall][2]=z; + ptL[nall][3]=w; ptL[nall][4]=u; ptL[nall][5]=v; + ptL[0][0]=nall; + return nall; +} + +void pull3(int num, double ptL[][3],double pt[3]){ + pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; pt[2]=ptL[num][2]; +} + +void pull2(int num, double ptL[][2],double pt[2]){ + pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; +} + +void pull4(int num, double ptL[][4],double pt[4]){ + pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; + pt[2]=ptL[num][2]; pt[3]=ptL[num][3]; +} + +void pull5(int num, double ptL[][5],double pt[5]){ + pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; + pt[2]=ptL[num][2]; pt[3]=ptL[num][3]; pt[4]=ptL[num][4]; +} + +void pull6(int num, double ptL[][6],double pt[6]){ + pt[0]=ptL[num][0]; pt[1]=ptL[num][1]; pt[2]=ptL[num][2]; + pt[3]=ptL[num][3]; pt[4]=ptL[num][4]; pt[5]=ptL[num][5]; +} + +int add2(double gdata[][2], double x, double y){ + int nall; + nall=push2(x, y, gdata); + return nall; +} + +int add3(double gdata[][3], double x, double y, double z){ + int nall; + nall=push3(x, y, z,gdata); + return nall; +} + +int addptd3(double gdata[][3], double pt[3]){ + int nall; + nall=push3(pt[0],pt[1],pt[2],gdata); + return nall; +} + +int addptd2(double gdata[][2], double pt[2]){ + int nall; + nall=push2(pt[0],pt[1],gdata); + return nall; +} + +int addptd4(double gdata[][4], double pt[4]){ + int nall; + nall=push4(pt[0],pt[1],pt[2],pt[3],gdata); + return nall; +} + +int addptd5(double gdata[][5], double pt[5]){ + int nall; + nall=push5(pt[0],pt[1],pt[2],pt[3],pt[4],gdata); + return nall; +} + +int addptd6(double gdata[][6], double pt[6]){ + int nall; + nall=push6(pt[0],pt[1],pt[2],pt[3],pt[4],pt[5],gdata); + return nall; +} + +int appendd2(int level, int from,int upto, double gdata[][2],double mat[][2]){ + int j, np, nall=length2(mat); + np=upto-from+1; + if(np<0){np=0;} + if(np>0){ + if(level>0){ + if(nall>0){ + nall=add2(mat,Inf, level); + } + } + for(j=1; j<=np; j++){ + mat[nall+j][0]=gdata[from+j-1][0]; + mat[nall+j][1]=gdata[from+j-1][1]; + } + } + mat[0][0]=nall+np; + return nall+np; +} + +int appendd3(int level, int from,int upto, double gdata[][3],double mat[][3]){ + int j, np, nall=length3(mat); + np=upto-from+1; + if(np<0){np=0;} + if(np>0){ + if(level>0){ + if(nall>0){ + nall=add3(mat,Inf, level, 0); + } + } + for(j=1; j<=np; j++){ + mat[nall+j][0]=gdata[from+j-1][0]; + mat[nall+j][1]=gdata[from+j-1][1]; + mat[nall+j][2]=gdata[from+j-1][2]; + } + } + mat[0][0]=nall+np; + return nall+np; +} + +int prependd2(int from,int upto, double gdata[][2],double mat[][2]){ + int j,np, nall=length2(mat); + if(from<=upto){ + np=upto-from+1; + for(j=nall; j>=1; j--){ + mat[np+j][0]=mat[j][0]; + mat[np+j][1]=mat[j][1]; + } + for(j=1; j<=np; j++){ + mat[j][0]=gdata[from+j-1][0]; + mat[j][1]=gdata[from+j-1][1]; + } + }else{ + np=from-upto+1; + for(j=nall; j>=1; j--){ + mat[np+j][0]=mat[j][0]; + mat[np+j][1]=mat[j][1]; + } + for(j=1; j<=np; j++){ + mat[nall+j][0]=gdata[from-j+1][0]; + mat[nall+j][1]=gdata[from-j+1][1]; + } + } + mat[0][0]=nall+np; + return nall+np; +} + +int appendd1(double numvec[], double num){ + int n; + n=floor(numvec[0]+0.5); + numvec[n+1]=num; + numvec[0]=n+1; + return n+1; +} + +int appendi1(int numvec[],int num){ + int n; + n=numvec[0]; + numvec[n+1]=num; + numvec[0]=n+1; + return n+1; +} + +int prependi1(int num, int numvec[]){ + int i,n; + n=numvec[0]; + for(i=n; i>=1; i--){ + numvec[i+1]=numvec[i]; + } + numvec[1]=num; + numvec[0]=n+1; + return n+1; +} + +int removei2(int nrmv,int mat[][2]){ + int j,nall=mat[0][0]; + if(nall>0){ + for(j=nrmv; j<=nall-1; j++){ + mat[j][0]=mat[j+1][0]; + mat[j][1]=mat[j+1][1]; + } + mat[0][0]=nall-1; + return nall-1; + } + else{ + mat[0][0]=0; + return 0; + } +} + +int dataindexd1(double out[],int din[][2]){ + int n1=1, j, ctr=0, nall=floor(out[0]+0.5); + if(nall==0){ + din[0][0]=0; + return 0; + } + for(j=1; j<=nall; j++){ + if(out[j]>Inf-1){ + ctr++; + din[ctr][0]=n1; + din[ctr][1]=j-1; + n1=j+1; + } + } + ctr++; + din[ctr][0]=n1; + din[ctr][1]=nall; + din[0][0]=ctr; + return ctr; +} + +int dataindexd2(int level,double out[][2],int din[][2]){ + int n1=1, j, ctr=0, nall=length2(out); + if(nall==0){ + din[0][0]=0; + return 0; + } + for(j=1; j<=nall; j++){ + if(out[j][0]>Inf-1){ + if((level==0) || (out[j][1]==level)){ + ctr++; + din[ctr][0]=n1; + din[ctr][1]=j-1; + n1=j+1; + } + } + } + ctr++; + din[ctr][0]=n1; + din[ctr][1]=nall; + din[0][0]=ctr; + return ctr; +} + +int dataindexd3(int level, double out[][3],int din[][2]){ + int n1=1, j, ctr=0, nall=length3(out); + if(nall==0){ + din[0][0]=0; + return 0; + } + for(j=1; j<=nall; j++){ + if(out[j][0]>Inf-1){ + if((level==0) || (out[j][1]==level)){ + ctr++; + din[ctr][0]=n1; + din[ctr][1]=j-1; + n1=j+1; + } + } + } + ctr++; + din[ctr][0]=n1; + din[ctr][1]=nall; + din[0][0]=ctr; + return ctr; +} + +int dataindexd4(int level, double out[][4],int din[][2]){ + int n1=1, j, ctr=0, nall=length4(out); + if(nall==0){ + din[0][0]=0; + return 0; + } + for(j=1; j<=nall; j++){ + if(out[j][0]>Inf-1){ + if((level==0) || (out[j][1]==level)){ + ctr++; + din[ctr][0]=n1; + din[ctr][1]=j-1; + n1=j+1; + } + } + } + ctr++; + din[ctr][0]=n1; + din[ctr][1]=nall; + din[0][0]=ctr; + return ctr; +} + +int dataindexi1(int out[],int din[][2]){ + int n1=1, j, ctr=0, nall=out[0], ninfp=0; + if(nall==0){ + din[0][0]=0; din[0][1]=-1; + return 0; + } + for(j=1; j<=nall; j++){ + if(out[j]==Infint){ + ctr++; + if(j==ninfp+1){ + din[ctr][0]=0; din[ctr][1]=-1; + } + else{ + din[ctr][0]=n1; din[ctr][1]=j-1; + } + ninfp=j; + n1=j+1; + } + } + ctr++; + if(out[nall]==Infint){ + din[ctr][0]=0; din[ctr][1]=0; + } + else{ + din[ctr][0]=n1; din[ctr][1]=nall; + } + din[0][0]=ctr; + return ctr; +} + +void dispvec(int n,double dt[]){ + int i; + for(i=0; i<n; i++){ + printf("%f ",dt[i]); + } + printf("\n"); +} + +void dispveci(int n,int dt[]){ + int i; + for(i=0; i<n; i++){ + printf("%d ",dt[i]); + } + printf("\n"); +} + +void dispmatd2(int from, int upto, double mat[][2]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %7.5f %7.5f\n",i,mat[i][0],mat[i][1]); + } +} + +void dispmatd3(int from, int upto, double mat[][3]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %7.5f %7.5f %7.5f\n",i,mat[i][0],mat[i][1],mat[i][2]); + } +} + +void dispmatd4(int from, int upto, double mat[][4]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]); + printf("%7.5f %7.5f\n",mat[i][2],mat[i][3]); + } +} + +void dispmatd5(int from, int upto, double mat[][5]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]); + printf("%7.5f %7.5f %7.5f\n",mat[i][2],mat[i][3],mat[i][4]); + } +} + +void dispmatd6(int from, int upto, double mat[][6]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %7.5f %7.5f %7.5f ",i,mat[i][0],mat[i][1],mat[i][2]); + printf("%7.5f %7.5f %7.5f\n",mat[i][3],mat[i][4],mat[i][5]); + } +} + +void dispmatd2all(double mat[][2]){ + int i; + if(mat[0][0]==0){printf("%s\n","no data");return;} + for(i=1; i<mat[0][0]+1; i++){ + printf("%3d %7.5f %7.5f\n",i,mat[i][0],mat[i][1]); + } +} + +void dispmatd3all(double mat[][3]){ + int i; + if(mat[0][0]==0){printf("%s\n","no data");return;} + for(i=1; i<mat[0][0]+1; i++){ + printf("%3d %7.5f %7.5f %7.5f\n",i,mat[i][0],mat[i][1],mat[i][2]); + } +} + +void dispmatd4all(double mat[][4]){ + int i; + if(mat[0][0]==0){printf("%s\n","no data");return;} + for(i=1; i<mat[0][0]+1; i++){ + printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]); + printf("%7.5f %7.5f\n",mat[i][2],mat[i][3]); + } +} + +void dispmatd5all(double mat[][5]){ + int i; + if(mat[0][0]==0){printf("%s\n","no data");return;} + for(i=1; i<mat[0][0]+1; i++){ + printf("%3d %7.5f %7.5f ",i,mat[i][0],mat[i][1]); + printf("%7.5f %7.5f %7.5f\n",mat[i][2],mat[i][3],mat[i][4]); + } +} + +void dispmatd6all(double mat[][6]){ + int i; + if(mat[0][0]==0){printf("%s\n","no data");return;} + for(i=1; i<mat[0][0]+1; i++){ + printf("%3d %7.5f %7.5f %7.5f ",i,mat[i][0],mat[i][1],mat[i][2]); + printf("%7.5f %7.5f %7.5f\n",mat[i][3],mat[i][4],mat[i][5]); + } +} + +void dispmati2(int from, int upto, int mat[][2]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %d %d\n",i,mat[i][0],mat[i][1]); + } +} + +void dispmati1(int from, int upto, int mat[]){ + int i; + for(i=from; i<upto+1; i++){ + printf("%3d %d\n",i,mat[i]); + } +} + +void readdata3(const char *fname, double out[][3]){ + double data[DsizeL][4]; + int i, j, ndata, nall; + FILE *fp; + fp=fopen(fname,"r"); + if(fp==NULL){ + printf("file not found\n"); + return; + } + ndata=0; + i=1; + while( ! feof(fp) && i<DsizeL){ + for(j=0;j<4;j++){ + fscanf(fp,"%lf",&data[i][j]); + } + i++; + ndata++; + } + ndata--; + fclose(fp); + out[0][0]=0; nall=0; + for(i=1; i<=ndata; i++){ + nall=add3(out, data[i][0], data[i][1], data[i][2]); + } + return; +} + +void readdata2(const char *fname, double out[][2]){ + double data[DsizeL][3]; + int i, j, ndata, nall; + FILE *fp; + fp=fopen(fname,"r"); + if(fp==NULL){ + printf("file not found\n"); + return; + } + ndata=0; + i=1; + while( ! feof(fp) && i<DsizeL){ + for(j=0;j<3;j++){ + fscanf(fp,"%lf",&data[i][j]); + } + i++; + ndata++; + } + ndata--; + fclose(fp); + out[0][0]=0; nall=0; + for(i=1; i<=ndata; i++){ + nall=add2(out, data[i][0], data[i][1]); + } + return; +} + +int output2(const char *var, const char *fname, int level, double out[][2]){ + int nall=length2(out), ndin, din[DsizeM][2],i,j; + double tmpd[2]; + FILE *fp; + fp=fopen(fname,"w"); + if (fp == NULL) { + printf("cannot open\n"); + return -1; + } + ndin=dataindexd2(level, out,din); + fprintf(fp,"%s//\n",var); + for(i=1; i<=ndin; i++){ + fprintf(fp,"start//\n"); + fprintf(fp,"["); + for(j=din[i][0]; j<=din[i][1]; j++){ + pull2(j,out,tmpd); + fprintf(fp,"[ %7.5f, %7.5f]",tmpd[0],tmpd[1]); + if(j<din[i][1]){ + fprintf(fp,","); + } + } + fprintf(fp,"]//\n"); + if(i<ndin){ + fprintf(fp,"end//\n"); + } + else{ + fprintf(fp,"end////\n"); + } + } + fclose(fp); + return 0; +} + +void simplesort(double number[]){ + int i,j, total=number[0]; + double tmp; + for (i=1; i<=total; i++) { + for (j=i+1; j<=total; j++) { + if (number[i] > number[j]) { + tmp = number[i]; + number[i] = number[j]; + number[j] = tmp; + } + } + } +} + +int memberi(int a, int list[]){ + int i, n=list[0]; + for(i=1; i<=n; i++){ + if(a==list[i]){ + return 1; + } + } + return 0; +} + +double dotprod(int dim, double p[], double q[]){ + double ans; + ans=p[0]*q[0]+p[1]*q[1]; + if(dim>2){ + ans=ans+p[2]*q[2]; + } + return ans; +} + +double norm(int dim, double p[]){ + double ans; + ans=dotprod(dim,p,p); + return sqrt(ans); +} + +double dist(int dim, double p[], double q[]){ + double tmp; + tmp=pow(q[0]-p[0],2)+pow(q[1]-p[1],2); + if(dim>2){ + tmp=tmp+pow(q[2]-p[2],2); + } + return sqrt(tmp); +} + +void crossprod(int dim,double a[3],double b[3], double out[3]){ + if(dim==3){ + out[0]=a[1]*b[2]-a[2]*b[1]; + out[1]=a[2]*b[0]-a[0]*b[2]; + out[2]=a[0]*b[1]-a[1]*b[0]; + }else{ + out[0]=a[0]*b[1]-a[1]*b[0]; + out[1]=Inf; + } +} + +void reflectpoint(double p1[2],double regard[][2],double p2[2]){ + double pa[2],pb[2],u,v,a,b,u2,v2; + pull2(1,regard,pa); + if(length2(regard)==1){ + p2[0]=2*pa[0]-p1[0]; p2[1]=2*pa[1]-p1[1]; + }else{ + pull2(2,regard,pb); + u=pb[0]-pa[0]; u2=pow(u,2.0); + v=pb[1]-pa[1]; v2=pow(v,2.0); + a=pa[0]; b=pa[1]; + p2[0]=(u2-v2)/(u2+v2)*p1[0]+2*u*v/(u2+v2)*p1[1]; + p2[0]=p2[0]-2*v*(u*b-v*a)/(u2+v2); + p2[1]=2*u*v/(u2+v2)*p1[0]-(u2-v2)/(u2+v2)*p1[1]; + p2[1]=p2[1]+2*u*(u*b-v*a)/(u2+v2); + } +} + +void pointoncurve(double t, double gdata[][2], double eps, double pt[]){ + double pa[2],pb[2], s; + int n, nall; + n=floor(t+eps); + s=fmax(t-n,0); + nall=length2(gdata); + if(n>=nall){ + pull2(nall,gdata,pt); +// pt[0]=gdata[nall][0]; pt[1]=gdata[nall][1]; + } + else{ + pull2(n,gdata,pa); + pull2(n+1,gdata,pb); + pt[0]=(1-s)*pa[0]+s*pb[0]; + pt[1]=(1-s)*pa[1]+s*pb[1]; + } +} + +int partcrv(double a,double b,double pkL[][2], double pL[][2]){ + int i, is, ie, nall; + double eps=pow(10,-3.2); + is=ceil(a); + ie=floor(b); + nall=0; + if(a<is-eps){ + nall++; + pL[nall][0]=(is-a)*pkL[is-1][0]+(1-is+a)*pkL[is][0]; + pL[nall][1]=(is-a)*pkL[is-1][1]+(1-is+a)*pkL[is][1]; + } + for(i=is; i<=ie; i++){ + nall++; + pL[nall][0]=pkL[i][0]; + pL[nall][1]=pkL[i][1]; + } + if(b>ie+eps){ + nall++; + pL[nall][0]=(1-b+ie)*pkL[ie][0]+(b-ie)*pkL[ie+1][0]; + pL[nall][1]=(1-b+ie)*pkL[ie][1]+(b-ie)*pkL[ie+1][1]; + } + pL[0][0]=nall; + return nall; +} + +int partcrv3(double a,double b,double pkL[][3], double pL[][3]){ + int i, is, ie, nall; + double eps=pow(10,-3.2); + is=ceil(a); + ie=floor(b); + nall=0; + if(a<is-eps){ + nall++; + pL[nall][0]=(is-a)*pkL[is-1][0]+(1-is+a)*pkL[is][0]; + pL[nall][1]=(is-a)*pkL[is-1][1]+(1-is+a)*pkL[is][1]; + pL[nall][2]=(is-a)*pkL[is-1][2]+(1-is+a)*pkL[is][2]; + } + for(i=is; i<=ie; i++){ + nall++; + pL[nall][0]=pkL[i][0]; + pL[nall][1]=pkL[i][1]; + pL[nall][2]=pkL[i][2]; + } + if(b>ie+eps){ + nall++; + pL[nall][0]=(1-b+ie)*pkL[ie][0]+(b-ie)*pkL[ie+1][0]; + pL[nall][1]=(1-b+ie)*pkL[ie][1]+(b-ie)*pkL[ie+1][1]; + pL[nall][2]=(1-b+ie)*pkL[ie][2]+(b-ie)*pkL[ie+1][2]; + } + pL[0][0]=nall; + return nall; +} + +int connectseg(double pdata[][2], double eps, double out[][2]){ + int din[DsizeM][2], nd, nq, i, j, flg, st, en, ntmp; + int nall=length2(pdata), nout=0; + double qd[DsizeM][2], ah[2], ao[2], p[2], q[2], tmp[2]; + out[0][0]=0; + nd=dataindexd2(0,pdata,din); + qd[0][0]=0; + nq=appendd2(0,din[1][0],din[1][1],pdata,qd); + nd=removei2(1,din); + while(nd>0){ + pull2(1,qd,ah); + pull2(nq,qd,ao); + flg=0; + for(j=1; j<=nd; j++){ + st=din[j][0]; + en=din[j][1]; + pull2(st,pdata,p); + pull2(en,pdata,q); + if(dist(2,p,ao)<eps){ + nq=appendd2(0,st+1,en,pdata,qd); + nd=removei2(j,din); + flg=1; + break; + } + if(dist(2,q,ao)<eps){ + nq=appendd2(0,en-1,st,pdata,qd); + nd=removei2(j,din); + flg=1; + break; + } + if(dist(2,p,ah)<eps){ + nq=prependd2(en,st+1,pdata,qd); + nd=removei2(j,din); + flg=1; + break; + } + if(dist(2,q,ah)<eps){ + nq=prependd2(st,en-1,pdata,qd); + nd=removei2(j,din); + flg=1; + break; + } + } + if(flg==0){ + nout=appendd2(2,1,nq,qd,out); + qd[0][0]=0; + nq=appendd2(0,din[1][0],din[1][1],pdata,qd); + nd=removei2(1,din); + } + } + if(nq>0){ + nout=appendd2(2,1,nq,qd,out); + } + return nout; +} + +void koutenseg(double a[2], double b[2], double c[2], double d[2], + double eps, double eps2, double out[4]){ + double v[2], epss,eps0=Eps, sv2, p1, p2, q1,q2,m1, m2, M1, M2, t; + double tmp1d2[2], tmp2d1, p[2], q[2]; + v[0]=b[0]-a[0]; v[1]=b[1]-a[1]; + sv2=pow(v[0],2.0)+pow(v[1],2.0); + if(sv2<pow(10,-6.0)){ + out[0]=Inf; + return; + } + p[0]=c[0]-a[0]; p[1]=c[1]-a[1]; + q[0]=d[0]-a[0]; q[1]=d[1]-a[1]; + epss=fmin(eps2, eps/sqrt(sv2)); + p1=(p[0]*v[0]+p[1]*v[1])/sv2; + p2=(-p[0]*v[1]+p[1]*v[0])/sv2; + q1=(q[0]*v[0]+q[1]*v[1])/sv2; + q2=(-q[0]*v[1]+q[1]*v[0])/sv2; + m1=-epss; M1=1+epss; + m2=-epss; M2=epss; + if((fmax(p1,q1)<m1) || (fmin(p1,q1)>M1)){ + out[0]=Inf; + return; + } + if((fmax(p2,q2)<m2) || (fmin(p2,q2)>M2)){ + out[0]=Inf; + return; + } + if(p2*q2<0){ + t=p1-(q1-p1)/(q2-p2)*p2; + if((t>m1)&& (t<M1)){ + if((t>-eps0) && (t<1+eps0)){ + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + tmp2d1=fmin(fmax(t,0),1); + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=tmp2d1; out[3]=0; + } + else{ + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + tmp2d1=fmin(fmax(t,0),1); + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=tmp2d1; out[3]=1; + } + return; + } + if((p1<m1)||(p1>M1)||(p2<m2)||(p2>M2)){ + if((q1<m1)||(q1>M1)||(q2<m2)||(q2>M2)){ + out[0]=Inf; + return; + } + t=fmin(fmax(q1,0),1); + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=1; + return; + } + t=fmin(fmax(p1,0),1); + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=1; + return; + } + if((p1>-eps0)&&(p1<1+eps0)&&(p2>-eps0)&&(p2<eps0)){ + t= p1; + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=0; + return; + } + if((q1>-eps0)&&(q1<1+eps0)&&(q2>-eps0)&&(q2<eps0)){ + t= q1; + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=0; + return; + } + if((p1<m1)||(p1>M1)||(p2<m2)||(p2>M2)){ + if((q1<m1)||(q1>M1)||(q2<m2)||(q2>M2)){ + out[0]=Inf; + return; + } + t=fmin(fmax(q1,0),1); + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=1; + return; + } + if((q1<m1)||(q1>M1)||(q2<m2)||(q2>M2)){ + t=fmin(fmax(p1,0),1); + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=1; + return; + } + if(fabs(p2)<fabs(q2)){ + t=fmin(fmax(p1,0),1); + } + else{ + t=fmin(fmax(q1,0),1); + } + tmp1d2[0]=a[0]+t*v[0]; tmp1d2[1]=a[1]+t*v[1]; + out[0]=tmp1d2[0]; out[1]=tmp1d2[1]; + out[2]=t; out[3]=1; + return; +} + +double paramoncurve(double p[2], int n, double ptL[][2]){ + double pa[2], pb[2], v[2], w[2], o[2], d2, s; + int nptL=length2(ptL); + if(n==nptL){ + return nptL; + } + else{ + pull2(n,ptL,pa); + pull2(n+1,ptL,pb); + v[0]=pb[0]-pa[0]; v[1]=pb[1]-pa[1]; + w[0]=p[0]-pa[0]; w[1]=p[1]-pa[1]; + o[0]=0; o[1]=0; + d2=pow(dist(2,v,o),2.0); + s=dotprod(2,v,w)/d2; + s=fmin(fmax(s,0),1); + return n+s; + } +} + +void nearestptpt(double pa[2], double pL[][2], double pt[4]){ + double eps0=pow(10,-6.0), p[2], pm[2], im, sm, s; + double ans[7], a1,a2, b1,b2, v1, v2, x1, x2, t, tmpd1; + int i, n, npL=length2(pL); + ans[0]=pa[0]; ans[1]=pa[1]; + ans[2]=1; + p[0]=pL[1][0]; p[1]=pL[1][1]; + ans[3]=p[0]; ans[4]=p[1]; + ans[5]=1; ans[6]=dist(2,pa,p); + pm[0]=pL[1][0]; pm[1]=pL[1][1]; + im=1; + sm=dist(2,pm,pa); + for(i=1; i<=npL-1; i++){ + a1=pL[i][0]; a2=pL[i][1]; + b1=pL[i+1][0]; b2=pL[i+1][1]; + v1=b1-a1; v2=b2-a2; + x1=pa[0]; x2=pa[1]; + tmpd1=pow(v2,2.0)+pow(v1,2.0); + if(fabs(tmpd1)<eps0){ + continue; + } + t=(-a2*v2-v1*a1+v1*x1+x2*v2)/tmpd1; + if(t<-eps0){ + p[0]=a1; p[1]=a2; + } + else if(t>1+eps0){ + p[0]=b1; p[1]=b2; + } + else{ + p[0]=a1+t*v1; p[1]=a2+t*v2; + } + s=dist(2,p,pa); + if(s<sm-eps0){ + tmpd1=paramoncurve(p,i,pL); + pm[0]=p[0]; pm[1]=p[1]; + im=tmpd1; sm=s; + } + } + if(sm<ans[6]){ + ans[0]=pa[0]; ans[1]=pa[1]; + ans[2]=n; + ans[3]=pm[0]; ans[4]=pm[1]; + ans[5]=im; ans[6]=sm; + } + pt[0]=ans[3]; pt[1]=ans[4]; + pt[2]=ans[5]; pt[3]=ans[6]; +} + +int intersectselfPp(double data1[][2],double eps,double eps1s,double kL[][4]){ + double a[2], b[2], tmpd2[2], tmpd4[4], tmpd5[5], p[2], q[2], t, s; + double eps0=pow(10,-5.0); + double tmpd1, tmp1d1, tmp2d1; + int i, j, k, n, is, ndata1,nkL, ntmp, ntmp1,flg, closed=0; + ndata1=length2(data1); + p[0]=data1[1][0]; p[1]=data1[1][1]; + q[0]=data1[ndata1][0]; q[1]=data1[ndata1][1]; + if(dist(2,p,q)<Eps){closed=1;} + kL[0][0]=0; nkL=0; + for(i=1; i<=ndata1; i++){ + ntmp=0; tmpd1=Inf; + pull2(i,data1,a); + nearestptpt(a, data1, tmpd4); + n=floor(tmpd4[2]+0.5); + if(n==i){continue;} + if(closed==1){ + if((i==1)&&(n==ndata1)){continue;} + if((i==ndata1)&&(n==1)){continue;} + } + ntmp1=floor(tmpd4[2]+Eps); + if((tmpd4[3]<eps1s)&&(ntmp1!=i)){ + if(ntmp<ntmp1){ + ntmp=ntmp1; + tmpd1=tmpd4[3]; + nkL=addptd4(kL,tmpd4); + } + else if(tmpd4[3]<tmpd1){ + tmpd1=tmpd4[3]; + kL[0][0]--; + nkL=addptd4(kL,tmpd4); + } + } + } + return nkL; +} + +/* +int intersectcrvsPptest(double data1[][2], double data2[][2], + double eps1,double eps2, double kL[][4]){ + double t, s, a[2], b[2], p[2], q[2]; + double sqeps=pow(10,-10.0), eps0=pow(10,-2.0); + double tmp1d2[2],tmp2d2[2],tmpd4[4],tmp1d4[4],tmpd5[5]; + double tmpd1,tmp1d1,tmp2d1, kL0[DsizeL][5], ds; + int i, j, k, n, is, nkL, nall, ndata1,ndata2,ntmp, ntmp1, flg; + ndata1=length2(data1); + ndata2=length2(data2); + if(ndata1==ndata2){ + tmp1d1=0; tmp2d1=0; + for(i=1; i<=ndata1; i++){ + tmp1d1=tmp1d1+pow(data1[i][0]-data2[i][0],2.0); + tmp1d1=tmp1d1+pow(data1[i][1]-data2[i][1],2.0); + tmp2d1=tmp2d1+pow(data1[i][0]-data2[ndata2-i+1][0],2.0); + tmp2d1=tmp2d1+pow(data1[i][1]-data2[ndata2-i+1][1],2.0); + } + tmp1d1=sqrt(tmp1d1); tmp2d1=sqrt(tmp2d1); + if((tmp1d1<Eps)||(tmp2d1<Eps)){ + kL[0][0]=0; nkL=0; + return nkL; + } + } + kL0[0][0]=0; nkL=0; + for(i=1; i<=ndata1; i++){ + ntmp=0; tmpd1=Inf; + pull2(i,data1,a); + nearestptpt(a, data2, tmpd4); + tmp1d1=floor(tmpd4[2]+Eps); //17.06.03 + a[0]=tmpd4[0]; a[1]=tmpd4[1]; + nearestptpt(a, data1, tmpd4); + tmpd5[0]=tmpd4[0]; tmpd5[1]=tmpd4[1]; + tmpd5[2]=tmpd4[2]; tmpd5[3]=tmp1d1; + tmpd5[4]=tmpd4[3]; + nkL=addptd5(kL0,tmpd5); + } + kL[0][0]=0; nall=0; + p[0]=Inf; p[1]=0; + for(i=1; i<=nkL; i++){ + if(kL0[i][4]<eps2){ + q[0]=kL0[i][0]; q[1]=kL0[i][1]; + if(dist(2,p,q)>eps2){ + nall=push4(kL0[i][0],kL0[i][1],kL0[i][2],kL0[i][3],kL); + p[0]=kL0[i][0]; p[1]=kL0[i][1]; + ds=kL0[i][4]; + continue; + } + if(kL0[i][4]<ds){ + kL[0][0]--; + nall=push4(kL0[i][0],kL0[i][1],kL0[i][2],kL0[i][3],kL); + p[0]=kL0[i][0]; p[1]=kL0[i][1]; + ds=kL0[i][4]; + } + } + } + return nall; +} +*/ + +int intersectcrvsPp(double data1[][2], double data2[][2], + double eps,double eps1s, double kL[][4]){ + double t, s, a[2], b[2], p[2], q[2]; + double sqeps=pow(10,-10.0), eps0=pow(10,-2.0); + double tmp1d2[2],tmp2d2[2],tmpd4[4],tmp1d4[4],tmpd5[5]; + double tmpd1,tmp1d1,tmp2d1; + int i, j, k, n, is, ndata1, ndata2, nkL,nkL1,nkL2,ntmp, ntmp1, flg; + ndata1=length2(data1); + ndata2=length2(data2); + if(ndata1==ndata2){ + tmp1d1=0; tmp2d1=0; + for(i=1; i<=ndata1; i++){ + tmp1d1=tmp1d1+pow(data1[i][0]-data2[i][0],2.0); + tmp1d1=tmp1d1+pow(data1[i][1]-data2[i][1],2.0); + tmp2d1=tmp2d1+pow(data1[i][0]-data2[ndata2-i+1][0],2.0); + tmp2d1=tmp2d1+pow(data1[i][1]-data2[ndata2-i+1][1],2.0); + } + tmp1d1=sqrt(tmp1d1); tmp2d1=sqrt(tmp2d1); + if((tmp1d1<Eps)||(tmp2d1<Eps)){ + kL[0][0]=0; nkL=0; + return nkL; + } + } + kL[0][0]=0; nkL=0; + for(i=1; i<=ndata1;i++){ + ntmp=0; tmpd1=Inf; + pull2(i,data1,a); + nearestptpt(a, data2, tmpd4); + tmp1d1=floor(tmpd4[2]+Eps); //17.06.03 + a[0]=tmpd4[0]; a[1]=tmpd4[1]; + nearestptpt(a, data1, tmpd4); + if(tmpd4[3]<eps1s){ + flg=0; + p[0]=tmpd4[0]; p[1]=tmpd4[1]; + for(j=1; j<=nkL; j++){ + q[0]=kL[j][0]; q[1]=kL[j][1]; + if(dist(2,p,q)<eps0){ + flg=1; + break; + } + } + if(flg==1){continue;} + ntmp1=floor(tmpd4[2]+Eps); + if(ntmp<ntmp1){ + ntmp=ntmp1; + tmpd1=tmpd4[3]; + tmp1d4[0]=tmpd4[0]; tmp1d4[1]=tmpd4[1]; + tmp1d4[2]=tmpd4[2]; tmp1d4[3]=tmp1d1; + nkL=addptd4(kL,tmp1d4); + } + else if(tmpd4[3]<tmpd1){ + tmpd1=tmpd4[3]; + kL[0][0]--; + tmp1d4[0]=tmpd4[0]; tmp1d4[1]=tmpd4[1]; + tmp1d4[2]=tmpd4[2]; tmp1d4[3]=tmp1d1; + nkL=addptd4(kL,tmp1d4); + } + } + } + return nkL; +} + +int intersectcrvsPpold(double data1[][2], double data2[][2], + double eps,double eps1s, double kL[][4]){ +// 17.05.22 + double kL1[DsizeL][5], kL2[DsizeL][5]; + double t, s, a[2], b[2], p[2], q[2]; + double sqeps=pow(10,-10.0), eps0=pow(10,-6.0); + double tmp1d2[2],tmp2d2[2],tmpd4[4],tmpd5[5]; + double tmpd1,tmp1d1,tmp2d1; + int i, j, k, n, is, ndata1, ndata2, nkL,nkL1,nkL2,ntmp,flg; + ndata1=length2(data1); + ndata2=length2(data2); + if(ndata1==ndata2){ + tmp1d1=0; tmp2d1=0; + for(i=1; i<=ndata1; i++){ + tmp1d1=tmp1d1+pow(data1[i][0]-data2[i][0],2.0); + tmp1d1=tmp1d1+pow(data1[i][1]-data2[i][1],2.0); + tmp2d1=tmp2d1+pow(data1[i][0]-data2[ndata2-i+1][0],2.0); + tmp2d1=tmp2d1+pow(data1[i][1]-data2[ndata2-i+1][1],2.0); + } + tmp1d1=sqrt(tmp1d1); tmp2d1=sqrt(tmp2d1); + if((tmp1d1<eps0)||(tmp2d1<eps0)){ + kL[0][0]=0; nkL=0; + return nkL; + } + } + kL1[0][0]=0; nkL1=0; + kL2[0][0]=0; nkL2=0; + for(i=1; i<=ndata1-1; i++){ + pull2(i,data1,a); + pull2(i+1,data1,b); + for(j=1; j<=ndata2-1; j++){ + pull2(j, data2, p); + pull2(j+1, data2, q); + koutenseg(a,b,p,q,eps,eps1s,tmpd4); + if(tmpd4[0]<Inf-1){ + ntmp=floor(tmpd4[3]+Eps); + if(ntmp==0){ + tmpd5[0]=tmpd4[0]; tmpd5[1]=tmpd4[1]; + tmpd5[2]=tmpd4[2]; tmpd5[3]=i; tmpd5[4]=j; + nkL1=addptd5(kL1,tmpd5); + } + else{ + tmpd5[0]=tmpd4[0]; tmpd5[1]=tmpd4[1]; + tmpd5[2]=tmpd4[2]; tmpd5[3]=i; tmpd5[4]=j; + nkL2=addptd5(kL2,tmpd5); + } + } + } + } + kL[0][0]=0; nkL=0; + if(nkL1>0){ + pull5(1,kL1,tmpd5); + p[0]=tmpd5[0]; p[1]=tmpd5[1]; + t=tmpd5[2]; + i=floor(tmpd5[3]+Eps); + j=floor(tmpd5[4]+Eps); + tmpd4[0]=p[0]; tmpd4[1]=p[1]; + tmpd4[2]=i+t; tmpd4[3]=j; + nkL=addptd4(kL,tmpd4); + } + for(n=2; n<=nkL1; n++){ + pull5(n,kL1,tmpd5); + p[0]=tmpd5[0]; p[1]=tmpd5[1]; + flg=0; + for(k=1; k<=nkL; k++){ + pull4(k,kL,tmpd4); + tmp1d2[0]=tmpd4[0]; tmp1d2[1]=tmpd4[1]; + if(pow(dist(2,p,tmp1d2),2.0)<sqeps){ + flg=1; + break; + } + } + if(flg==0){ + tmpd4[0]=p[0]; tmpd4[1]=p[1]; + t=tmpd5[2]; + i=floor(tmpd5[3]+Eps); + j=floor(tmpd5[4]+Eps); + tmpd4[2]=t+i; tmpd4[3]=j; + nkL=addptd4(kL,tmpd4); + } + } + for(n=1; n<=nkL2; n++){ + pull5(n,kL2,tmpd5); + p[0]=tmpd5[0]; p[1]=tmpd5[1]; + flg=0; + for(k=1; k<=nkL; k++){ + pull4(k,kL,tmpd4); + tmp1d2[0]=tmpd4[0]; tmp1d2[1]=tmpd4[1]; + if(pow(dist(2,p,tmp1d2),2.0)<sqeps){ + flg=1; + break; + } + } + if(flg==0){ + tmpd4[0]=p[0]; tmpd4[1]=p[1]; + t=tmpd5[2]; + i=floor(tmpd5[3]+Eps); + j=floor(tmpd5[4]+Eps); + tmpd4[2]=i+t; tmpd4[3]=j; + nkL=addptd4(kL,tmpd4); + } + } + return nkL; +} + +int dropnumlistcrv(double qd[][2], double eps1, int out[]){ + int i,j,k, nout,nall=length2(qd), nd, se, en, npd, nptL; + int din[DsizeM][2],ptL[DsizeL]; + double eps=pow(10.0,-6.0), pd[DsizeL][2], p[2], tmp2d[2]; + out[0]=0; nout=0; + nd=dataindexd2(2,qd,din); + for(j=1; j<=nd; j++){ + se=din[j][0]; en=din[j][1]; + pd[0][0]=0; npd=0; + npd=appendd2(0,se,en,qd,pd); + ptL[0]=0; nptL=0; + nptL=appendi1(ptL,1); + pull2(1,pd,p); + for(k=2; k<=npd-1; k++){ + pull2(k,pd,tmp2d); + if(dist(2,p,tmp2d)>eps1){ + nptL=appendi1(ptL,k); + pull2(k,pd,p); + } + } + pull2(npd-1,pd,p); + pull2(npd,pd,tmp2d); + if(dist(2,p,tmp2d)>eps){ // eps -> eps1 ? + nptL=appendi1(ptL,npd); + } + if(nptL==1){ptL[0]=0; nptL=0;} + if(nout>0){ + nout=appendi1(out,Infint); + } + for(i=1; i<=nptL; i++){ + nout=appendi1(out,ptL[i]); + } + } + out[0]=nout; + return nout; +} + +int increasept2(double ptL[][2], int nn, double out[][2]){ + int ii, jj, npt, nall; + double tmpd2[2]; + npt=length2(ptL); + out[0][0]=0; + for(ii=1; ii<=npt-1; ii++){ + for(jj=0; jj<nn; jj++){ + tmpd2[0]=ptL[ii][0]+1.0*jj/nn*(ptL[ii+1][0]-ptL[ii][0]); + tmpd2[1]=ptL[ii][1]+1.0*jj/nn*(ptL[ii+1][1]-ptL[ii][1]); + addptd2(out,tmpd2); + } + } + nall=appendd2(0,npt,npt,ptL,out); + return nall; +} + +int increasept3(double ptL[][3], int nn, double out[][3]){ + int ii, jj, npt, nall; + double tmpd3[3]; + npt=length3(ptL); + out[0][0]=0; + for(ii=1; ii<=npt-1; ii++){ + for(jj=0; jj<nn; jj++){ + tmpd3[0]=ptL[ii][0]+1.0*jj/nn*(ptL[ii+1][0]-ptL[ii][0]); + tmpd3[1]=ptL[ii][1]+1.0*jj/nn*(ptL[ii+1][1]-ptL[ii][1]); + tmpd3[2]=ptL[ii][2]+1.0*jj/nn*(ptL[ii+1][2]-ptL[ii][2]); + addptd3(out,tmpd3); + } + } + nall=appendd3(0,npt,npt,ptL,out); + return nall; +} + +void bezierpt(double t, double p0[], double p3[], double p1[], double p2[], + double pout[]){ + double p4[2], p5[2], p6[2], p7[2], p8[2]; + int i, flg3=1; + if(p2[0]>Inf-1){ + for(i=0; i<2; i++){ + p2[i]=p3[i]; + } + flg3=0; + } + for(i=0; i<2; i++){ + p4[i]=(1-t)*p0[i]+t*p1[i]; + } + for(i=0; i<2; i++){ + p5[i]=(1-t)*p1[i]+t*p2[i]; + } + for(i=0; i<2; i++){ + p6[i]=(1-t)*p2[i]+t*p3[i]; + } + for(i=0; i<2; i++){ + p7[i]=(1-t)*p4[i]+t*p5[i]; + } + if(flg3==0){ + for(i=0; i<2; i++){ + pout[i]=p7[i]; + } + } + else{ + for(i=0; i<2; i++){ + p8[i]=(1-t)*p5[i]+t*p6[i]; + } + for(i=0; i<2; i++){ + pout[i]=(1-t)*p7[i]+t*p8[i]; + } + } +} + +int bezier(double ptL[][2], double ctrL[][4], int num, double out[][2]){ + double p0[2], p3[2], p1[2], p2[2], pt[2], dt=1.0/num; + int i, j, n, nall; + out[0][0]=0; nall=0; + nall=appendd2(0,1,1,ptL,out); + n=length2(ptL); + for(i=1; i<=n-1; i++){ + pull2(i,ptL, p0); + pull2(i+1,ptL, p3); + p1[0]=ctrL[i][0]; p1[1]=ctrL[i][1]; + p2[0]=ctrL[i][2]; p2[1]=ctrL[i][3]; + for(j=1; j<=num; j++){ + bezierpt(j*dt,p0,p3,p1,p2,pt); + nall=push2(pt[0],pt[1], out); + } + } + return nall; +} + +int ospline(int level, double ptL[][2], int num, double out[][2]){ + double ctrlist[DsizeM][4], eps0=pow(10,-6.0), tmp, cc; + double p[2], q[2], p0[2], p3[2], p1[2], p2[2], pQ[2], pR[2]; + double ptlist[DsizeM][2], tmpmd2[DsizeM][2]; + int nd, i, j, npt, nall, nctr, flg, ndin, closed, din[DsizeM][2], ntmp; + out[0][0]=0; + if(num==1){ + npt=length2(ptL); + for(i=0; i<=npt; i++){ + out[i][0]=ptL[i][0]; out[i][1]=ptL[i][1]; + } + nall=length2(out); + return nall; + } + ndin=dataindexd2(level,ptL,din); + for(nd=1; nd<=ndin; nd++){ + ptlist[0][0]=0; + npt=appendd2(0,din[nd][0],din[nd][1],ptL,ptlist); + if(npt==2){ + appendd2(2, 1,npt,ptlist,out); + continue; + } + closed=0; + pull2(1,ptlist,p); + pull2(npt,ptlist,q); + if(dist(2,p,q)<eps0){ + closed=1; + } + ctrlist[0][0]=0; nctr=0; + for(i=1; i<=npt-1; i++){ + flg=0; + pull2(i, ptlist,p1); + pull2(i+1, ptlist,p2); + if((i==1)||(i==npt-1)){ + flg=1; + if(closed==0){ + pQ[0]=(p1[0]+2*p2[0])/3.0; + pQ[1]=(p1[1]+2*p2[1])/3.0; + pR[0]=(2*p1[0]+p2[0])/3.0; + pR[1]=(2*p1[1]+p2[1])/3.0; + nctr++; + ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; + ctrlist[nctr][2]=pR[0]; ctrlist[nctr][3]=pR[1]; + flg=2; + } + else{ + if(i==1){ + pull2(npt-1, ptlist,p0); + pull2(i+2, ptlist,p3); + } + else{ + pull2(i-1, ptlist, p0); + pull2(2, ptlist, p3); + } + } + } + if(flg<=1){ + if(flg==0){ + pull2(i-1, ptlist, p0); + pull2(i+2, ptlist, p3); + } + p[0]=p2[0]-p0[0]; p[1]=p2[1]-p0[1]; + q[0]=p3[0]-p1[0]; q[1]=p3[1]-p1[1]; + tmp=1+sqrt((1+dotprod(2,p,q)/dist(2,p0,p2)/dist(2,p1,p3))/2); + cc=4*dist(2,p1,p2)/3/(dist(2,p0,p2)+dist(2,p1,p3))/tmp; + pQ[0]=p1[0]+cc*(p2[0]-p0[0]); pQ[1]=p1[1]+cc*(p2[1]-p0[1]); + pR[0]=p2[0]+cc*(p1[0]-p3[0]); pR[1]=p2[1]+cc*(p1[1]-p3[1]); + nctr++; + ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; + ctrlist[nctr][2]=pR[0]; ctrlist[nctr][3]=pR[1]; + } + } + if(closed==0){ + p1[0]=ctrlist[2][0]; p1[1]=ctrlist[2][1]; + pull2(2, ptlist, p2); + pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]); + ctrlist[1][0]=pQ[0]; ctrlist[1][1]=pQ[1]; + ctrlist[1][2]=Inf; ctrlist[1][3]=0; + p1[0]=ctrlist[nctr-1][2]; p1[1]=ctrlist[nctr-1][3]; + pull2(nctr,ptlist,p2); + pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]); + ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; + ctrlist[nctr][2]=Inf; ctrlist[1][3]=0; + } + ctrlist[0][0]=nctr; + if(length2(ptlist)>2){ + ntmp=bezier(ptlist,ctrlist,num, tmpmd2); + nall=appendd2(2,1,ntmp,tmpmd2,out); + } + } + return nall; +} + +void bezierpt3(double t, double p0[], double p3[], double p1[], double p2[], + double pout[]){ + double p4[3], p5[3], p6[3], p7[3], p8[3]; + int i, flg3=1; + if(p2[0]>Inf-1){ + for(i=0; i<3; i++){ + p2[i]=p3[i]; + } + flg3=0; + } + for(i=0; i<3; i++){ + p4[i]=(1-t)*p0[i]+t*p1[i]; + } + for(i=0; i<3; i++){ + p5[i]=(1-t)*p1[i]+t*p2[i]; + } + for(i=0; i<3; i++){ + p6[i]=(1-t)*p2[i]+t*p3[i]; + } + for(i=0; i<3; i++){ + p7[i]=(1-t)*p4[i]+t*p5[i]; + } + if(flg3==0){ + for(i=0; i<3; i++){ + pout[i]=p7[i]; + } + } + else{ + for(i=0; i<3; i++){ + p8[i]=(1-t)*p5[i]+t*p6[i]; + } + for(i=0; i<3; i++){ + pout[i]=(1-t)*p7[i]+t*p8[i]; + } + } +} + +int bezier3(double ptL[][3], double ctrL[][6], int num, double out[][3]){ + double p0[3], p3[3], p1[3], p2[3], pt[3], dt=1.0/num; + int i, j, n, nall; + out[0][0]=0; nall=0; + nall=appendd3(0,1,1,ptL,out); + n=length3(ptL); + for(i=1; i<=n-1; i++){ + pull3(i,ptL, p0); + pull3(i+1,ptL, p3); + p1[0]=ctrL[i][0]; p1[1]=ctrL[i][1]; p1[2]=ctrL[i][2]; + p2[0]=ctrL[i][3]; p2[1]=ctrL[i][4]; p2[2]=ctrL[i][5]; + for(j=1; j<=num; j++){ + bezierpt3(j*dt,p0,p3,p1,p2,pt); + nall=push3(pt[0],pt[1],pt[2], out); + } + } + return nall; +} + +int ospline3(int level, double ptL[][3], int num, double out[][3]){ + double ctrlist[DsizeM][6], eps0=pow(10,-6.0), tmp, cc; + double p[3], q[3], p0[3], p3[3], p1[3], p2[3], pQ[3], pR[3]; + double ptlist[DsizeM][3], tmpmd2[DsizeM][3]; + int nd, i, j, n, npt, nall, nctr, closed, flg, ndin, din[DsizeM][2], ntmp; + out[0][0]=0; nall=0; + if(num==1){ + npt=length3(ptL); + for(i=0; i<=npt; i++){ + out[i][0]=ptL[i][0]; out[i][1]=ptL[i][1]; out[i][2]=ptL[i][2]; + } + nall=length3(out); + return nall; + } + ndin=dataindexd3(2,ptL,din); + for(nd=1; nd<=ndin; nd++){ + ptlist[0][0]=0; + npt=din[nd][1]-din[nd][0]; + if(npt<=3){ + n=appendd3(0, din[nd][0], din[nd][0], ptL, ptlist); + for(i=1; i<=npt-1; i++){ + pull3(din[nd][0]+i,ptL,p); + pull3(din[nd][0]+i+1,ptL,q); + for(j=1; j<=num; j++){ + p1[0]=p[0]+1.0*j/num*(q[0]-p[0]); + p1[1]=p[1]+1.0*j/num*(q[1]-p[1]); + p1[2]=p[2]+1.0*j/num*(q[2]-p[2]); + n=addptd3(ptlist,p1); + } + } + nall=appendd3(2, 1,n,ptlist,out); + continue; + } + appendd3(0,din[nd][0],din[nd][1],ptL,ptlist); + closed=0; + pull3(1,ptlist,p); + pull3(npt,ptlist,q); + if(dist(3,p,q)<eps0){closed=1;} + ctrlist[0][0]=0; nctr=0; + for(i=1; i<=npt-1; i++){ + flg=0; + pull3(i, ptlist,p1); + pull3(i+1, ptlist,p2); + if((i==1)||(i==npt-1)){ + flg=1; + if(closed==0){ + pQ[0]=(p1[0]+2*p2[0])/3.0; + pQ[1]=(p1[1]+2*p2[1])/3.0; + pQ[2]=(p1[2]+2*p2[2])/3.0; + pR[0]=(2*p1[0]+p2[0])/3.0; + pR[1]=(2*p1[1]+p2[1])/3.0; + pR[2]=(2*p1[2]+p2[2])/3.0; + nctr++; + ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; ctrlist[nctr][2]=pQ[2]; + ctrlist[nctr][3]=pR[0]; ctrlist[nctr][4]=pR[1]; ctrlist[nctr][5]=pR[2]; + flg=2; + } + else{ + if(i==1){ + pull3(npt-1, ptlist,p0); + pull3(i+2, ptlist,p3); + } + else{ + pull3(i-1, ptlist, p0); + pull3(2, ptlist, p3); + } + } + } + if(flg<=1){ + if(flg==0){ + pull3(i-1, ptlist, p0); + pull3(i+2, ptlist, p3); + } + p[0]=p2[0]-p0[0]; p[1]=p2[1]-p0[1]; p[2]=p2[2]-p0[2]; + q[0]=p3[0]-p1[0]; q[1]=p3[1]-p1[1]; q[2]=p3[2]-p1[2]; + tmp=1+sqrt((1+dotprod(3,p,q)/dist(3,p0,p2)/dist(3,p1,p3))/2); + cc=4*dist(3,p1,p2)/3/(dist(3,p0,p2)+dist(3,p1,p3))/tmp; + pQ[0]=p1[0]+cc*(p2[0]-p0[0]); pQ[1]=p1[1]+cc*(p2[1]-p0[1]); + pQ[2]=p1[2]+cc*(p2[2]-p0[2]); + pR[0]=p2[0]+cc*(p1[0]-p3[0]); pR[1]=p2[1]+cc*(p1[1]-p3[1]); + pR[2]=p2[2]+cc*(p1[2]-p3[2]); + nctr++; + ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; ctrlist[nctr][2]=pQ[2]; + ctrlist[nctr][3]=pR[0]; ctrlist[nctr][4]=pR[1]; ctrlist[nctr][5]=pR[2]; + } + } + if(closed==0){ + p1[0]=ctrlist[2][0]; p1[1]=ctrlist[2][1]; p1[2]=ctrlist[2][2]; + pull3(2, ptlist, p2); + pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]); + pQ[2]=p2[2]+(p2[2]-p1[2]); + ctrlist[1][0]=pQ[0]; ctrlist[1][1]=pQ[1]; ctrlist[1][2]=pQ[2]; + ctrlist[1][3]=Inf; ctrlist[1][4]=0; ctrlist[1][5]=0; + p1[0]=ctrlist[nctr-1][3]; p1[1]=ctrlist[nctr-1][4]; p1[2]=ctrlist[nctr-1][5]; + pull3(nctr,ptlist,p2); + pQ[0]=p2[0]+(p2[0]-p1[0]); pQ[1]=p2[1]+(p2[1]-p1[1]); + pQ[2]=p2[2]+(p2[2]-p1[2]); + ctrlist[nctr][0]=pQ[0]; ctrlist[nctr][1]=pQ[1]; ctrlist[nctr][2]=pQ[2]; + ctrlist[nctr][3]=Inf; ctrlist[1][4]=0; ctrlist[1][5]=0; + } + ctrlist[0][0]=nctr; + if(length3(ptlist)>2){ + ntmp=bezier3(ptlist,ctrlist,num, tmpmd2); + nall=appendd3(2,1,ntmp,tmpmd2,out); + } + } + return nall; +} + +void segplot(double x1,double y1,double x2,double y2,double seg[][2]){ + seg[0][0]=0; + push2(x1,y1,seg); + push2(x2,y2,seg); +} + +void addsegplot(double x1,double y1,double seg[][2]){ + push2(x1,y1,seg); +} + +void circledata(double cx,double cy,double r, int num,double out[][2]){ + double t,dt,pt[2]; + int i; + dt=2*M_PI/num; + out[0][0]=0; + for(i=0; i<=num; i++){ + push2(cx+r*cos(i*dt),cy+r*sin(i*dt),out); + } +} + +void intersectline(double p1[2],double v1[2], double p2[2],double v2[2],double out[4]){ + double tmp,d,dt,ds,t,s,tmp1[2],tmp2[2],tmp3[2],pt[2]; + tmp=dotprod(2,v1,v2); + tmp1[0]=tmp; tmp1[1]=pow(norm(2,v1),2.0); + tmp2[0]=-pow(norm(2,v2),2.0);tmp2[1]=-tmp; + pt[0]=p2[0]-p1[0];pt[1]=p2[1]-p1[1]; + tmp3[0]=dotprod(2,pt,v2);tmp3[1]=dotprod(2,pt,v1); + crossprod(2,tmp1,tmp2,pt); d=pt[0]; + crossprod(2,v1,v2,pt); tmp=fabs(pt[0]); + if(tmp>Eps){ + crossprod(2,tmp3,tmp2,pt); dt=pt[0]; + crossprod(2,tmp1,tmp3,pt); ds=pt[0]; + t=dt/d; + s=ds/d; + pt[0]=p1[0]+v1[0]*t; pt[1]=p1[1]+v1[1]*t; + out[0]=pt[0];out[1]=pt[1];out[2]=t; out[3]=s; + }else{ + pt[0]=p2[0]-p1[0];pt[1]=p2[1]-p1[1]; + crossprod(2,pt,v1,tmp1); + tmp=tmp1[0]/norm(2,v1); + out[0]=fabs(tmp); out[1]=Inf; + } +} + +void intersectseg(double seg1[][2],double seg2[][2], double out[5]){ + double p1[2],q1[2],p2[2],q2[2],pt[2],pt1[2],nv[2],pts[DsizeS][4]; + double q1p1[2],q2p2[2],p1q2[2],q1p2[2],p1p2[2],q1q2[2]; + double tmp[2],tmp1[2],tmp2[2],tmp3[2],tmp4[2]; + double tmpd1,x,y,t,s,t1,s1,distance,tmpd[4],tmp3d[7]; + int ctr=0, i; + pull2(1,seg1,p1); pull2(2,seg1,q1); + pull2(1,seg2,p2); pull2(2,seg2,q2); + q1p1[0]=q1[0]-p1[0];q1p1[1]=q1[1]-p1[1]; + q2p2[0]=q2[0]-p2[0];q2p2[1]=q2[1]-p2[1]; + p1q2[0]=p1[0]-q2[0];p1q2[1]=p1[1]-q2[1]; + q1p2[0]=q1[0]-p2[0];q1p2[1]=q1[1]-p2[1]; + p1p2[0]=p1[0]-p2[0];p1p2[1]=p1[1]-p2[1]; + p1q2[0]=p1[0]-q2[0];p1q2[1]=p1[1]-q2[1]; + q1q2[0]=q1[0]-q2[0];q1q2[1]=q1[1]-q2[1]; + if((norm(2,q1p1)<Eps)||(norm(2,q2p2)<Eps)){ + out[0]=-1; out[1]=Inf; + }else{ + intersectline(p1,q1p1,p2,q2p2,tmpd); + if(tmpd[1]<Inf-Eps){ + pt[0]=tmpd[0];pt[1]=tmpd[1];t=tmpd[2];s=tmpd[3]; + if((t*(t-1)<Eps)&&(s*(s-1)<Eps)){ + out[0]=0;out[1]=pt[0];out[2]=pt[1];out[3]=t;out[4]=s; + }else{ + if(t<0){t=0;}; if(t>1){t=1;}; + if(s<0){s=0;}; if(s>1){s=1;}; + tmp3d[ctr]=norm(2,p1p2);ctr++; + tmp3d[ctr]=norm(2,p1q2);ctr++; + tmp3d[ctr]=norm(2,q1p2);ctr++; + tmp3d[ctr]=norm(2,q1q2);ctr++; + tmp1[0]=q2p2[1];tmp1[1]=-q2p2[0]; + intersectline(p1,tmp1,p2,q2p2,tmpd); + if(tmpd[1]!=Inf){ + pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3]; + if(s1*(s1-1)<Eps){ + tmp3d[ctr]=dist(2,pt1,p1);ctr++; + } + } + intersectline(q1,tmp1,p2,q2p2,tmpd); + if(tmpd[1]!=Inf){ + pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3]; + if(s1*(s1-1)<Eps){ + tmp3d[ctr]=dist(2,pt1,q1);ctr++; + } + } + tmp1[0]=q1p1[1];tmp1[1]=-q1p1[0]; + intersectline(p2,tmp1,p1,q1p1,tmpd); + if(tmpd[1]!=Inf){ + pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3]; + if(s1*(s1-1)<Eps){ + tmp3d[ctr]=dist(2,pt1,p2);ctr++; + } + } + intersectline(q2,tmp1,p1,q1p1,tmpd); + if(tmpd[1]!=Inf){ + pt1[0]=tmpd[0];pt1[1]=tmpd[1];s1=tmpd[3]; + if(s1*(s1-1)<Eps){ + tmp3d[ctr]=dist(2,pt1,q2);ctr++; + } + } + tmpd1=tmp3d[0]; + for(i=1;i<ctr;i++){ + if(tmp3d[i]<tmpd1){tmpd1=tmp3d[i];}; + } + out[0]=tmpd1;out[1]=pt[0];out[2]=pt[1]; + out[3]=t;out[4]=s; + } + }else{ + distance=tmpd[0]; + nv[0]=q1p1[1]/norm(2,q1p1);nv[1]=-q1p1[0]/norm(2,q1p1); + pts[0][0]=0; + intersectline(p1,nv,p2,q2p2,tmpd); + s1=tmpd[3]; + if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){ + x=(1-s1)*p2[0]+s1*q2[0];y=(1-s1)*p2[1]+s1*q2[1]; + push4(x,y,0,s1,pts); + } + intersectline(q1,nv,p2,q2p2,tmpd); + s1=tmpd[3]; + if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){ + x=(1-s1)*p2[0]+s1*q2[0];y=(1-s1)*p2[1]+s1*q2[1]; + push4(x,y,1,s1,pts); + } + intersectline(p2,nv,p1,q1p1,tmpd); + s1=tmpd[3]; + if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){ + x=(1-s1)*p1[0]+s1*q1[0];y=(1-s1)*p1[1]+s1*q1[1]; + push4(x,y,s1,0,pts); + } + intersectline(q2,nv,p1,q1p1,tmpd); + if((tmpd[1]!=Inf)&&(s1*(s1-1)<Eps)){ + x=(1-s1)*p1[0]+s1*q1[0];y=(1-s1)*p1[1]+s1*q1[1]; + push4(x,y,s1,2,pts); + } + if(length4(pts)==0){ + tmpd1=norm(2,p1p2); + if(norm(2,p1q2)<tmpd1){tmpd1=norm(2,p1q2);} + if(norm(2,q1p2)<tmpd1){tmpd1=norm(2,q1p2);} + if(norm(2,q1q2)<tmpd1){tmpd1=norm(2,q1q2);} + out[0]=tmpd1;out[1]=Inf; + }else{ + if(distance>Eps1){ + out[0]=distance;out[1]=Inf; + }else{ + x=0; y=0; + for(i=1; i<=length4(pts); i++){ + x=x+pts[i][0]; + y=y+pts[i][1]; + } + tmp3[0]=x/length4(pts); tmp3[1]=y/length4(pts); + nearestptpt(tmp3,seg1,tmpd); + t=tmpd[2]; + nearestptpt(tmp3,seg2,tmpd); + s=tmpd[2]; + out[0]=distance; out[1]=tmp3[0]; out[2]=tmp3[1]; + out[3]=t;out[4]=s; + } + } + } + } +} + +int osplineseg(double ptlist[5][2],int num,double out[][2]){ + double p0[2],p1[2],p2[2],p3[2],p2p0[2],p3p1[2],p2p1[2]; + double tmp,cc,pq[2],pr[2],knots[3][2],ctrlist[2][4]; + pull2(1,ptlist,p0); pull2(2,ptlist,p1); + pull2(3,ptlist,p2); pull2(4,ptlist,p3); + p2p0[0]=p2[0]-p0[0];p2p0[1]=p2[1]-p0[1]; + p3p1[0]=p3[0]-p1[0];p3p1[1]=p3[1]-p1[1]; + p2p1[0]=p2[0]-p1[0];p2p1[1]=p2[1]-p1[1]; + tmp=norm(2,p2p0)*norm(2,p3p1); + tmp=1+sqrt((1+dotprod(2,p2p0,p3p1)/tmp)/2); + cc=4*norm(2,p2p1)/3/(norm(2,p2p0)+norm(2,p3p1))/tmp; + pq[0]=p1[0]+cc*p2p0[0]; pq[1]=p1[1]+cc*p2p0[1]; + pr[0]=p2[0]-cc*p3p1[0];pr[1]=p2[1]-cc*p3p1[1]; + ctrlist[0][0]=0; + push4(pq[0],pq[1],pr[0],pr[1],ctrlist); + knots[0][0]=0; + addptd2(knots,p1);addptd2(knots,p2); + bezier(knots,ctrlist,num,out); + return length2(out); +} + +void intersectpartseg(double crv1[][2],double crv2[][2],int ii, int jj,int num,double out[6]){ + double distance=10*Eps2,seg1[3][2],seg2[3][2],tmpd[2],tmp1d[2],tmp2d[2]; + double tmp,snang,dst,tmpd5[5],os1[DsizeS][2],os2[DsizeS][2]; + double p0[2],p1[2],p2[2],p3[2],ptlist[5][2],regard[3][2]; + double tmp1md[20][2],tmpd3[3],tmp3md3[20][3],tmp2md3[20][3]; + double tmpd4[4],Eps00=pow(10,-8.0); + double sx,sy; + int kk,ll,nn; + out[0]=Inf; + for(kk=1;kk<6;kk++){out[kk]=0;} + seg1[0][0]=0;seg2[0][0]=0; + appendd2(1,ii,ii+1,crv1,seg1); + appendd2(1,jj,jj+1,crv2,seg2); + pull2(2,seg1,tmp1d); + pull2(1,seg1,tmpd); + tmp1d[0]=tmp1d[0]-tmpd[0]; tmp1d[1]=tmp1d[1]-tmpd[1]; + pull2(2,seg2,tmp2d); + pull2(1,seg2,tmpd); + tmp2d[0]=tmp2d[0]-tmpd[0]; tmp2d[1]=tmp2d[1]-tmpd[1]; + crossprod(2,tmp1d,tmp2d,tmpd); + snang=fabs(tmpd[0])/(norm(2,tmp1d)*norm(2,tmp2d)); + intersectseg(seg1,seg2,tmpd5); + dst=tmpd5[0]; + if(dst<Eps){ + if(dst>-Eps){ + out[0]=tmpd5[1];out[1]=tmpd5[2];out[2]=ii+tmpd5[3]; + out[3]=jj+tmpd5[4];out[4]=dst;out[5]=snang; + } + }else{ + if(dst<Eps2){ + os1[0][0]=0;os2[0][0]=0; + pull2(1,seg1,tmp1d); + pull2(2,seg1,tmp2d); + tmpd[0]=tmp2d[0]-tmp1d[0];tmpd[1]=tmp2d[1]-tmp1d[1]; + if((length2(crv1)==2)||(norm(2,tmpd)>distance-Eps)){ + appendd2(1,1,2,seg1,os1); + }else{ + pull2(1,seg1,p1);pull2(2,seg1,p2); + if(ii==1){ + pull2(3,crv1,p3); + tmpd[0]=p2[0]-p1[0];tmpd[1]=p2[1]-p1[1]; + tmp1d[0]=(p1[0]+p2[0])/2; + tmp1d[1]=(p1[1]+p2[1])/2; + tmp2d[0]=tmp1d[0]+tmpd[1]; + tmp2d[1]=tmp1d[1]-tmpd[0]; + regard[0][0]=0; + addptd2(regard,tmp1d); addptd2(regard,tmp2d); + reflectpoint(p3,regard,p0); + }else{ + if(ii==length2(crv1)-1){ + pull2(ii-1,crv1,p0); + tmpd[0]=p2[0]-p1[0]; tmpd[1]=p2[1]-p1[1]; + tmp1d[0]=(p1[0]+p2[0])/2; + tmp1d[1]=(p1[1]+p2[1])/2; + tmp2d[0]=tmp1d[0]+tmpd[1]; + tmp2d[1]=tmp1d[1]-tmpd[0]; + regard[0][0]=0; + addptd2(regard,tmp1d); addptd2(regard,tmp2d); + reflectpoint(p0,regard,p3); + }else{ + pull2(ii-1,crv1,p0);pull2(ii+2,crv1,p3); + } + } + ptlist[0][0]=0; + addptd2(ptlist,p0);addptd2(ptlist,p1); + addptd2(ptlist,p2);addptd2(ptlist,p3); + osplineseg(ptlist,num,os1); + } + pull2(1,seg2,tmp1d); + pull2(2,seg2,tmp2d); + if((length2(crv2)==2)||(dist(2,tmp2d,tmp1d)>distance-Eps)){ + appendd2(1,1,2,seg2,os2); + }else{ + pull2(1,seg2,p1);pull2(2,seg2,p2); + if(jj==1){ + pull2(3,crv2,p3); + tmpd[0]=p2[0]-p1[0];tmpd[1]=p2[1]-p1[1]; + tmp1d[0]=(p1[0]+p2[0])/2; + tmp1d[1]=(p1[1]+p2[1])/2; + tmp2d[0]=tmp1d[0]+tmpd[1]; + tmp2d[1]=tmp1d[1]-tmpd[0]; + regard[0][0]=0; + addptd2(regard,tmp1d); addptd2(regard,tmp2d); + reflectpoint(p3,regard,p0); + }else{ + if(jj==length2(crv2)-1){ + pull2(jj-1,crv2,p0); + tmpd[0]=p2[0]-p1[0]; tmpd[1]=p2[1]-p1[1]; + tmp1d[0]=(p1[0]+p2[0])/2; + tmp1d[1]=(p1[1]+p2[1])/2; + tmp2d[0]=tmp1d[0]+tmpd[1]; + tmp2d[1]=tmp1d[1]-tmpd[0]; + regard[0][0]=0; + addptd2(regard,tmp1d); addptd2(regard,tmp2d); + reflectpoint(p0,regard,p3); + }else{ + pull2(jj-1,crv2,p0);pull2(jj+2,crv2,p3); + } + } + ptlist[0][0]=0; + addptd2(ptlist,p0);addptd2(ptlist,p1); + addptd2(ptlist,p2);addptd2(ptlist,p3); + osplineseg(ptlist,num,os2); + } + tmp2md3[0][0]=0; + for(kk=1;kk<length2(os1);kk++){ + for(ll=1;ll<length2(os2);ll++){ + seg1[0][0]=0; seg2[0][0]=0; + pull2(kk,os1,tmpd);addptd2(seg1,tmpd); + pull2(kk+1,os1,tmpd);addptd2(seg1,tmpd); + pull2(ll,os2,tmpd);addptd2(seg2,tmpd); + pull2(ll+1,os2,tmpd);addptd2(seg2,tmpd); + intersectseg(seg1,seg2,tmpd5); + if((tmpd5[0]<Eps1)&&(tmpd5[1]<Inf-Eps)){ + if(tmpd5[0]<dst+Eps00){ + dst=tmpd5[0]; + tmp3md3[0][0]=0; + for(nn=1;nn<=length3(tmp2md3);nn++){ + pull3(nn,tmp2md3,tmpd3); + if(tmpd3[0]<dst+Eps00){ + addptd3(tmp3md3,tmpd3); + } + } + tmp2md3[0][0]=0; + for(nn=1;nn<=length3(tmp3md3);nn++){ + pull3(nn,tmp3md3,tmpd3); + addptd3(tmp2md3,tmpd3); + } + push3(dst,tmpd5[1],tmpd5[2],tmp2md3); + } + } + } + } + if(length3(tmp2md3)>0){ + tmp1md[0][0]=0; + sx=0; sy=0; + for(nn=1;nn<=length3(tmp2md3);nn++){ + pull3(nn,tmp2md3,tmpd3); + push2(tmpd3[1],tmpd3[2],tmp1md); + sx=sx+tmpd3[1]; + sy=sy+tmpd3[2]; + } + sx=sx/length3(tmp2md3); + sy=sy/length3(tmp2md3); + p0[0]=sx; p0[1]=sy; + pull2(ii,crv1,p1); pull2(ii+1,crv1,p2); + tmp1d[0]=p2[0]-p1[0];tmp1d[1]=p2[1]-p1[1]; + tmpd[0]=tmp1d[1];tmpd[1]=-tmp1d[0]; + intersectline(p0,tmpd,p1,tmp1d,tmpd4); + tmp=tmpd4[3]; + if(tmp<0){tmp=0;} + if(tmp>1){tmp=1;} + out[0]=sx;out[1]=sy;out[2]=ii+tmp; + pull2(jj,crv2,p1); pull2(jj+1,crv2,p2); + tmp1d[0]=p2[0]-p1[0];tmp1d[1]=p2[1]-p1[1]; + tmpd[0]=tmp1d[1];tmpd[1]=-tmp1d[0]; + intersectline(p0,tmpd,p1,tmp1d,tmpd4); + tmp=tmpd4[3]; + if(tmp<0){tmp=0;} + if(tmp>1){tmp=1;} + out[3]=jj+tmp;out[4]=dst;out[5]=snang; + } + } + } +} + +void collectsameseg(double ptdL[][6],double rL[][6]){ + double Eps00=pow(10,-8.0),tmp1d[2],tmp2d[2],tmp1d6[6],tmp2d6[6]; + int ii,jj,kk,nr,ctr,numL[40],diffL[40]; + double s1,e1,s2,e2,dst,tmp1,tmp2,tmp1md6[40][6]; + rL[0][0]=0; + if(length6(ptdL)==0){ + return; + } + for(ii=2;ii<=length6(ptdL);ii++){ + pull6(ii,ptdL,tmp1d6);addptd6(rL,tmp1d6); + } + tmp1md6[0][0]=0; + pull6(1,ptdL,tmp1d6); + addptd6(tmp1md6,tmp1d6); + kk=tmp1d6[2]; + if(tmp1d6[2]<kk+Eps00){ + s1=kk-1-Eps00;e1=s1+2+2*Eps00; + }else{ + s1=kk-Eps00; e1=s1+1+2*Eps00; + } + kk=tmp1d6[3]; + if(tmp1d6[3]<kk+Eps00){ + s2=kk-1-Eps00;e2=s2+2+2*Eps00; + }else{ + s2=kk-Eps00; e2=s2+1+2*Eps00; + } + ctr=1;numL[0]=0; + for(ii=1;ii<=length6(rL);ii++){ + pull6(ii,rL,tmp1d6); + tmp1=tmp1d6[2]; tmp2=tmp1d6[3]; + if((tmp1>s1)&&(tmp1<e1)&&(tmp2>s2)&&(tmp2<e2)){ + addptd6(tmp1md6,tmp1d6); + numL[ctr]=ii;ctr++; + } + } + dst=100; + for(ii=1;ii<=length6(tmp1md6);ii++){ + pull6(ii,tmp1md6,tmp1d6); + if(tmp1d6[4]<dst){ + dst=tmp1d6[4]; + } + } + ptdL[0][0]=0; + for(ii=1;ii<=length6(tmp1md6);ii++){ + pull6(ii,tmp1md6,tmp1d6); + if(tmp1d6[4]<dst+Eps00){ + addptd6(ptdL,tmp1d6); + } + } + for(ii=0;ii<=length6(rL);ii++){ + diffL[ii]=0; + } + for(ii=1;ii<=ctr;ii++){ + kk=numL[ii]; + diffL[kk]=1; + } + nr=length6(rL); + rL[0][0]=0; + for(ii=1;ii<=nr;ii++){ + if(diffL[ii]==0){ + pull6(ii,rL,tmp1d6);addptd6(rL,tmp1d6); + } + } +} + +void intersectcurvesPp(double crv1s[][2],double crv2s[][2],int num,double out[][6]){ + double Eps00=pow(10,-8.0),Dist=10*Eps2, crv1[DsizeL][2],crv2[DsizeL][2]; + double tmp1md6[DsizeS][6],tmp2md6[DsizeS][6]; + double tmpd[2],tmp1d[2],tmp2d[2],tmpd6[6],tmp1d6[6]; + double sx,sy,dst,ans[4]; + int ii,jj,kk,nall,self,js,je; + crv1[0][0]=0; + pull2(1,crv1s,tmpd);addptd2(crv1,tmpd); + for(ii=2;ii<=length2(crv1s);ii++){ + pull2(length2(crv1),crv1,tmp1d); + pull2(ii,crv1s,tmp2d); + if(dist(2,tmp1d,tmp2d)>Eps){ + addptd2(crv1,tmp2d); + } + } + crv2[0][0]=0; + pull2(1,crv2s,tmpd); + addptd2(crv2,tmpd); + for(ii=2;ii<=length2(crv2s);ii++){ + pull2(length2(crv2),crv2,tmp1d); + pull2(ii,crv2s,tmp2d); + if(dist(2,tmp1d,tmp2d)>Eps){ + addptd2(crv2,tmp2d); + } + } + if(length2(crv1)!=length2(crv2)){ + self=0; + }else{ + self=1; + for(ii=1;ii<=length2(crv1);ii++){ + pull2(ii,crv1,tmp1d);pull2(ii,crv2,tmp2d); + if(dist(2,tmp1d,tmp2d)>0){ + self=0; + break; + } + } + } + tmp1md6[0][0]=0; + for(ii=1;ii<=length2(crv1)-1;ii++){ + if(self==0){ + js=1; je=length2(crv2)-1; + }else{ + js=ii+2; je=length2(crv2)-1; + } + for(jj=js;jj<=je;jj++){ + intersectpartseg(crv1,crv2,ii,jj,num,tmpd6); + if(tmpd6[0]<Inf-Eps){ + tmpd[0]=tmpd6[0];tmpd[1]=tmpd6[1]; + if(length6(tmp1md6)==0){ + addptd6(tmp1md6,tmpd6); + }else{ + pull6(length6(tmp1md6),tmp1md6,tmp1d6); + tmp1d[0]=tmp1d6[0];tmp1d[1]=tmp1d6[1]; + if(dist(2,tmp1d,tmpd)>Eps1){ + addptd6(tmp1md6,tmpd6); + } + } + if(self==1){ + push6(tmpd6[0],tmpd6[1],tmpd6[3],tmpd6[2],tmpd6[4],tmpd6[5],tmp1md6); + } + } + } + } + out[0][0]=0; + tmp2md6[0][0]=0; + while(length6(tmp1md6)>0){ + collectsameseg(tmp1md6,tmp2md6); + dst=100; + for(ii=1;ii<=length6(tmp1md6);ii++){ + pull6(ii,tmp1md6,tmpd6); + if(tmpd6[4]<dst){ + dst=tmpd6[4]; + } + } + sx=0; sy=0;nall=0; + for(ii=1;ii<=length6(tmp1md6);ii++){ + pull6(ii,tmp1md6,tmpd6); + if(tmpd6[4]<dst+Eps00){ + nall++; + sx=sx+tmpd6[0]; sy=sy+tmpd6[1]; + } + tmp1d[0]=sx/nall; tmp1d[1]=sy/nall; + tmp1d6[0]=tmp1d[0];tmp1d6[1]=tmp1d[1]; + nearestptpt(tmp1d,crv1,ans); + tmp1d6[2]=ans[2]; + nearestptpt(tmp1d,crv2,ans); + tmp1d6[3]=ans[2]; + tmp1d6[4]=dst; + tmp1d6[5]=tmpd6[5]; + addptd6(out,tmp1d6); + } + tmp1md6[0][0]=0; + for(ii=1;ii<=length6(tmp2md6);ii++){ + pull6(ii,tmp2md6,tmpd6);addptd6(tmp1md6,tmpd6); + } + } +} + +void intersectcurves(double crv1[][2],double crv2[][2],int num,double out[][2]){ + double out6[DsizeS][6],tmpd[2],tmpd6[6]; + int jj; + intersectcurvesPp(crv1,crv2,num,out6); + out[0][0]=0; + for(jj=1;jj<=length6(out6);jj++){ + pull6(jj,out6,tmpd6);push2(tmpd6[0],tmpd6[1],out); + } +} |