Different results with libpng and numpy for average vector computation

Question

I'm trying to compute the average vector in RGB tristimulus space with libpng in C, and NumPy in Python, but I'm getting different results with each. I'm quite confident Python is giving the correct result with this image of [ 127.5 127.5 0. ]. However, with the following block of C I get the preposterous result of [ 38.406494 38.433670 38.459641 ]. I've been staring at my code for weeks without any give, so I thought I'd see if others may have an idea.

Also, I've tested this code with other images and it gives similar preposterous results. It's quite curious because all three numbers usually match for the first 4 or so digits. I'm not sure what may be causing this.

/* See if our average vector matches that of Python's */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <png.h>

// For getting the PNG data and header/information back
typedef struct 
{
    uint32_t width;         // width of image
    uint32_t height;        // height of image
    int bit_depth;          // bits/pixel component (should be 8 in RGB)
    png_bytep datap;        // data
} rTuple;

#define PNG_BYTES_TO_CHECK 8
#define CHANNELS 3

int
check_PNG_signature(unsigned char *buffer)
{
    unsigned i;
    const unsigned char signature[8] = { 0x89, 0x50, 0x4e, 0x47, 
                                         0x0d, 0x0a, 0x1a, 0x0a };
    for (i = 0; i < PNG_BYTES_TO_CHECK; ++i) 
    {
        if (buffer[i] != signature[i]) 
        {
            fprintf(stderr, "** File sig does not match PNG, received ");
            for (i = 0; i < PNG_BYTES_TO_CHECK; ++i)
                fprintf(stderr, "%.2X ", buffer[i]);
            fprintf(stderr, "\n");
            abort();
        }   
    }
    return 1;
}

rTuple 
read_png_file(char *file_name)
{
    /* Get PNG data - I've pieced this together by reading `example.c` from
       beginning to end */
    printf("** Reading data from %s\n", file_name);

    png_uint_32 width, height;  // holds width and height of image

    uint32_t row;  // for iteration later
    int bit_depth, color_type, interlace_type;

    unsigned char *buff = malloc(PNG_BYTES_TO_CHECK * sizeof(char));
    memset(buff, 0, PNG_BYTES_TO_CHECK * sizeof(char));

    FILE *fp = fopen(file_name, "rb");
    if (fp == NULL) abort();

    if (fread(buff, 1, PNG_BYTES_TO_CHECK, fp) != PNG_BYTES_TO_CHECK) {
        fprintf(stderr, "** Could not read %d bytes\n", PNG_BYTES_TO_CHECK);
        abort();
    }

    check_PNG_signature(buff);
    rewind(fp);

    // create and initialize the png_struct, which will be destroyed later
    png_structp png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING
        , NULL  /* Following 3 mean use stderr & longjump method */
        , NULL
        , NULL
    );
    if (!png_ptr) abort();

    png_infop info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr) abort();

    // following I/O initialization method is required
    png_init_io(png_ptr, fp);
    png_set_sig_bytes(png_ptr, 0);  // libpng has this built in too

    // call to png_read_info() gives us all of the information from the
    // PNG file before the first IDAT (image data chunk)
    png_read_info(png_ptr, info_ptr);

    // Get header metadata now
    png_get_IHDR(png_ptr, info_ptr, &width, &height, &bit_depth, &color_type, 
        &interlace_type, NULL, NULL);

    // Scale 16-bit images to 8-bits as accurately as possible (shouldn't be an
    // issue though, since we're working with RGB data)
#ifdef PNG_READ_SCALE_16_TO_8_SUPPORTED
    png_set_scale_16(png_ptr);
#else
    png_set_strip_16(png_ptr);
#endif

    png_set_packing(png_ptr);

    // PNGs we're working with should have a color_type RGB
    if (color_type == PNG_COLOR_TYPE_PALETTE)
        png_set_palette_to_rgb(png_ptr);

    // Required since we selected the RGB palette
    png_read_update_info(png_ptr, info_ptr);

    // Allocate memory to _hold_ the image data now (lines 547-)
    png_bytep row_pointers[height];

    for (row = 0; row < height; ++row)
        row_pointers[row] = NULL;

    for (row = 0; row < height; ++row)
        row_pointers[row] = png_malloc(png_ptr,\
            png_get_rowbytes(png_ptr, info_ptr)
        );

    png_read_image(png_ptr, row_pointers);
    png_read_end(png_ptr, info_ptr);

    // Now clean up - the image data is in memory
    png_destroy_read_struct(&png_ptr, &info_ptr, NULL); 
    fclose(fp);

    rTuple t = { width, height, bit_depth, *row_pointers };

    return t;
}

int 
main(int argc, char *argv[])
{
    if (argc != 2) {
        printf("** Provide filename\n");
        abort();
    }

    char *fileName = argv[1];

    // get data read
    rTuple data = read_png_file(fileName);

    /* let's try computing the absolute average vector */
    uint32_t i, j, k;
    double *avV = malloc(CHANNELS * sizeof(double));
    memset(avV, 0, sizeof(double) * CHANNELS);

    double new_px[CHANNELS];
    png_bytep row, px;
    for (i = 0; i < data.height; ++i)
    {
        row = &data.datap[i];
        for (j = 0; j < data.width; ++j) 
        {
            px = &(row[j * sizeof(int)]);

            for (k = 0; k < CHANNELS; ++k) {
                new_px[k] = (double)px[k];
                avV[k] += new_px[k];
            }   
        }
    }

    double size = (double)data.width * (double)data.height;

    for (k = 0; k < CHANNELS; ++k) {
        avV[k] /= size;
        printf("channel %d: %lf\n", k + 1, avV[k]);
    }

    printf("\n");

    return 0;
}

Now with Python I'm just opening an image with a simple context manager and computing np.mean(image_data, axis=(0, 1)), which yields my result above.

Answer 1

Basically, you had a couple bugs (libpng side and pointer arithmetic) that I try to find them by comparing your code with this Github gist. The followings are the list of changes that I have made to produce the same image mean as Python NumPy.

1. In rTuple struct, you need to change the png_bytep datap to a pointer of type png_byte using: png_bytep *datap;.

2. In read_png_file, use png_set_filler to add the filler byte after reading the image. See here for more information about it.

   if(color_type == PNG_COLOR_TYPE_RGB  ||
      color_type == PNG_COLOR_TYPE_GRAY ||
      color_type == PNG_COLOR_TYPE_PALETTE)
   png_set_filler(png_ptr, 0xFF, PNG_FILLER_AFTER);
   

3. In read_png_file, update the changes before allocating the row_pointers using png_read_update_info(png_ptr, info_ptr);

4. Again, in read_png_file, change the way you are mallocing memory for image pixels using:

   png_bytep *row_pointers = (png_bytep*)malloc(sizeof(png_bytep) * height);
   for(row = 0; row < height; row++)
   {
       row_pointers[row] = malloc(png_get_rowbytes(png_ptr,info_ptr));
   }
   

5. In main, change row = &data.datap[i]; to row = data.datap[i]; as your accessing a pointer here.

I did not want to populate the answer with a code that is barely the same as the question, so if you want to just copy and paste the answer, this is the link to the complete code.