Subversion Repositories wimsdev

Rev

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