home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / Software / ASCENDER / ascendMar8.tar / UMass / BoldtNew / LLVS / signal2D.c < prev    next >
C/C++ Source or Header  |  1996-01-31  |  9KB  |  368 lines
  1. /***********************************************************************
  2. *  Signal2D.c - convolution with a 2d template
  3. *  Author: Xiao-Guang Wang
  4. *  Modifications:
  5. *    2/20/95 - Bob Collins, added new convolution functions for
  6. *      different input and output types.
  7. ************************************************************************/
  8.  
  9.  
  10. #include <stdio.h>
  11. #include <stdlib.h>
  12. #include "signal2D.h"
  13.  
  14. void free_signal_2D(Signal_2D * psignal)
  15. {
  16.  if((psignal->xwidth > 0) && (psignal->xwidth > 0)) free(psignal->value);
  17.  psignal->xwidth = 0;
  18.  psignal->ywidth = 0;
  19.  psignal->xbase = 0;
  20.  psignal->ybase = 0;
  21. }
  22.  
  23.  
  24. Bool read_signal_2D(char * filename, Signal_2D * psignal)
  25. /*
  26. Read the 2D signal from a file.
  27. */
  28. {
  29.  int j, i;
  30.  FILE * fp;
  31.  double * dpvalue;
  32.  
  33.  fp = fopen(filename, "r");
  34.  if(! fp)
  35.  {
  36.   printf("template file open error.\n");
  37.   return(FALSE);
  38.  }
  39.  
  40.  fscanf(fp, "%d  %d", (& psignal->xwidth), (& psignal->ywidth));
  41.  psignal->xbase = psignal->xwidth / 2;
  42.  psignal->ybase = psignal->ywidth / 2;
  43.  psignal->value = (double *)malloc(  psignal->xwidth * psignal->ywidth
  44.                                    * sizeof(double));
  45.  dpvalue = psignal->value;
  46.  for(i = 0; i < psignal->ywidth; i ++)
  47.  {
  48.   for(j = 0; j < psignal->xwidth; j ++)
  49.   {
  50.    fscanf(fp, "%lf", dpvalue);
  51.    dpvalue ++;
  52.   }
  53.  }
  54.  fclose(fp);
  55.  return(TRUE);
  56. }
  57.  
  58.  
  59. Bool write_signal_2D(char * filename, Signal_2D * psignal)
  60. /*
  61. Write the 2D signal into a file.
  62. */
  63. {
  64.  int j, i;
  65.  FILE * fp;
  66.  double * dpvalue;
  67.  
  68.  fp = fopen(filename, "w");
  69.  if(! fp)
  70.  {
  71.   printf("template file open error.\n");
  72.   return(FALSE);
  73.  }
  74.  
  75.  dpvalue = psignal->value;
  76.  for(i = 0; i < psignal->ywidth; i ++)
  77.  {
  78.   for(j = 0; j < psignal->xwidth; j ++)
  79.   {
  80.    dpvalue ++;
  81.   }
  82.   fprintf(fp, "\n");
  83.  }
  84.  fclose(fp);
  85.  return(TRUE);
  86. }
  87.  
  88.  
  89. #define Signal_value_significance 0.1
  90. /*
  91. void get_signal_2D_gaussian(double sd, Signal_2D * pgaussian)
  92. {
  93.  int i;
  94.  double * list;
  95.  int list_length = 0;
  96.  double numerator, denominator;
  97.  double next_value;
  98.  double sum;
  99.  
  100. /the middle value/
  101.  append_list((char **) & list, & list_length, sizeof(double));
  102.  list[0] = 1.0;
  103.  
  104. /if sd is too small, the gaussian becomes an impulse/
  105.  if(sd > Small_double)
  106.  {
  107.   denominator = 2.0 * sd * sd;
  108.   numerator = (double)list_length * list_length;
  109.   next_value = exp(- (numerator / denominator));
  110.  
  111. /other values/
  112.   while(next_value > (list[0] * Signal_value_significance))
  113.   {
  114.    append_list((char **) & list, & list_length, sizeof(double));
  115.    list[list_length - 1] = next_value;
  116.    numerator = (double)list_length * list_length;
  117.    next_value = exp(- (numerator / denominator));
  118.   }
  119.  }
  120.  
  121. /write to the template/
  122.  pgaussian->width = (2 * list_length) - 1;
  123.  pgaussian->base = list_length - 1;
  124.  pgaussian->value = (double *)malloc(pgaussian->width * sizeof(double));
  125.  
  126.  pgaussian->value[(pgaussian->width - 1) / 2] = list[0];
  127.  for(i = 1; i < (pgaussian->width + 1) / 2; i ++)
  128.  {
  129.   pgaussian->value[((pgaussian->width - 1) / 2) - i] = list[i];
  130.   pgaussian->value[((pgaussian->width - 1) / 2) + i] = list[i];
  131.  }
  132.  
  133.  free_list((char **) & list, & list_length);
  134.  
  135. /normalize/
  136.  sum = 0;
  137.  for(i = 0; i < pgaussian->width; i ++)
  138.  {
  139.   sum += pgaussian->value[i];
  140.  }
  141.  for(i = 0; i < pgaussian->width; i ++)
  142.  {
  143.   pgaussian->value[i] /= sum;
  144.  }
  145. }
  146. */
  147.  
  148.  
  149. Bool signal_2D_convolve_edgecut(short * psignal, int xsize, int ysize,
  150.                                 Signal_2D * ptemplate, short ** ppresult)
  151. /*
  152. This function convolve the psignal with the ptemplate, and write the
  153. result to presult. For image borders, this function cuts the edge. That is,
  154. only part of the template (the part that is in the signal) are multiplied
  155. with the signal. presult has the same size as psignal, and is allocated in
  156. this function. Note: the template is not inverted in this convolution.
  157. */
  158. {
  159.  int l, k, j, i;
  160.  int n, m;
  161.  double sum;
  162.  short * dpsignal;
  163.  double * dptempvalue;
  164.  short * presult;
  165.  short * dpresult;
  166.  
  167.  presult = (short *)malloc(xsize * ysize * sizeof(short));
  168.  if(presult == NULL)
  169.  {
  170.   printf("presult allocation error in convolution.\n");
  171.   return(FALSE);
  172.  }
  173.  
  174.  dpresult = presult;
  175.  for(i = 0; i < ysize; i ++)
  176.  {
  177.   for(j = 0; j < xsize; j ++)
  178.   {
  179.    sum = 0;
  180.    for(k = 0; k < ptemplate->ywidth; k ++)
  181.    {
  182.     m = i + (k - ptemplate->ybase);
  183.     if((m >= 0) && (m < ysize))
  184.     {
  185.      dpsignal = psignal + (m * xsize);
  186.      dptempvalue = ptemplate->value + (k * ptemplate->xwidth);
  187.      for(l = 0; l < ptemplate->xwidth; l ++)
  188.      {
  189.       n = j + (l - ptemplate->xbase);
  190.       if((n >= 0) && (n < xsize))
  191.       {
  192.        sum += ((* (dptempvalue + l)) * (* (dpsignal + n)));
  193.       }
  194.      } /*l*/
  195.     }
  196.    } /*k*/
  197.    (* dpresult) = (short)sum;
  198.    dpresult ++;
  199.   } /*j*/
  200.  } /*i*/
  201.  
  202.  (* ppresult) = presult;
  203.  return(TRUE);
  204. }
  205.  
  206.  
  207. Bool shortIn_floatOut_signal_2D_convolve_edgecut(short * psignal, 
  208.       int xsize, int ysize, Signal_2D * ptemplate, float * ppresult)
  209. /*
  210. Same as signal_2D_convolve_edgecut, but produces floating point output.
  211. */
  212. {
  213.  int l, k, j, i;
  214.  int n, m;
  215.  double sum;
  216.  short * dpsignal;
  217.  double * dptempvalue;
  218.  float * dpresult;
  219.  
  220. /* Memory allocation for result is done prior to this function --BobC
  221.  presult = (float *)malloc(xsize * ysize * sizeof(float));
  222.  if(presult == NULL)
  223.  {
  224.   printf("presult allocation error in convolution.\n");
  225.   return(FALSE);
  226.  }
  227. */
  228.  
  229.  dpresult = ppresult;
  230.  for(i = 0; i < ysize; i ++)
  231.  {
  232.   for(j = 0; j < xsize; j ++)
  233.   {
  234.    sum = 0;
  235.    for(k = 0; k < ptemplate->ywidth; k ++)
  236.    {
  237.     m = i + (k - ptemplate->ybase);
  238.     if((m >= 0) && (m < ysize))
  239.     {
  240.      dpsignal = psignal + (m * xsize);
  241.      dptempvalue = ptemplate->value + (k * ptemplate->xwidth);
  242.      for(l = 0; l < ptemplate->xwidth; l ++)
  243.      {
  244.       n = j + (l - ptemplate->xbase);
  245.       if((n >= 0) && (n < xsize))
  246.       {
  247.        sum += ((* (dptempvalue + l)) * (* (dpsignal + n)));
  248.       }
  249.      } /*l*/
  250.     }
  251.    } /*k*/
  252.    (* dpresult) = (float)sum;
  253.    dpresult ++;
  254.   } /*j*/
  255.  } /*i*/
  256.  
  257.  return(TRUE);
  258. }
  259.  
  260. Bool byteIn_floatOut_signal_2D_convolve_edgecut(Byte * psignal, 
  261.       int xsize, int ysize, Signal_2D * ptemplate, float * ppresult)
  262. /*
  263. Same as signal_2D_convolve_edgecut, but byte input and floating point output
  264. */
  265. {
  266.  int l, k, j, i;
  267.  int n, m;
  268.  double sum;
  269.  Byte * dpsignal;
  270.  double * dptempvalue;
  271.  float * dpresult;
  272.  
  273. /* Memory allocation for result is done prior to this function --BobC
  274.  presult = (float *)malloc(xsize * ysize * sizeof(float));
  275.  if(presult == NULL)
  276.  {
  277.   printf("presult allocation error in convolution.\n");
  278.   return(FALSE);
  279.  }
  280. */
  281.  
  282.  dpresult = ppresult;
  283.  for(i = 0; i < ysize; i ++)
  284.  {
  285.   for(j = 0; j < xsize; j ++)
  286.   {
  287.    sum = 0;
  288.    for(k = 0; k < ptemplate->ywidth; k ++)
  289.    {
  290.     m = i + (k - ptemplate->ybase);
  291.     if((m >= 0) && (m < ysize))
  292.     {
  293.      dpsignal = psignal + (m * xsize);
  294.      dptempvalue = ptemplate->value + (k * ptemplate->xwidth);
  295.      for(l = 0; l < ptemplate->xwidth; l ++)
  296.      {
  297.       n = j + (l - ptemplate->xbase);
  298.       if((n >= 0) && (n < xsize))
  299.       {
  300.        sum += ((* (dptempvalue + l)) * (* (dpsignal + n)));
  301.       }
  302.      } /*l*/
  303.     }
  304.    } /*k*/
  305.    (* dpresult) = (float)sum;
  306.    dpresult ++;
  307.   } /*j*/
  308.  } /*i*/
  309.  
  310.  return(TRUE);
  311. }
  312.  
  313.  
  314. Bool signal_2D_convolve_rotation(short * psignal, int xsize, int ysize,
  315.                                  Signal_2D * ptemplate, short ** ppresult)
  316. /*
  317. This function convolve the psignal with the ptemplate, and write the
  318. result to presult. For the two ends, this function uses the rotation
  319. method. psignal must have a length not shorter than ptemplate.
  320. presult has the same width as psignal, and is allocated in this function.
  321. Note: the template is not inverted in this convolution.
  322. */
  323. {
  324. /*
  325.  int j, i;
  326.  
  327.  presult->width = psignal->width;
  328.  presult->base = psignal->base;
  329.  presult->value = (double *)malloc(presult->width * sizeof(double));
  330.  
  331.  for(i = 0; i < presult->width; i ++)
  332.  {
  333.   presult->value[i] = 0;
  334.   for(j = 0; j < ptemplate->width; j ++)
  335.   {
  336.    presult->value[i]
  337.        += (  ptemplate->value[j]
  338.            * psignal->value[mod((i + (j - ptemplate->base)), psignal->width)]
  339.           );
  340.   }
  341.  }
  342. */
  343.  return(TRUE);
  344. }
  345.  
  346.  
  347. void signal_2D_difference(Signal_2D * psignal, Signal_2D * presult)
  348. /*
  349. In this function psignal is differenced and the result is written
  350. into presult. psignal must have a width longer than 1.
  351. presult->value[0] is meaningless.
  352. */
  353. {
  354. /*
  355.  int i;
  356.  
  357.  presult->width = psignal->width;
  358.  presult->base = psignal->base;
  359.  presult->value = (double *)malloc(presult->width * sizeof(double));
  360.  
  361.  presult->value[0] = 0;
  362.  for(i = 1; i < presult->width; i ++)
  363.  {
  364.   presult->value[i] = psignal->value[i] - psignal->value[i - 1];
  365.  }
  366. */
  367. }
  368.