Subversion Repositories wimsdev

Rev

Rev 8177 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. /*    Copyright (C) 1998-2003 XIAO, Gang of Universite de Nice - Sophia Antipolis
  2.  *
  3.  *  This program is free software; you can redistribute it and/or modify
  4.  *  it under the terms of the GNU General Public License as published by
  5.  *  the Free Software Foundation; either version 2 of the License, or
  6.  *  (at your option) any later version.
  7.  *
  8.  *  This program is distributed in the hope that it will be useful,
  9.  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10.  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11.  *  GNU General Public License for more details.
  12.  *
  13.  *  You should have received a copy of the GNU General Public License
  14.  *  along with this program; if not, write to the Free Software
  15.  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  16.  */
  17.  
  18.  /* Finds the shortest paths linking given points. */
  19.  
  20. /* fast algorithm is not ready. */
  21.  
  22. #define MAX_POINTS 32 /* cannot exceed number of bits in long */
  23. #define MAX_DIM 10 /* dimension limit */
  24. #define MAX_SHORT 20 /* limit of registered shortest paths */
  25.  
  26.  
  27. #include <stdio.h>
  28. #include <stdlib.h>
  29. #include <string.h>
  30. #include <math.h>
  31. #include <ctype.h>
  32. #include "../config.h"
  33.  
  34. int  pointcnt,dimension,first,last,count,shortcnt;
  35. /* style = 0: loop to the start
  36.  *    1: arbitrary open path
  37.  *    2: open path with fixed start
  38.  *    3: open path with fixed end
  39.  *    4: open path with fixed start and end */
  40. int style=0;
  41. int fast=0; /* fast algo if 1, no rigor */
  42. int optimize=0;
  43. char *input;
  44. double dist[MAX_POINTS][MAX_POINTS];
  45. double ds[MAX_POINTS];
  46. double coord[MAX_POINTS][MAX_DIM];
  47. double shortest;
  48. double tolerance; /* to compensate floating rounding error in comparison */
  49. double distmin; /* minimum of distance between two points */
  50. double firstlast; /* distance between first and last */
  51. int shpath[MAX_SHORT][MAX_POINTS];
  52. int chain[MAX_POINTS];
  53.  
  54. /* reverse a subchain */
  55. void reverse(int i, int j)
  56. {
  57.   int k,t;
  58.   for(k=0;k<(j-i+1)/2;k++) {
  59.     t=chain[i+k]; chain[i+k]=chain[j-k];
  60.     chain[j-k]=t;
  61.   }
  62. }
  63.  
  64. /* move point i to position j */
  65. void reinsert(int i, int j)
  66. {
  67.   int t;
  68.   t=chain[i];
  69.   if(j>i) {
  70.     j--;
  71.     memmove(chain+i,chain+i+1,(j-i)*sizeof(chain[0]));
  72.     chain[j]=t;
  73.   }
  74.   else {
  75.     memmove(chain+j+1,chain+j,(i-j)*sizeof(chain[0]));
  76.     chain[j]=t;
  77.   }
  78. }
  79.  
  80. int _fcalc(void)
  81. {
  82.   int modif;
  83.   int i,j;
  84.  
  85.   modif=0;
  86. /* four-configuration optimize */
  87.   for(i=1;i<pointcnt-2;i++) {
  88.     for(j=i+1;j<pointcnt-1;j++) {
  89.         if(dist[chain[i-1]][chain[i]]+dist[chain[j]][chain[j+1]]
  90.            >dist[chain[i-1]][chain[j]]+dist[chain[i]][chain[j+1]]
  91.            +tolerance) {
  92.           reverse(i,j); modif=1;
  93.         }
  94.     }
  95.   }
  96. /* Point relink */
  97.   for(i=1;i<pointcnt-1;i++) {
  98.     double dd;
  99.     int th;
  100.     th=chain[i];
  101.     dd=dist[chain[i-1]][th]+dist[th][chain[i+1]]
  102.       -dist[chain[i-1]][chain[i+1]];
  103.     for(j=0;j<pointcnt-1;j++) {
  104.       if(j==i || j==i-1) continue;
  105.       if(dd>dist[chain[j]][th]+dist[chain[j+1]][th]
  106.            -dist[chain[j]][chain[j+1]]+tolerance) {
  107.         reinsert(i,j+1); modif=1;
  108.       }
  109.     }
  110.   }
  111.   for(i=0;i<pointcnt;i++) shpath[shortcnt][i]=chain[i];
  112.   shortcnt+=modif; return modif;
  113. }
  114.  
  115. void fastcalc(void)
  116. {
  117.   int i,j;
  118.  
  119.   for(i=0;i<pointcnt;i++) chain[i]=i;
  120.   shortcnt=1;
  121.   for(i=0,j=1;i<10 && j!=0;i++) j=_fcalc();
  122.  
  123.   for(shortest=0,i=0;i<pointcnt-1;i++) shortest+=dist[chain[i]][chain[i+1]];
  124.   if(style==0) shortest+=dist[chain[0]][chain[pointcnt-1]];
  125.   for(i=0;i<pointcnt;i++) shpath[0][i]=chain[i];
  126. }
  127.  
  128. void _calc(int level, double dis, long used, int utab[])
  129. {
  130.   int current[MAX_POINTS];
  131.   double dd,d1,d3,d2,dthis;
  132.   int i;
  133.   long l;
  134.  
  135.   dd=dis; d1=d2=d3=0;
  136.   if(level<pointcnt-2) {
  137.     int j,b,prec;
  138.     double btab[MAX_POINTS];
  139.     if(level>0) prec=utab[level-1]; else prec=0;
  140.     if(level>2) {
  141.       for(i=0;i<level-2;i++)
  142.         btab[i]=dist[utab[i]][prec];
  143.       if(level>3) {
  144.         d1=ds[level-4]+ds[level-3]+ds[level-2];
  145.         d3=dist[utab[level-4]][utab[level-2]]
  146.           +ds[level-2]
  147.           +dist[prec][utab[level-3]];
  148.         d2=dist[utab[level-4]][prec]
  149.           +dist[prec][utab[level-3]]
  150.           +ds[level-3];
  151.       }
  152.   }
  153.   if(level==0 || (level==1 && style==0)) b=last+1; else b=0;
  154.   for(i=b,l=(long) 1<<b;i<pointcnt;i++,l<<=1) {
  155.     if(used&l) continue;
  156.     dthis=dist[prec][i];
  157.     if(optimize) {
  158.       if(style==0) {
  159.         if(level>1 &&
  160.              firstlast+dthis>dist[prec][last]+dist[first][i])
  161.           continue;
  162.         }
  163.         else if(level>0) {
  164. /* cut differently */
  165.           if(style==1 && dthis>firstlast) continue;
  166. /* switch with first */
  167.           if((style&1) && level>1 && dist[first][i]<dthis) continue;
  168. /* switch with last */
  169.           if(style<3 && dist[prec][last]<dthis) continue;
  170.         }
  171. /* already too long */
  172.         if(level>0 && shortcnt>0 &&
  173.              dthis+dis+distmin*(pointcnt-level-1)>shortest)
  174.           continue;
  175. /* four-config optimize */
  176.         for(j=0;j<=level-3;j++) {
  177.           if(ds[j]+dthis>btab[j]+dist[utab[j+1]][i])
  178.             goto next;
  179.         }
  180. /* five-config optimize */
  181.         if(level>3 &&
  182.              (d1+dthis>d2+dist[utab[level-2]][i] ||
  183.             d1+dthis>d3+dist[utab[level-3]][i]))
  184.           continue;
  185. /* try prec elsewhere in the chain */
  186.         if(level>3) {
  187.           double dt;
  188.           dt=dist[prec][i]+ds[level-2]-dist[utab[level-2]][i];
  189.           for(j=0;j<level-3;j++) {
  190.               if(dt>dist[prec][utab[j]]+dist[prec][utab[j+1]]-ds[j])
  191.                 goto next;
  192.           }
  193.         }
  194.       }
  195.       memmove(current,utab,level*sizeof(utab[0]));
  196.       current[level]=i;
  197.       if(level>0) {dd=dis+dthis; ds[level-1]=dthis;}
  198.       else {
  199.         first=i; firstlast=dist[first][last];
  200.       }
  201.       _calc(level+1,dd,used|l,current);
  202.       next: ;
  203.     }
  204.  
  205.   }
  206.   else {
  207.     count++;
  208.     for(i=0,l=(long) 1;i<pointcnt && used&l; i++,l<<=1);
  209.     if(i>=pointcnt) {
  210.       printf("ERROR: Internal error.\n"); exit(1);
  211.     }
  212.     utab[pointcnt-2]=i; utab[pointcnt-1]=last;
  213.     dis+=dist[utab[pointcnt-3]][i]+dist[i][last];
  214.     if(shortest==0 || dis<shortest-tolerance) {
  215.       shortest=dis; shortcnt=0;
  216.       memmove(shpath[shortcnt++],utab,MAX_POINTS*sizeof(utab[0]));
  217.     }
  218.     else if(dis<shortest+tolerance && shortcnt<MAX_SHORT)
  219.       memmove(shpath[shortcnt++],utab,MAX_POINTS*sizeof(utab[0]));
  220.   }
  221. }
  222.  
  223. /* main calculation loops */
  224. void calc(void)
  225. {
  226.   int i;
  227.   long used;
  228.   int utab[MAX_POINTS];
  229.  
  230.   if(fast) {
  231.     fastcalc(); return;
  232.   }
  233.   count=0; shortest=0; shortcnt=0;
  234.   memset(utab,0,sizeof(utab));
  235.   switch(style) {
  236.     case 0: { /* loop to the start */
  237.       first=0; utab[0]=0;
  238.       for(i=1;i<pointcnt;i++) {
  239.         last=i; used=(long) 1<<i; used|=(long) 1;
  240.         firstlast=dist[first][last];
  241.         _calc(1,firstlast,used,utab);
  242.       }
  243.       break;
  244.     }
  245.     case 1: { /* arbitrary open path */
  246.       for(i=0;i<pointcnt;i++) {
  247.         last=i; used=(long) 1<<i;
  248.         _calc(0,0,used,utab);
  249.       }
  250.       break;
  251.     }
  252.     case 2: { /* open path with fixed start */
  253.       first=0; utab[0]=0;
  254.         for(i=1;i<pointcnt;i++) {
  255.           last=i; used=(long) 1|(1<<i);
  256.           firstlast=dist[0][last];
  257.           _calc(1,0,used,utab);
  258.         }
  259.         break;
  260.     }
  261.     case 3: { /* open path with fixed end */
  262.       last=0; used=(long) 1;
  263.       _calc(0,0,used,utab);
  264.       break;
  265.     }
  266.     case 4: { /* open path with fixed start and end */
  267.       first=0; utab[0]=0;
  268.       last=pointcnt-1; used=(long) 1|(1<<last);
  269.       firstlast=dist[first][last];
  270.       _calc(1,0,used,utab);
  271.       break;
  272.     }
  273.   }
  274. }
  275.  
  276. /* Print the result */
  277. void printresult(void)
  278. {
  279.   int i,j;
  280.  
  281.   printf("%f %d %d\n",shortest,shortcnt,count);
  282.   for(j=0;j<shortcnt;j++) {
  283.     for(i=0;i<pointcnt;i++) printf("%d ",shpath[j][i]+1);
  284.     printf("\n");
  285.   }
  286. }
  287.  
  288. /* prepare distance table */
  289. void makedist(void)
  290. {
  291.   int i,j,k;
  292.   double d;
  293.  
  294.   distmin=-1;
  295.   for(i=0;i<pointcnt-1;i++) {
  296.     for(j=i+1;j<pointcnt;j++) {
  297.       d=0;
  298.       for(k=0;k<dimension;k++) {
  299.         d+=pow(coord[i][k]-coord[j][k],2);
  300.       }
  301.       d=sqrt(d); dist[i][j]=dist[j][i]=d;
  302.       if(distmin<0 || distmin>d) distmin=d;
  303.     }
  304.   }
  305.   tolerance=distmin*0.000001;
  306. }
  307.  
  308. int main(int argc, char *argv[])
  309. {
  310.   char *p,*p2,*p3;
  311.   int dim;
  312.   input=getenv("wims_exec_parm");
  313. /* nothing to do if no problem */
  314.   if(input==NULL || *input==0) return 0;
  315.  
  316.   for(p=input+strlen(input);p>input && isspace(*(p-1));p--);
  317.   *p=0; dimension=0;
  318.   for(pointcnt=0,p=input;*p!=0 && pointcnt<MAX_POINTS;p=p3) {
  319.     p2=strchr(p,'\n');
  320.     if(p2==NULL) p3=p+strlen(p);
  321.     else {p3=p2+1; *p2=0;}
  322.     dim=0;
  323.     for(p2=strchr(p,','); p2!=NULL && dim<MAX_DIM;
  324.         p=p2+1,dim++,p2=strchr(p,',')) {
  325.       *p2=0; coord[pointcnt][dim]=strtod(p,NULL);
  326.     }
  327.     if(dim<MAX_DIM) coord[pointcnt][dim++]=strtod(p,NULL);
  328.     if(dim<2) continue; /* bad line */
  329.     if(dimension==0) dimension=dim;
  330.     pointcnt++;
  331.   }
  332.   if(pointcnt<3) {
  333.     printf("ERROR: Not enough points.\n"); exit(1);
  334.   }
  335.   p=getenv("w_shortpath_style");
  336.   if(p!=NULL) style=atoi(p);
  337.   if(style<0 || style>4) style=0;
  338.   makedist();
  339.   optimize=1; calc();printresult();
  340. /*    optimize=0; calc();printresult();
  341. */
  342.   return 0;
  343. }
  344.  
  345.