gauss.c Source File

00001 #include "sadie.h"
00002 #include "proto.h"
00003 
00004 /*----------------------------------------------------------------------------*/
00005 /*-General Information--------------------------------------------------------*/
00006 /*                                                                            */
00007 /* This function computes a (size x size) Gaussian filter with standard   */
00008 /* deviation sigma.                                                           */
00009 /*                                                                            */
00010 /*----------------------------------------------------------------------------*/
00011 /*-Interface Information------------------------------------------------------*/
00012 void CREATEGAUSS (
00013 int size,       /*  I   The size of the image.                                */
00014 double sigma,   /*  I   The standard deviation of the gaussian.               */
00015 IMAGE **out     /*  O   Address of a pointer to the output image              */
00016 /*----------------------------------------------------------------------------*/
00017 ) { register short i, j, k;
00018     char msg[SLEN];
00019     double x, y, *elin=0, *epix=0;
00020 
00021     /* if (TIMES) TIMING(T_GAUSS); */
00022     if (NAMES) {
00023         MESSAGE('I',"");
00024         MESSAGE('I',"GAUSS");
00025         MESSAGE('I',"");
00026         sprintf(msg," Size:                       %i x %i",size,size);
00027         MESSAGE('I',msg);
00028         sprintf(msg," Std. Deviation:             %f",sigma);
00029         MESSAGE('I',msg);
00030         MESSAGE('I'," ...............");
00031     }
00032 
00033     /* check parameters */
00034     if (size < 2  ||  size%2 == 0) {
00035         MESSAGE('E'," Size of gaussian must be an odd number greater than two.");
00036         goto the_end;
00037     } else if (sigma <= 0.0) {
00038         MESSAGE('E'," Standard deviation of gaussian must be positive.");
00039         goto the_end;
00040     }
00041     
00042     /* create image of appropriate size */
00043     if (!CHECKIMG(*out)) GETMEM(1,size,size,out);
00044     if (!*out) goto the_end;
00045 
00046     /* allocate memory for temporary storage */
00047     if (!(elin = (double *)malloc(((size/2)+1)*sizeof(double)))) {
00048         MESSAGE('E'," Memory request failed.");
00049         goto the_end;
00050     }
00051     if (!(epix = (double *)malloc(((size/2)+1)*sizeof(double)))) {
00052         MESSAGE('E'," Memory request failed.");
00053         goto the_end;
00054     }
00055 
00056 
00057     /* store one dimensional components */
00058     for (j=0; j<size/2; j++) {
00059         y = ((double)(j-size/2)*(double)(j-size/2))/(2*sigma*sigma);
00060         elin[j] = exp(-y);
00061     }
00062     for (k=0; k<size/2; k++) {
00063         x = ((double)(k-size/2)*(double)(k-size/2))/(2*sigma*sigma);
00064         epix[k] = exp(-x);
00065     }
00066     /* symmetric quadrants */
00067     for (j=0; j<size/2; j++) {
00068         for (k=0; k<size/2; k++) {
00069             (*out)->data[0][j][k] =
00070             (*out)->data[0][j][size-k-1] =
00071             (*out)->data[0][size-j-1][k] =
00072             (*out)->data[0][size-j-1][size-k-1] =
00073             (PIXEL)(epix[k]*elin[j]);
00074         }
00075     }
00076     /* zero coordinates */
00077     for (j=0; j<size; j++) {
00078         y = ((double)(j-size/2)*(double)(j-size/2))/(2*sigma*sigma);
00079         (*out)->data[0][j][size/2] = (PIXEL)(exp(-y));
00080         (*out)->data[0][size/2][j] = (PIXEL)(exp(-y));
00081     }
00082 
00083         
00084     the_end:
00085     if (elin)  free(elin);
00086     if (epix)  free(epix);
00087     /* if (TIMES) TIMING(T_EXIT); */
00088 }