Back to index

wims  3.65+svn20090927
curvecomp.c
Go to the documentation of this file.
00001 /*    Copyright (C) 1998-2003 XIAO, Gang of Universite de Nice - Sophia Antipolis
00002  *
00003  *  This program is free software; you can redistribute it and/or modify
00004  *  it under the terms of the GNU General Public License as published by
00005  *  the Free Software Foundation; either version 2 of the License, or
00006  *  (at your option) any later version.
00007  *
00008  *  This program is distributed in the hope that it will be useful,
00009  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00010  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011  *  GNU General Public License for more details.
00012  *
00013  *  You should have received a copy of the GNU General Public License
00014  *  along with this program; if not, write to the Free Software
00015  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00016  */
00017 
00018        /* Compare two curves */
00019 /* Input parameters: environment.
00020  * 
00021  * w_curvecomp_1 and w_curvecomp_2: curves to compare, as lists of points.
00022  * w_curvecomp_xrange and w_curvecomp_yrange: list of 2 integers each.
00023  * w_curvecomp_tolerance: Maximal tolerance of distances.
00024  * 
00025  * Output: 10 double numbers separated by white spaces.
00026  * 
00027  * Number 1: Average distance of curve 1 with respect to curve 2.
00028  * Number 2: Average distance of curve 2 with respect to curve 1.
00029  * Number 3: Maximal distance of curve 1 with respect to curve 2.
00030  * Number 4: Maximal distance of curve 2 with respect to curve 1.
00031  * Number 5: Proportion of curve 1 close to curve 2.
00032  * Number 6: Proportion of curve 2 close to curve 1.
00033  * Number 7: Maximal jump of curve 1.
00034  * Number 8: Maximal jump of curve 2.
00035  * Number 9: Ratio of repetitions found in curve 1.
00036  * Number 10: Ratio of repetitions found in curve 2.
00037  * Furthermore, words "fnofx" and/or "fnofy" will appear if curve 2
00038  *   represents the graph of a function of x (and/or y).
00039  * 
00040  * Returns empty if one of the curves is degenerated.
00041 */
00042 
00043 /*************** Customization: change values hereafter ****************/
00044 
00045        /* limit of point buffers */
00046 #define pointlim     4096
00047 #define default_tolerance 10
00048 
00049 /***************** Nothing should need change hereafter *****************/
00050 
00051 #include "../wims.h"
00052 char curve1[2*MAX_LINELEN+2], curve2[2*MAX_LINELEN+2];
00053 int bx[2], by[2];
00054 int tol=default_tolerance;
00055 int xfn=0, yfn=0;
00056 int discrete;
00057 double dist1, dist2, max1, max2, ratio1, ratio2, rep1, rep2;
00058 double djump1, djump2;
00059 struct cv {
00060     int x,y, closest, rep;
00061     double dist;
00062 } cv1[pointlim], cv2[pointlim];
00063 int points1,points2;
00064 
00065 int Abs(int t) {if(t>=0) return t; else return -t;}
00066 int Min(int x,int y) {if(x>y) return y; else return x;}
00067 int Max(int x,int y) {if(x<y) return y; else return x;}
00068 
00069 int listbuf[pointlim*2], listcnt;
00070 
00071        /* Strips leading spaces */
00072 char *find_word_start(char *p)
00073 {
00074     int i;
00075     for(i=0; isspace(*p) && i<MAX_LINELEN; p++,i++);
00076     return p;
00077 }
00078 
00079 void reverse(struct cv *cvbuf, int cnt)
00080 {
00081     int i;
00082     struct cv cvt;
00083     for(i=0;i<cnt/2;i++) {
00084        memmove(&cvt,cvbuf+i,sizeof(cvt));
00085        memmove(cvbuf+i,cvbuf+(cnt-i-1),sizeof(cvt));
00086        memmove(cvbuf+(cnt-i-1),&cvt,sizeof(cvt));
00087     }
00088 }
00089 
00090 int str2list(char *p, int lim)
00091 {
00092     char *p2;
00093     for(p2=strchr(p,';'); p2; p2=strchr(p2+1,';')) *p2=',';
00094     for(listcnt=0; p && listcnt<lim; p=p2) {
00095        p2=strchr(p,','); if(p2!=NULL) *p2++=0;
00096        p=find_word_start(p); if(*p==0) continue;
00097        listbuf[listcnt++]=atoi(p); 
00098     }
00099     return listcnt;
00100 }
00101 
00102 int list2curve(struct cv *cvbuf)
00103 {
00104     int i, j, m, t, st, xx, yy, x2, y2, x3, y3, ll;
00105 
00106     ll=listcnt/2; if(ll<2) return 0;
00107     xx=listbuf[0]; yy=listbuf[1];
00108     j=0; if(xx>=bx[0] && xx<=bx[1] && yy>=by[0] && yy<=by[1]) {
00109        cvbuf[0].x=xx; cvbuf[0].y=yy; j++;
00110     }
00111     for(i=1; i<ll && j<pointlim; i++) {
00112        x2=listbuf[2*i]; y2=listbuf[2*i+1];
00113        m=Max(Abs(x2-xx),Abs(y2-yy)); if(m<=0) continue;
00114        if(discrete==1) st=m; else st=1;
00115        for(t=st; t<=m && j<pointlim; t++) {
00116            x3=(double) (x2*t+xx*(m-t))/m+0.5;
00117            y3=(double) (y2*t+yy*(m-t))/m+0.5;
00118            if(x3>=bx[0] && x3<=bx[1] && y3>=by[0] && y3<=by[1]) {
00119               cvbuf[j].x=x3; cvbuf[j].y=y3; cvbuf[j].dist=-1;
00120               cvbuf[j].rep=0; j++;
00121            }
00122        }
00123        xx=x2; yy=y2;
00124     }
00125     return j;
00126 }
00127 
00128 void compare(void)
00129 {
00130     int i, j, cl;
00131     double d, dt;
00132     
00133     d=2*pointlim; cl=-1;
00134     for(i=0,djump1=1;i<points1-1;i++) {
00135        dt=sqrt(pow(cv1[i].x-cv1[i+1].x,2)+pow(cv1[i].y-cv1[i+1].y,2));
00136        if(dt>djump1) djump1=dt;
00137     }
00138     for(i=0,djump2=1;i<points2-1;i++) {
00139        dt=sqrt(pow(cv2[i].x-cv2[i+1].x,2)+pow(cv2[i].y-cv2[i+1].y,2));
00140        if(dt>djump2) djump2=dt;
00141     }
00142     for(i=0;i<points1;i++) {
00143        for(j=0;j<points2;j++) {
00144            dt=sqrt(pow(cv2[j].x-cv1[i].x,2)+pow(cv2[j].y-cv1[i].y,2));
00145            if(dt<d) {d=dt; cl=j;}
00146            else {dt=(dt-d)/djump2; if(dt>2) j+=dt-1;}
00147        }
00148        cv1[i].dist=d; cv1[i].closest=cl;
00149        if(i<points1)
00150          d+=sqrt(pow(cv1[i].x-cv1[i+1].x,2)+pow(cv1[i].y-cv1[i+1].y,2))+0.1;
00151     }
00152     d=2*pointlim; for(i=0;i<points2;i++) {
00153        for(j=0;j<points1;j++) {
00154            dt=sqrt(pow(cv1[j].x-cv2[i].x,2)+pow(cv1[j].y-cv2[i].y,2));
00155            if(dt<d) {d=dt; cl=j;}
00156            else {dt=(dt-d)/djump1; if(dt>2) j+=dt-1;}
00157        }
00158        cv2[i].dist=d; cv2[i].closest=cl;
00159        if(i<points2)
00160          d+=sqrt(pow(cv2[i].x-cv2[i+1].x,2)+pow(cv2[i].y-cv2[i+1].y,2))+0.1;
00161     }
00162     for(i=1, cl=cv1[0].closest;i<points1;i++) {
00163        j=cv1[i].closest; if(j!=cl) cv2[j].rep++;
00164        cl=j;
00165     }
00166     for(i=1, cl=cv2[0].closest;i<points2;i++) {
00167        j=cv2[i].closest; if(j!=cl) cv1[j].rep++;
00168        cl=j;
00169     }
00170     for(i=cl=0; i<points1; i++) if(cv1[i].rep>1) cl+=cv1[i].rep-1;
00171     rep1=(double) cl/points1;
00172     for(i=cl=0; i<points1; i++) if(cv2[i].rep>1) cl+=cv2[i].rep-1;
00173     rep2=(double) cl/points2;
00174 }
00175 
00176 void check(void)
00177 {
00178     int i,j,xret1,yret1,xret2,yret2;
00179     int xx,yy,xmin,xmax,ymin,ymax;
00180     double d;
00181     max1=max2=0;
00182     for(i=j=0,d=0;i<points1;i++) {
00183        if(cv1[i].dist<=tol) j++;
00184        d+=pow(cv1[i].dist,4);
00185        if(max1<cv1[i].dist) max1=cv1[i].dist;
00186     }
00187     ratio1=(double) j/points1; dist1=pow(d/points1,0.25)*1.8;
00188     for(i=j=0,d=0;i<points2;i++) {
00189        if(cv2[i].dist<=tol) j++;
00190        d+=pow(cv2[i].dist,4);      
00191        if(max2<cv2[i].dist) max2=cv2[i].dist;
00192     }
00193     ratio2=(double) j/points2; dist2=pow(d/points2,0.25)*1.8;
00194     xret1=xret2=yret1=yret2=0;
00195     xmin=xmax=cv2[0].x; ymin=ymax=cv2[0].y;
00196     for(i=1;i<points2;i++) {
00197        xx=cv2[i].x; yy=cv2[i].y;
00198        xret1=Max(xret1,xmax-xx); xret2=Max(xret2,xx-xmin);
00199        yret1=Max(yret1,ymax-yy); yret2=Max(yret2,yy-ymin);
00200        xmin=Min(xx,xmin);xmax=Max(xx,xmax);
00201        ymin=Min(yy,ymin);ymax=Max(yy,ymax);
00202     }
00203     if(Min(xret1,xret2)<=2) xfn=1;
00204     if(Min(yret1,yret2)<=2) yfn=1;
00205 }
00206 
00207 void output(void)
00208 {
00209     printf("%.2f %.2f %.2f %.2f \
00210 %.3f %.3f %.1f %.1f \
00211 %.3f %.3f",
00212           dist1, dist2, max1, max2,
00213           ratio1, ratio2, djump1, djump2,
00214           rep1, rep2);
00215     if(xfn) printf(" fx");
00216     if(yfn) printf(" fy");
00217 }
00218 
00219 void test(void)
00220 {
00221     int i;
00222     for(i=0;i<points1;i++) printf("%d,%d,%d,%.2f\n",
00223                               cv1[i].x,cv1[i].y,cv1[i].closest,cv1[i].dist);
00224     printf("\n");
00225     for(i=0;i<points2;i++) printf("%d,%d,%d,%.2f\n",
00226                               cv2[i].x,cv2[i].y,cv2[i].closest,cv2[i].dist);
00227     printf("\n");
00228 }
00229 
00230 int main(int argc, char *argv[])
00231 {
00232     int i;
00233     char *c1, *c2, *op;
00234     c1=getenv("w_curvecomp_1"); c2=getenv("w_curvecomp_2");
00235        /* nothing to do */
00236     if(c1==NULL || c2==NULL || *c1==0 || *c2==0) return 0;
00237     snprintf(curve1,sizeof(curve1),"%s",c1);
00238     snprintf(curve2,sizeof(curve2),"%s",c2);
00239     bx[0]=by[0]=bx[1]=by[1]=-1;
00240     c1=getenv("w_curvecomp_xrange"); c2=getenv("w_curvecomp_yrange");
00241     if(c1!=NULL && *c1!=0) {
00242        str2list(c1,2); for(i=0;i<listcnt;i++) bx[i]=listbuf[i];
00243     }
00244     if(c2!=NULL && *c2!=0) {
00245        str2list(c2,2); for(i=0;i<listcnt;i++) by[i]=listbuf[i];
00246     }
00247     op=getenv("w_curvecomp_option"); if(op==NULL) op="";
00248     if(bx[0]==-1) bx[0]=0; if(bx[1]==-1) bx[1]=pointlim;
00249     if(by[0]==-1) by[0]=0; if(by[1]==-1) by[1]=pointlim;
00250     c1=getenv("w_curvecomp_tolerance");
00251     if(c1!=NULL && *c1!=0) tol=atoi(c1);
00252     tol=Min(30,Max(5,tol));
00253     
00254     if(strstr(op,"discrete1")!=NULL) discrete=1; else discrete=0;
00255     str2list(curve1,pointlim*2); points1=list2curve(cv1);
00256     if(strstr(op,"discrete2")!=NULL) discrete=1; else discrete=0;
00257     str2list(curve2,pointlim*2); points2=list2curve(cv2);
00258     if(points1<2 || points2<2) return 0;
00259     compare(); check(); output();
00260     return 0;
00261 }
00262