Two dimensional gaussian image generator in C

Question

This is a follow-up question for Image Processing Median Filter in C. I am attempting to create a two dimensional gaussian image like below in C.

The formula is as follows.

The experimental implementation

• base.h: Contains the basic type implementation

  /* Develop by Jimmy Hu */
  
  #ifndef BASE_H
  #define BASE_H
  
  #include <math.h>
  #include <stdbool.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  #include <unistd.h>
  
  #define MAX_PATH 256
  #define FILE_ROOT_PATH "./"
  
  #define True true
  #define False false
  
  typedef struct RGB
  {
      unsigned char channels[3];
  } RGB;
  
  typedef struct HSV
  {
      long double channels[3];    //  Range: 0 <= H < 360, 0 <= S <= 1, 0 <= V <= 255
  }HSV;
  
  typedef struct BMPIMAGE
  {
      char FILENAME[MAX_PATH];
  
      unsigned int XSIZE;
      unsigned int YSIZE;
      unsigned char FILLINGBYTE;
      unsigned char *IMAGE_DATA;
  } BMPIMAGE;
  
  typedef struct RGBIMAGE
  {
      unsigned int XSIZE;
      unsigned int YSIZE;
      RGB *IMAGE_DATA;
  } RGBIMAGE;
  
  typedef struct HSVIMAGE
  {
      unsigned int XSIZE;
      unsigned int YSIZE;
      HSV *IMAGE_DATA;
  } HSVIMAGE;
  
  #endif
  

• imageio.h contains the function declarations which deal with the bmp image read / write operations

  /* Develop by Jimmy Hu */
  
  #ifndef IMAGEIO_H
  #define IMAGEIO_H
  
  #include "base.h"
  
  RGB *raw_image_to_array(const int xsize, const int ysize, const unsigned char *image);
  
  unsigned long bmp_read_x_size(const char *filename, const bool extension);
  
  unsigned long bmp_read_y_size(const char *filename, const bool extension);
  
  char bmp_read(unsigned char *image, const int xsize, const int ysize, const char *filename, const bool extension);
  
  BMPIMAGE bmp_file_read(const char *filename, const bool extension);
  
  int bmp_write(const char *filename, const int xsize, const int ysize, const unsigned char *image);
  
  unsigned char *array_to_raw_image(const int xsize, const int ysize, const RGB* input_data);
  
  unsigned char bmp_filling_byte_calc(const unsigned int xsize);
  
  #endif
  

• imageio.c contains the function definitions which deal with the bmp image read / write operations

  /* Develop by Jimmy Hu */
  
  #include "imageio.h"
  
  RGB *raw_image_to_array(const int xsize, const int ysize, const unsigned char * const image)
  {
      RGB *output;
      output = malloc(sizeof *output * xsize * ysize);
      if(output == NULL)
      {    
          printf("Memory allocation error!");
          return NULL;
      }        
      unsigned char FillingByte;
      FillingByte = bmp_filling_byte_calc(xsize);
      for(int y = 0; y < ysize; y++)
      {
          for(int x = 0; x < xsize; x++)
          {
              for (int color = 0; color < 3; color++) {
                  output[y * xsize + x].channels[color] =
                      image[3 * (y * xsize + x) + y * FillingByte + (2 - color)];
              }
          }
      }
      return output;
  }
  
  //----bmp_read_x_size function----
  unsigned long bmp_read_x_size(const char *filename, const bool extension)
  {
      char fname_bmp[MAX_PATH];
      if(extension == false)
      {    
          sprintf(fname_bmp, "%s.bmp", filename);
      }        
      else
      {    
          strcpy(fname_bmp,filename);
      }    
      FILE *fp;
      fp = fopen(fname_bmp, "rb");
      if (fp == NULL) 
      {     
          printf("Fail to read file!\n");
          return -1;
      }             
      unsigned char header[54];
      fread(header, sizeof(unsigned char), 54, fp);
      unsigned long output;
      output = header[18] + 
          ((unsigned long)header[19] << 8) +
          ((unsigned long)header[20] << 16) +
          ((unsigned long)header[21] << 24);
      fclose(fp);
      return output;
  }
  
  //---- bmp_read_y_size function ----
  unsigned long bmp_read_y_size(const char * const filename, const bool extension)
  {
      char fname_bmp[MAX_PATH];
      if(extension == false)
      {    
          sprintf(fname_bmp, "%s.bmp", filename);
      }        
      else
      {    
          strcpy(fname_bmp,filename);
      }    
      FILE *fp;
      fp = fopen(fname_bmp, "rb");
      if (fp == NULL)
      {
          printf("Fail to read file!\n");
          return -1;
      }             
      unsigned char header[54];
      fread(header, sizeof(unsigned char), 54, fp);
      unsigned long output; 
      output= header[22] + 
          ((unsigned long)header[23] << 8) +
          ((unsigned long)header[24] << 16) +
          ((unsigned long)header[25] << 24);
      fclose(fp);
      return output;
  }
  
  //---- bmp_file_read function ---- 
  char bmp_read(unsigned char * const image, const int xsize, const int ysize, const char * const filename, const bool extension)
  {
      char fname_bmp[MAX_PATH];
      if(extension == false)
      {    
          sprintf(fname_bmp, "%s.bmp", filename);
      }        
      else
      {    
          strcpy(fname_bmp,filename);
      }    
      unsigned char filling_bytes;
      filling_bytes = bmp_filling_byte_calc(xsize);
      FILE *fp;
      fp = fopen(fname_bmp, "rb");
      if (fp==NULL)
      {     
          printf("Fail to read file!\n");
          return -1;
      }             
      unsigned char header[54];
      fread(header, sizeof(unsigned char), 54, fp);
      fread(image, sizeof(unsigned char), (size_t)(long)(xsize * 3 + filling_bytes)*ysize, fp);
      fclose(fp); 
      return 0;
  }
  
  BMPIMAGE bmp_file_read(const char * const filename, const bool extension)
  {
      BMPIMAGE output;
      stpcpy(output.FILENAME, "");
      output.XSIZE = 0;
      output.YSIZE = 0;
      output.IMAGE_DATA = NULL;
      if(filename == NULL)
      {    
          printf("Path is null\n");
          return output;
      }
      char fname_bmp[MAX_PATH];
      if(extension == false)
      {
          sprintf(fname_bmp, "%s.bmp", filename);
      }
      else
      {    
          strcpy(fname_bmp,filename);
      }    
      FILE *fp;
      fp = fopen(fname_bmp, "rb");
      if (fp == NULL)
      {     
          printf("Fail to read file!\n");
          return output;
      }             
      stpcpy(output.FILENAME, fname_bmp);
      output.XSIZE = (unsigned int)bmp_read_x_size(output.FILENAME,true);
      output.YSIZE = (unsigned int)bmp_read_y_size(output.FILENAME,true);
      if( (output.XSIZE == -1) || (output.YSIZE == -1) )
      {     
          printf("Fail to fetch information of image!");
          return output;
      }        
      else
      {     
          printf("Width of the input image: %d\n",output.XSIZE);
          printf("Height of the input image: %d\n",output.YSIZE);
          printf("Size of the input image(Byte): %d\n", (size_t)output.XSIZE * output.YSIZE * 3);
          output.FILLINGBYTE = bmp_filling_byte_calc(output.XSIZE);
          output.IMAGE_DATA = malloc(sizeof *output.IMAGE_DATA * (output.XSIZE * 3 + output.FILLINGBYTE) * output.YSIZE);
          if (output.IMAGE_DATA == NULL)
          { 
              printf("Memory allocation error!");
              return output;
          }     
          else
          {
              for(int i = 0; i < ((output.XSIZE * 3 + output.FILLINGBYTE) * output.YSIZE);i++)
              {
                  output.IMAGE_DATA[i] = 255;
              }
              bmp_read(output.IMAGE_DATA, output.XSIZE, output.YSIZE, output.FILENAME, true);
          }
      }    
      return output;
  }
  
  //----bmp_write function---- 
  int bmp_write(const char * const filename, const int xsize, const int ysize, const unsigned char * const image) 
  {
      unsigned char FillingByte;
      FillingByte = bmp_filling_byte_calc(xsize);
      unsigned char header[54] =
      {
      0x42, 0x4d, 0, 0, 0, 0, 0, 0, 0, 0,
      54, 0, 0, 0, 40, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 24, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0                        
      };                                
      unsigned long file_size = (long)xsize * (long)ysize * 3 + 54;
      unsigned long width, height;
      char fname_bmp[MAX_PATH];
      header[2] = (unsigned char)(file_size &0x000000ff);
      header[3] = (file_size >> 8) & 0x000000ff;        
      header[4] = (file_size >> 16) & 0x000000ff;        
      header[5] = (file_size >> 24) & 0x000000ff;        
  
      width = xsize;
      header[18] = width & 0x000000ff;
      header[19] = (width >> 8) &0x000000ff;
      header[20] = (width >> 16) &0x000000ff;
      header[21] = (width >> 24) &0x000000ff;
  
      height = ysize;
      header[22] = height &0x000000ff;
      header[23] = (height >> 8) &0x000000ff;
      header[24] = (height >> 16) &0x000000ff;
      header[25] = (height >> 24) &0x000000ff;
      sprintf(fname_bmp, "%s.bmp", filename);
      FILE *fp; 
      if (!(fp = fopen(fname_bmp, "wb")))
      {    
          return -1;
      }        
      fwrite(header, sizeof(unsigned char), 54, fp);
      fwrite(image, sizeof(unsigned char), (size_t)(long)(xsize * 3 + FillingByte) * ysize, fp);
      fclose(fp);
      return 0;
  }
  
  unsigned char *array_to_raw_image(const int xsize, const int ysize, const RGB* const input_data)
  {
      unsigned char FillingByte;
      FillingByte = bmp_filling_byte_calc(xsize);
      unsigned char *output;
      output = malloc(sizeof *output * (xsize * 3 + FillingByte) * ysize);
      if(output == NULL)
      {    
          printf("Memory allocation error!");
          return NULL;
      }
      for(int y = 0;y < ysize;y++)
      {
          for(int x = 0;x < (xsize * 3 + FillingByte);x++)
          {
              output[y * (xsize * 3 + FillingByte) + x] = 0;
          }
      }
      for(int y = 0;y<ysize;y++)
      {
          for(int x = 0;x<xsize;x++)
          {
              for (int color = 0; color < 3; color++) {
                  output[3 * (y * xsize + x) + y * FillingByte + (2 - color)]
                  = input_data[y*xsize + x].channels[color];
              }
          }
      }
      return output;
  }
  
  unsigned char bmp_filling_byte_calc(const unsigned int xsize)
  {
      unsigned char filling_bytes;
      filling_bytes = ( xsize % 4);
      return filling_bytes;
  }
  

• main.c

  /* Develop by Jimmy Hu */
  
  #include "base.h"
  #include "imageio.h"
  
  RGB* GaussianFigure2D(const int xsize, const int ysize, const long double standard_deviation);
  
  long double NormalDistribution2D(long double xlocation, long double ylocation, long double standard_deviation);
  
  int main(int argc, char** argv)
  {
      printf("BMP image file name:(ex:test): ");
      char *FilenameString;
      FilenameString = malloc( sizeof *FilenameString * MAX_PATH);
      scanf("%s",FilenameString);
      int xsize = 512, ysize = 512;
      RGBIMAGE RGBImage1;
      RGBImage1.XSIZE = xsize;
      RGBImage1.YSIZE = ysize;
      RGBImage1.IMAGE_DATA = GaussianFigure2D(RGBImage1.XSIZE, RGBImage1.YSIZE, 100);
      bmp_write(FilenameString, RGBImage1.XSIZE, RGBImage1.YSIZE, array_to_raw_image(RGBImage1.XSIZE, RGBImage1.YSIZE, RGBImage1.IMAGE_DATA));
      free(FilenameString);
      return 0;
  }
  
  RGB* GaussianFigure2D(const int xsize, const int ysize, const long double standard_deviation)
  {
      RGB* output;
      output = malloc(sizeof *output * xsize * ysize);
      if (output == NULL)
      {
          printf("Memory allocation error!");
          return NULL;
      }
      int Centerx = xsize / 2, Centery = ysize / 2;
      long double normalize_factor;
      normalize_factor = (long double) 127.0 / NormalDistribution2D(0, 0, standard_deviation);
  
      for(int x = 0; x < xsize; x ++)
      {
          for(int y = 0; y < ysize; y ++)
          {
              long double Data = 
              NormalDistribution2D(fabs(Centerx - x), fabs(Centery - y), standard_deviation) * normalize_factor + 128;
              for (int channel_index = 0; channel_index < 3; channel_index++) {
                  output[ y * xsize + x ].channels[channel_index] = Data;
              }
          }
      }
      return output;
  }
  
  long double NormalDistribution2D(long double xlocation, long double ylocation, long double standard_deviation)
  {
      return pow(M_E, -(pow(xlocation,2) + pow(ylocation,2)) / (2 * pow(standard_deviation, 2))) / (2 * M_PI * pow(standard_deviation, 2));
  }
  

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Image Processing Median Filter in C

• What changes has been made in the code since last question?

  I am attempting to create a two dimensional gaussian image in C in this post.

• Why a new review is being asked for?

  If there is any possible improvement, please let me know.

Reference

Multivariate Normal Distribution

Answer 1

Don't use long double

Almost no code needs long double; it has much more precision than you need, and it is much less efficient than a float. Really, stop using it.

Use size_t for sizes, counts and indices

This is something that has been mentioned multiple times in answers to earlier questions posted by you. An int is not guaranteed to be large enough to represent all sizes of things that fit in memory. The proper type is size_t.

Print errors to stderr

Error messages should be printed to stderr, not to stdout. This prevents them from being mixed with the regular output (which is especially important if you are redirecting standard output to a file or a pipe).

Guard against buffer overflows

Your code contains multiple possible buffer overflows. For example, in bmp_read_x_size(), if the input filename is longer than MAX_PATH, you will write past the end of fname_bmp. Use functions that guard against buffer overflows, like snprintf(). Also, even if those functions might prevent the buffer overflow by truncating the output, consider if that truncation might be dangerous in itself.

Check the return value of all calls to I/O functions

It's not enough to check if fopen() failed or not. Errors can happen even after opening a file. Make sure you also check the return value of fread(), fwrite() and fclose(), and return an error in case of any error.

Wrong order of nested loop

In GuassianFigure2D, the outer loop loops over x, the inner one over y. This results in a less efficient memory access pattern. Switch the loops to ensure memory is accessed sequentially.

Possible optimizations

Compilers might not always be able to optimize inefficient expressions when floating point math is involved. Try to write them as efficient as possible yourself. Use expf() instead of pow(M_E, ...), write x * x instead of pow(x, 2), and try to move constants out of loops where possible. For example, instead of NormalDistribution2D() having to divide by \$2\pi\sigma^2\$, move that into normalize_factor instead.

You also don't need to call fabs() when passing coordinates to NormalDistribution2D(), the coordinates are squared anyway.

Another optimization is realizing that there are four quadrants in the output that are just mirror copies of each other. Use this to reduce the number of times you need to calculate the Gaussian by a factor four.

Loose ends

In main() you take care of freeing FilenameString, but you forget to free RGBImage1. Make it a habit to clean up properly, even if you know the program will exit immediately afterwards.

Consider dithering the output

There are only 128 possible intensity levels in your output. However, the human eye is quite sensitive to small changes in intensity, and your output has faint but noticable color banding. The usual trick is to apply dithering to the output. There are various ways to do dithering, but even a very simple 1D error diffusion algorithm would already be a nice improvement:

for (...) {
    float error = 0;
    for (...) {
        float Data = NormalDistribution2D(...) * normalize_factor + 128;
        Data += error;
        for (...) {
            output[y * xsize + x].channels[channel_index] = Data;
            error = Data - output[y * xsize + x].channels[channel_index];
        }
    }
}

Answer 2

Note that the multivariate Gaussian function is separable: G(x,y) = G(x) G(y). You can use this to significantly reduce the computational cost of your generator:

1. Create an array G of size max(sizex, sizey) / 2 + 1, and fill it with the output of NormalDistribution1D(index, standard_deviation) (where index is the array index).

2. Iterate over your image, for each pixel (x,y) write

   lround(G[abs(Centerx - x)] * G[abs(Centery - y)] * 127) + 128
   

NormalDistribution1D() should return a non-normalized Gaussian (i.e. a value in the range [0,1]). Your current NormalDistribution2D() divides by 2 * M_PI * pow(standard_deviation, 2), which you then have to undo by multiplying by normalize_factor. The non-normalized 1D Gaussian is simply

exp(-x * x / (2 * s * s))

with x and s your xlocation and standard_deviation respectively. Or use expf if you plan to use single-precision floating-point computations as suggested in the other answer, I don't think that improves computation time on modern hardware (see comments below).