home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / libip / Copy of Descript.c.083098 < prev    next >
Encoding:
Text File  |  1999-09-11  |  8.7 KB  |  373 lines
  1. /* 
  2. ** descript.c
  3. ** 
  4. ** Practical Algorithms for Image Analysis
  5. ** 
  6. ** Copyright (c) 1997 ????
  7. */
  8.  
  9. /*
  10. ** DESCRIPT
  11. **
  12. **  evaluation of Fourier descriptors for a plane curve
  13. **    encoded in a set {delta_phi(k), l(k)};
  14. **
  15. **  references:
  16. **  C. T. Zahn, 
  17. **  Proc. Joint Intern. Conf. on Artif. Intelligence, May 1969, pp. 621 - 628;
  18. **  SLAC Report 70, Stanford Linear Accelerator Center, 1968;
  19. **  C. T. Zahn & R. Z. Roskies, IEEE Trans. Comp., Vol C-21, 269 - 281 (1972);
  20. **
  21. **  -->the truncated Fourier expansion derived by Z & R is computed 
  22. **     in form of inner products using array processor functions
  23. **
  24. **  -->test shapes (see Z & R): 
  25. **        four (xdrawfour.c), inverted L (fdzahnl.c)
  26. **
  27. ** ---------------------------------------------------------------------------
  28. **
  29. */
  30. #include <stdio.h>
  31. #include <stdlib.h>
  32. #include <math.h>
  33. #include "ip.h"
  34.  
  35. #define    EPS        0.00000001
  36. #define    SQ2        1.414213562
  37.  
  38. #define    ON        1
  39. #define    OFF        0
  40. #define    FD_DEBUG    ON
  41.  
  42.  
  43. /*
  44.  * descriptors()
  45.  *   DESCRIPTION:
  46.  *     evaluate ZAHN-ROSKIES Fourier descriptors
  47.  *   ARGUMENTS:
  48.  *     dphik: delta phik for polygon (float *)
  49.  *     dlk: delta lk for polygon (float *)
  50.  *     mcp: length of dphik and dlk (long)
  51.  *     a_n: Fourier coefficients array (float *)
  52.  *     b_n: Fourier coefficients array (float *)
  53.  *     n_order: order of Fourier descriptors to compute (long)
  54.  *   RETURN VALUE:
  55.  *     none
  56.  */
  57.  
  58. void 
  59. descriptors (float *dphik, float *dlk, long mcp, float *a_n, float *b_n, long n_order)
  60. {
  61.     int k, n_ord;
  62.     int fund_flag = OFF;
  63.     
  64.     float real_length, length;
  65.     float dc_term, a, b, part_sum;
  66.     float sne1k, snen, csn1k, csnn;
  67.     float *angle;
  68.     float *sne1, *sne, *csn1, *csn;
  69.     float *mod, *arg;
  70.  
  71.  
  72. /*  
  73. **  compute (in place) running sums of the elements in dlk 
  74. **  and the corresponding phase angles required in the evaluation
  75. **  of the fundamental
  76. */
  77.     angle = (float *)calloc((int)mcp, sizeof(float));
  78.     sne1 = (float *)calloc((int)mcp, sizeof(float));
  79.     csn1 = (float *)calloc((int)mcp, sizeof(float));
  80.     sne = (float *)calloc((int)mcp, sizeof(float));
  81.     csn = (float *)calloc((int)mcp, sizeof(float));
  82.     if((angle == NULL) || (sne == NULL) || (csn1 == NULL) || (sne == NULL)
  83.         || (csn == NULL))
  84.     {
  85.         printf("memory allocation error in function descriptors\n");
  86.         exit(1);
  87.     }
  88.  
  89. /*
  90. ** evaluate contour length, based on curvature points
  91. */
  92.     part_sum = *(dlk+0);
  93.     for(k=1; k<mcp; k++)
  94.     {
  95.         part_sum += *(dlk+k);
  96.         *(dlk+k) = part_sum;        /* array of partial sums */
  97.     }
  98.     length = *(dlk+mcp-1);
  99.         
  100.     for(k=0; k<mcp; k++)
  101.         *(angle+k) = (float)(2*PI*(*(dlk+k))/length);
  102.  
  103.     real_length = (float)(0.5*length);
  104.     printf("\n...length = %f\n", real_length);
  105.  
  106. /*
  107. **  dphik contains elements of delta_phik multiplied by PI/4;
  108. **  to speed up computation, factor out PI 
  109. */
  110.     for(k=0; k<mcp; k++)
  111.         *(dphik+k) *= (float)0.25;
  112.  
  113. #ifdef FD_DEBUG
  114.     printf("...elements of dphik and dlk:\n");
  115.     for(k=0; k<mcp; k++)
  116.         printf("%7.4f  %7.4f\n", *(dphik+k), *(dlk+k) );
  117. #endif
  118.     
  119. /*
  120.  *  compute DC term
  121.  */
  122.     vdot(dphik, dlk, &dc_term, mcp);
  123.     dc_term = (float)(-PI*(1.0 + dc_term/length));
  124.     printf("...and we have a result for the DC term: %f   \n", dc_term);
  125.  
  126. /*
  127. **  execute routine to compute sine and
  128. **  cosine terms for fundamental, then fetch results
  129. */
  130.     for(n_ord=0; n_ord<(int)n_order; n_ord++)
  131.     {
  132.  
  133.     /*  
  134.      *  compute coefficients for fundamental;
  135.      */
  136.     /*    if(n_ord == 0)
  137.         {
  138.         */
  139.             for(k = 0; k < mcp; k++)
  140.                 *(angle + k) = *(angle + k) * (n_ord + 1);
  141.             vsin(angle, sne1, mcp);
  142.             vcos(angle, csn1, mcp);
  143.             for (k=0; k<mcp; k++)
  144.             {
  145.                 csn1k = *(csn1+k);
  146.                 sne1k = *(sne1+k);
  147.                 if( fabs(csn1k) < EPS )  csn1k = (float)0.0;
  148.                 if( fabs(sne1k) < EPS )  sne1k = (float)0.0;
  149.             
  150.                 *(sne+k) = *(sne1+k) = sne1k;
  151.                 *(csn+k) = *(csn1+k) = csn1k;
  152.             }
  153.             fund_flag = ON;
  154.         /*}*/
  155.  
  156.  
  157. /*
  158. **  employ double angle formulas to compute higher order coefficients
  159. */
  160. #if 0
  161.         else
  162.         {
  163.             if( fund_flag != ON )
  164.             {
  165.                 printf("\n");
  166.                 printf("hey, evaluate fundamental first!\n");
  167.                 printf(" something wrong!!\n");
  168.                 exit(1);
  169.             }
  170.             for (k=0; k<mcp; k++)
  171.             {
  172.                 snen = *(sne+k);
  173.                 csnn = *(csn+k);
  174.                 *(sne+k) = snen*(*(csn1+k))+csnn*(*(sne1+k));
  175.                 *(csn+k) = csnn*(*(csn1+k))-snen*(*(sne1+k));
  176.             }
  177. #ifdef FD_DEBUG
  178.             printf("\n...n_ord = %d:\n", n_ord);
  179.             for(k=0; k<10; k++)
  180.                 printf("......sne[%d] = %f   csn[%d] = %f\n",
  181.                         k, *(sne+k), k, *(csn+k) );
  182. #endif
  183.         }
  184. #endif
  185. #ifdef FD_DEBUG
  186.             printf("\n...n_ord = %d:\n", n_ord);
  187.             for(k=0; k<10; k++)
  188.                 printf("......sne[%d] = %f   csn[%d] = %f\n",
  189.                         k, *(sne+k), k, *(csn+k) );
  190. #endif
  191.  
  192.  
  193. /*
  194. ** compute coefficients a_n, b_n
  195. */
  196.  
  197.         vdot(dphik, sne, &a, mcp);
  198.         vdot(dphik, csn, &b, mcp);
  199.  
  200.         *(a_n+n_ord) = -a/(float)(n_ord+1); 
  201.         *(b_n+n_ord) = b/(float)(n_ord+1);
  202.     }
  203.  
  204.     printf("...descriptors up to %ld-th order are:\n", n_order);
  205.     for(n_ord=0; n_ord<(int)n_order; n_ord++)
  206.         printf ("%7.3f   %7.3f\n", *(a_n+n_ord), *(b_n+n_ord));
  207.  
  208.  
  209.     mod = (float *)calloc((int)n_order, sizeof(float));
  210.     arg = (float *)calloc((int)n_order, sizeof(float));
  211. /*
  212. **  calculate xy to polar coordinate transformation
  213. */
  214.     vrtop(a_n, b_n, mod, arg, n_order);
  215.     printf("...polar descriptors up to %ld-th order are:\n", n_order);
  216.     for(n_ord=0; n_ord<(int)n_order; n_ord++)
  217.     {
  218.         if( *(arg+n_ord) < 0.0)     *(arg+n_ord) += (float)(2.0*PI);
  219.         printf("%7.3f   %7.3f\n", *(mod+n_ord), *(arg+n_ord));
  220.     }
  221.  
  222.     free(angle);
  223.     free(sne1);
  224.     free(csn1);
  225.     free(sne);
  226.     free(csn);
  227. }
  228.  
  229. /*
  230.  * vdot()
  231.  *   DESCRIPTION:
  232.  *     take the dot product of 2 vectors
  233.  *     and compute the sum
  234.  *   ARGUMENTS:
  235.  *     v1: first vector (float *)
  236.  *     v2: second vector (float *)
  237.  *     vsum: vector sum (float *)
  238.  *     vlen: length of vectors (long)
  239.  *   RETURN VALUE:
  240.  *     none
  241.  */
  242.  
  243. void
  244. vdot(float *v1, float *v2, float *vsum, long vlen)
  245. {
  246.     float rsum = (float)0.0;
  247.     int i;
  248.  
  249.     for(i = 0; i < vlen; i++)
  250.         rsum = rsum + *(v1 + i) * (*(v2 + i));
  251.     *vsum = rsum;
  252. }
  253.  
  254. /*
  255.  * vcos()
  256.  *   DESCRIPTION:
  257.  *     takes the cosine of a vector
  258.  *   ARGUMENTS:
  259.  *     source: input vector (float *) NOTE: radians!!
  260.  *     dest: output vector (float *)
  261.  *     vlen: vector length (long)
  262.  *   RETURN VALUE:
  263.  *     none
  264.  */
  265.  
  266. void
  267. vcos(float *source, float *dest, long vlen)
  268. {
  269.     int i;
  270.  
  271.     for(i = 0; i < vlen; i++)
  272.         *(dest + i) = (float)cos((double)*(source + i));
  273. }
  274.  
  275. /*
  276.  * vsin()
  277.  *   DESCRIPTION:
  278.  *     takes the sine of a vector
  279.  *   ARGUMENTS:
  280.  *     source: input vector (float *) NOTE: radians!!
  281.  *     dest: output vector (float *)
  282.  *     vlen: vector length (long)
  283.  *   RETURN VALUE:
  284.  *     none
  285.  */
  286.  
  287. void
  288. vsin(float *source, float *dest, long vlen)
  289. {
  290.     int i;
  291.  
  292.     for(i = 0; i < vlen; i++)
  293.         *(dest + i) = (float)sin((double)*(source + i));
  294. }
  295.  
  296. /*
  297.  * vrtop()
  298.  *   DESCRIPTION:
  299.  *     converts rectangular cartesian
  300.  *     coordinates to polar
  301.  *   ARGUMENTS:
  302.  *     x: input x vector (float *)
  303.  *     y: input y vector (float *)
  304.  *     mod: modulus output vector (float *)
  305.  *     arg: atan2 output vector (float *)
  306.  *     vlen: vector lengths (long)
  307.  *   RETURN VALUE:
  308.  *     none
  309.  */
  310.  
  311. void
  312. vrtop(float *x, float *y, float *mod, float *arg, long vlen)
  313. {
  314.     int i;
  315.  
  316.     for(i = 0; i < vlen; i++) {
  317.         *(mod + i) = (*(x + i))*(*(x + i)) + (*(y + i))*(*(y + i));
  318.         *(arg + i) = (float)atan2(*(y + i), *(x + i));
  319.     }
  320. }
  321.  
  322. /*
  323.  * msdescriptors()
  324.  *   DESCRIPTION:
  325.  *     wrapper for call to descriptors()
  326.  *   ARGUMENTS:
  327.  *     delta_phik: delta phik for polygon (float *)
  328.  *     delta_lk: delta lk for polygon (float *)
  329.  *     mcp: length of dphik and dlk (long)
  330.  *     a_n: Fourier coefficients array (float *)
  331.  *     b_n: Fourier coefficients array (float *)
  332.  *     n_order: order of Fourier descriptors to compute (long)
  333.  *   RETURN VALUE:
  334.  *     none
  335.  */
  336.  
  337. void
  338. msdescriptors(float *delta_phik, float *delta_lk, long mcp, float *a_n, float *b_n, long n_order)
  339. {
  340.     int k;
  341.  
  342. #ifdef FD_DEBUG
  343.     printf("\nevaluation of Fourier descriptors (Zahn & Roskies)\n");
  344.     printf("--------------------------------------------------\n");
  345.  
  346.     printf("...given: two vectors, delta_phik[] and delta_lk[]\n");
  347.  
  348.     printf(" delta_phik:\n");
  349.     for(k=0; k<(int)mcp; k++)
  350.     {
  351.         printf( "%7.2f  ",*(delta_phik+k) );
  352.         if( (k+1)%8 == 0)
  353.             printf ("\n");
  354.     }   
  355.     printf ("\n   delta_lk:\n");
  356.     for (k=0; k<(int)mcp; k++)
  357.     {
  358.         printf( "%7.2f  ",*(delta_lk+k) );
  359.         if( (k+1)%8 == 0)
  360.             printf ("\n");
  361.     }
  362. #endif
  363.  
  364.     if(n_order > mcp)
  365.     {
  366.         printf("...polygon only has %ld vertices;");
  367.         printf(" cannot yield %ld harmonics!\n", mcp, n_order);
  368.         exit(1);
  369.     }
  370.     descriptors(delta_phik, delta_lk, mcp, a_n, b_n, n_order);
  371.  
  372. }
  373.