Two dimensional bicubic interpolation implementation in C

Question

This is a follow-up question for Two dimensional bicubic interpolation implementation in Matlab and Two dimensional gaussian image generator in C. Besides the Matlab version code, I am attempting to make a C version two dimensional bicubic interpolation function BicubicInterpolation here.

The experimental implementation

• BicubicInterpolation function implementation:

  RGB* BicubicInterpolation(const RGB* const image, const int originSizeX, const int originSizeY, const int newSizeX, const int newSizeY)
  {
      RGB* output;
      output = malloc(sizeof *output * newSizeX * newSizeY);
      if (output == NULL)
      {
          printf(stderr, "Memory allocation error!");
          return NULL;
      }
  
      float ratiox = (float)originSizeX / (float)newSizeX;
      float ratioy = (float)originSizeY / (float)newSizeY;
  
      for (size_t y = 0; y < newSizeY; y++)
      {
          for (size_t x = 0; x < newSizeX; x++)
          {
              for (size_t channel_index = 0; channel_index < 3; channel_index++) {
                  float xMappingToOrigin = (float)x * ratiox;
                  float yMappingToOrigin = (float)y * ratioy;
                  float xMappingToOriginFloor = floor(xMappingToOrigin);
                  float yMappingToOriginFloor = floor(yMappingToOrigin);
                  float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                  float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
  
                  unsigned char* ndata;
                  ndata = malloc(sizeof *ndata * 4 * 4);
                  if (ndata == NULL)
                  {
                      printf(stderr, "Memory allocation error!");
                      return NULL;
                  }
                  for (int ndatay = -1; ndatay < 2; ndatay++)
                  {
                      for (int ndatax = -1; ndatax < 2; ndatax++)
                      {
                          ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image[
                              clip(yMappingToOriginFloor + ndatay, 0, originSizeY - 1) * originSizeX + 
                              clip(xMappingToOriginFloor + ndatax, 0, originSizeX - 1)
                              ].channels[channel_index];
                      }
  
                  }
  
                  unsigned char result = BicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
                  output[ y * newSizeX + x ].channels[channel_index] = result;
                  free(ndata);
              }
          }
      }
      return output;
  }
  

• The other used functions:

  unsigned char BicubicPolate(const unsigned char* const ndata, const float fracx, const float fracy)
  {
      float x1 = CubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
      float x2 = CubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
      float x3 = CubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
      float x4 = CubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );
  
      float output = clip_float(CubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
      return (unsigned char)output;
  }
  
  float CubicPolate(const float v0, const float v1, const float v2, const float v3, const float fracy)
  {
      float A = (v3-v2)-(v0-v1);
      float B = (v0-v1)-A;
      float C = v2-v0;
      float D = v1;
      return D + fracy * (C + fracy * (B + fracy * A));
  }
  
  size_t clip(const size_t input, const size_t lowerbound, const size_t upperbound)
  {
      if (input < lowerbound)
      {
          return lowerbound;
      }
      if (input > upperbound)
      {
          return upperbound;
      }
      return input;
  }
  
  float clip_float(const float input, const float lowerbound, const float upperbound)
  {
      if (input < lowerbound)
      {
          return lowerbound;
      }
      if (input > upperbound)
      {
          return upperbound;
      }
      return input;
  }
  

• base.h

  /* 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
  

The full testing code

/* Develop by Jimmy Hu */

#include "base.h"
#include "imageio.h"

RGB* BicubicInterpolation(const RGB* const image, const int originSizeX, const int originSizeY, const int newSizeX, const int newSizeY);

unsigned char BicubicPolate(const unsigned char* ndata, const float fracx, const float fracy);

float CubicPolate(const float v0, const float v1, const float v2, const float v3, const float fracy);

size_t clip(const size_t input, const size_t lowerbound, const size_t upperbound);

float clip_float(const float input, const float lowerbound, const float upperbound);

int main(int argc, char** argv)
{
    char *FilenameString;
    FilenameString = malloc( sizeof *FilenameString * MAX_PATH);
    
    printf("BMP image input file name:(ex:test): ");
    scanf("%s", FilenameString);
    BMPIMAGE BMPImage1 = bmp_file_read(FilenameString, false);
    RGBIMAGE RGBImage1;
    RGBImage1.XSIZE = BMPImage1.XSIZE;
    RGBImage1.YSIZE = BMPImage1.YSIZE;
    RGBImage1.IMAGE_DATA = raw_image_to_array(BMPImage1.XSIZE, BMPImage1.YSIZE, BMPImage1.IMAGE_DATA);

    RGBIMAGE RGBImage2;
    RGBImage2.XSIZE = 1024;
    RGBImage2.YSIZE = 1024;
    RGBImage2.IMAGE_DATA = BicubicInterpolation(RGBImage1.IMAGE_DATA, RGBImage1.XSIZE, RGBImage1.YSIZE, RGBImage2.XSIZE, RGBImage2.YSIZE);
    
    printf("file name for saving:(ex:test): ");
    scanf("%s", FilenameString);
    bmp_write(FilenameString, RGBImage2.XSIZE, RGBImage2.YSIZE, array_to_raw_image(RGBImage2.XSIZE, RGBImage2.YSIZE, RGBImage2.IMAGE_DATA));

    free(FilenameString);
    free(RGBImage1.IMAGE_DATA);
    free(RGBImage2.IMAGE_DATA);
    return 0;
}

RGB* BicubicInterpolation(const RGB* const image, const int originSizeX, const int originSizeY, const int newSizeX, const int newSizeY)
{
    RGB* output;
    output = malloc(sizeof *output * newSizeX * newSizeY);
    if (output == NULL)
    {
        printf(stderr, "Memory allocation error!");
        return NULL;
    }
    
    float ratiox = (float)originSizeX / (float)newSizeX;
    float ratioy = (float)originSizeY / (float)newSizeY;
    
    for (size_t y = 0; y < newSizeY; y++)
    {
        for (size_t x = 0; x < newSizeX; x++)
        {
            for (size_t channel_index = 0; channel_index < 3; channel_index++) {
                float xMappingToOrigin = (float)x * ratiox;
                float yMappingToOrigin = (float)y * ratioy;
                float xMappingToOriginFloor = floor(xMappingToOrigin);
                float yMappingToOriginFloor = floor(yMappingToOrigin);
                float xMappingToOriginFrac = xMappingToOrigin - xMappingToOriginFloor;
                float yMappingToOriginFrac = yMappingToOrigin - yMappingToOriginFloor;
                
                unsigned char* ndata;
                ndata = malloc(sizeof *ndata * 4 * 4);
                if (ndata == NULL)
                {
                    printf(stderr, "Memory allocation error!");
                    return NULL;
                }
                for (int ndatay = -1; ndatay < 2; ndatay++)
                {
                    for (int ndatax = -1; ndatax < 2; ndatax++)
                    {
                        ndata[(ndatay + 1) * 4 + (ndatax + 1)] = image[
                            clip(yMappingToOriginFloor + ndatay, 0, originSizeY - 1) * originSizeX + 
                            clip(xMappingToOriginFloor + ndatax, 0, originSizeX - 1)
                            ].channels[channel_index];
                    }
                    
                }

                unsigned char result = BicubicPolate(ndata, xMappingToOriginFrac, yMappingToOriginFrac);
                output[ y * newSizeX + x ].channels[channel_index] = result;
                free(ndata);
            }
        }
    }
    return output;
}

unsigned char BicubicPolate(const unsigned char* const ndata, const float fracx, const float fracy)
{
    float x1 = CubicPolate( ndata[0], ndata[1], ndata[2], ndata[3], fracx );
    float x2 = CubicPolate( ndata[4], ndata[5], ndata[6], ndata[7], fracx );
    float x3 = CubicPolate( ndata[8], ndata[9], ndata[10], ndata[11], fracx );
    float x4 = CubicPolate( ndata[12], ndata[13], ndata[14], ndata[15], fracx );

    float output = clip_float(CubicPolate( x1, x2, x3, x4, fracy ), 0.0, 255.0);
    return (unsigned char)output;
}

float CubicPolate(const float v0, const float v1, const float v2, const float v3, const float fracy)
{
    float A = (v3-v2)-(v0-v1);
    float B = (v0-v1)-A;
    float C = v2-v0;
    float D = v1;
    return D + fracy * (C + fracy * (B + fracy * A));
}

size_t clip(const size_t input, const size_t lowerbound, const size_t upperbound)
{
    if (input < lowerbound)
    {
        return lowerbound;
    }
    if (input > upperbound)
    {
        return upperbound;
    }
    return input;
}

float clip_float(const float input, const float lowerbound, const float upperbound)
{
    if (input < lowerbound)
    {
        return lowerbound;
    }
    if (input > upperbound)
    {
        return upperbound;
    }
    return input;
}

All suggestions are welcome.

The summary information:

• Which question it is a follow-up to?

  Two dimensional bicubic interpolation implementation in Matlab and

  Two dimensional gaussian image generator in C.

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

  I am attempting to make a C version two dimensional bicubic interpolation function in this post.

• Why a new review is being asked for?

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

Answer 1

Avoid unnecesary allocation of temporary storage

In the innermost loop, you do this:

unsigned char* ndata;
ndata = malloc(sizeof *ndata * 4 * 4);

This is slow and completely unnecessary; you can just declare an array on the stack like so:

unsigned char ndata[4 * 4];

Possible improvements to the algorithm

It is likely that many of the intermediate values you are calculating in BicubicPolate() might be the same as those for neighbouring pixels. Also in CubicPolate(), none of the values of A to D depend on fracy, and some preprocessing of the image might allow you to avoid many of the operations.

Also consider that the ratio between the source and destination can be larger than 1 or smaller than 1, and different algorithms might be better for each case, and ratios of the form n or 1 / n, where n is an integer, might especially be candidates for algorithmic improvements.

Answer 2

Ok let's try it, using this test image. Frame 0 was resized (from a width of 512 to a width of 300) by other software, frame 1 by this code.

That regular pattern of darker pixels is not supposed to appear.

It looks like BicubicPolate reads some entries from ndata (16 bytes) that were never written to (9 bytes are written to it).

I'm not exactly sure how that is supposed to work, but changing the loops that fill ndata to go up to and including 2 seems to improve the output.

Answer 3

I know I might be a little late on this one. But I am currently working on a similar algorithm and was greatly inspired by yours. One thing I found surprising, is that you recalculate the values A B C and D for every pixel. I think they should be the same for pixels with the same x value, that lie in the same "square" of original data. You could probably optimize on this.