Fast loading of png files

Question

I need to process through 600 images that are stored as png-files. Each file has 8.3-ish Mb and 2400x2400 pixel (Real32). For a single processor it takes 0.93sec per image to be loaded with Import (Dell Latitude E6430, i7, windows 7, 64Bit). The processing takes 0.617sec. The processing is subject to optimization and should become even faster. So the bottle neck is the image loading. I want to increase speed as in future we will need to process 2000 images per run. Using a dual core maschine and ParallelTable I can increase speed by a factor of 1.23. Then going through the images takes ~15min. Also directly loading the png data with Import["... .png",{"Data"}] is not faster.

Is there a way to speed up reading/loading of png images? Can the loading of images be better distributed on the cores?

Answer 1

This might be faster (you will need to install libpng):

Needs["CCompilerDriver`"]

loadpng$source="
// Load a PNG image from a file into memory.
// Placed into the public domain by Mark Adler, January 24, 2016

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <setjmp.h>
#include <png.h>            // libpng

// Image returned by image_load().  For samples, Y means grey scale, RGB is
// red-green-blue, and A is an alpha channel.  The samples in each pixel are in
// the listed order in the row (e.g. RGBRGBRGB...).  Each of rows[0..height-1]
// and the rows vector itself are allocated.  image_free() should be used to
// release the allocated memory.  For one-byte samples, the samples are
// accessed using rows.b.  For two-byte samples, the samples are accessed using
// rows.w (the samples will have the correct endianess).  The original image
// may have been palettized, had less than 8-bits per pixel, and/or had a tRNS
// transparent pixel value, all of which are expanded for the returned image.
// The pre-expansion image bit depth and color type are saved in depth and
// type.
typedef struct {
    uint32_t width;             // image width in pixels (not zero)
    uint32_t height;            // image height in pixels (not zero)
    uint16_t depth;             // original (pre-expansion) bit depth
    uint16_t type;              // original (pre-expansion) color type
    uint16_t samples;           // 1 = Y, 2 = YA, 3 = RGB, 4 = RGBA
    uint16_t bytes;             // bytes per sample = 1 or 2
    size_t row_size;            // bytes per row = width * samples * bytes
    union {
        unsigned char **b;      // array of height rows, each row_size long
        uint16_t **w;           // same array for when bytes == 2
    } rows;
} image_t;

// Avoid writing a warning message to stderr -- passed to create function.
static void no_warn(png_structp png, const char *msg) {
    (void)png;
    (void)msg;
}

// Avoid writing an error message to stderr -- passed to create function.
static void no_errmsg(png_structp png, const char *msg) {
    (void)msg;
    longjmp(png_jmpbuf(png), 1);        // 1 = libpng error
}

// Free the allocations used by image.
static void image_free(image_t *image) {
    if (image->rows.b == NULL)
        return;
    for (size_t n = image->height; n;)
        free(image->rows.b[--n]);
    free(image->rows.b);
    image->rows.b = NULL;
}

// Load the PNG image from path, returning the result in image.  The return
// value is 0 on success, -1 if there is a memory allocation error, 1 if the
// input is invalid, or 2 if path could not be opened.
static int image_load(const char *path, image_t *image) {
    // create libpng read structure
    png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING,
                                             NULL, no_errmsg, no_warn);
    if (png == NULL)
        return -1;

    // set up error handling
    FILE * volatile fp = NULL;
    image->rows.b = NULL;
    png_infop info = NULL;
    {
        int jmp = setjmp(png_jmpbuf(png));
        if (jmp) {                      // someone longjmp'ed here
            image_free(image);
            png_destroy_read_struct(&png, &info, NULL);
            if (fp != NULL)
                fclose(fp);
            return jmp;
        }
    }

    // create libpng image information structure
    info = png_create_info_struct(png);
    if (info == NULL)
        longjmp(png_jmpbuf(png), -1);

    // open input file
    fp = fopen(path, \"rb\");
    if (fp == NULL)
        longjmp(png_jmpbuf(png), 2);
    png_init_io(png, fp);

    // read png header with the image information, save header data
    png_read_info(png, info);
    image->width = png_get_image_width(png, info);
    image->height = png_get_image_height(png, info);
    image->depth = png_get_bit_depth(png, info);
    image->type = png_get_color_type(png, info);

    // expand palettes, expand samples to 8 or 16 bits, and expand tRNS to an
    // alpha channel if present -- save the expanded image information for how
    // the image will be returned
    png_set_expand(png);
    png_read_update_info(png, info);
    {
        unsigned type = png_get_color_type(png, info);
        image->samples = (type >> 2) + (type & 2) + 1;
    }
    image->bytes = png_get_bit_depth(png, info) >> 3;
    image->row_size = png_get_rowbytes(png, info);

    // allocate space for the png image (use calloc() to initialize as NULLs)
    image->rows.b = calloc(image->height, sizeof(unsigned char *));
    if (image->rows.b == NULL)
        longjmp(png_jmpbuf(png), -1);
    for (size_t n = 0; n < image->height; n++)
        if ((image->rows.b[n] = malloc(image->row_size)) == NULL)
            longjmp(png_jmpbuf(png), -1);

    // match the endianess of 16-bit samples to this machine
    if (image->bytes == 2) {
        const unsigned little = 1;
        if (*(const unsigned char *)&little)
            png_set_swap(png);
    }

    // read in the png image and return success (if no longjmp's)
    png_read_image(png, image->rows.b);
    png_read_end(png, info);
    png_destroy_read_struct(&png, &info, NULL);
    fclose(fp);
    return 0;
}

#include \"mathlink.h\"
#include \"WolframLibrary.h\"

DLLEXPORT mint WolframLibrary_getVersion( ) {
   return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
   return LIBRARY_NO_ERROR;
}

DLLEXPORT void WolframLibrary_uninitialize( ) {
   return;
}

static void report(WolframLibraryData libData, const char *name,
                   const char *err, const char *text) {
   MLINK link = libData->getMathLink(libData);
      /* name::err = \"`1`\" */
   MLPutFunction(link, \"EvaluatePacket\", 1);
   MLPutFunction(link, \"Set\", 2);
   MLPutFunction(link, \"MessageName\", 2);
   MLPutSymbol(link, name);
   MLPutString(link, err);
   MLPutString(link, \"`1`\");
   libData->processMathLink(link);
   if (MLNextPacket(link) == RETURNPKT)
      MLNewPacket(link);
      /* Message[name::err, text] */
   MLPutFunction(link, \"EvaluatePacket\", 1);
   MLPutFunction(link, \"Message\", 2);
   MLPutFunction(link, \"MessageName\", 2);
   MLPutSymbol(link, name);
   MLPutString(link, err);
   MLPutString(link, text);
   libData->processMathLink(link);
   if (MLNextPacket(link) == RETURNPKT)
       MLNewPacket(link);
}

DLLEXPORT mint loadpng(WolframLibraryData libData, MLINK link) {
   const char *path;
   int len;
   const char *errmsg[] = {
      \"Out of memory\",                           /* -1 */
      \"Unexpected return code\",                  /* < -1 or > 2 */
      \"Invalid PNG input\",                       /* 1 */
      \"Could not open file\"                      /* 2 */
   };

   if (MLTestHead(link, \"List\", &len) && len == 1 &&
       MLGetType(link) == MLTKSTR && MLGetString(link, &path)) {
      MLNewPacket(link);
      image_t image;
      int ret = image_load(path, &image);
      MLReleaseString(link, path);
      int *a;
      if (ret == 0 && (ret = -1,
           (a = malloc(((sizeof(int) * image.height) * image.width) *
                      image.samples)) != NULL)) {
         size_t n = 0;
         for (size_t i = 0; i < image.height; i++) {
             size_t m = 0;
             for (size_t j = 0; j < image.width; j++)
                 for (size_t k = 0; k < image.samples; k++)
                     a[n++] = image.bytes == 1 ? image.rows.b[i][m++] :
                                                 image.rows.w[i][m++];
         }
         image_free(&image);
         const int dims[3] = {image.height, image.width, image.samples};
         MLPutFunction(link, \"Image\", 3);
         MLPutInteger32Array(link, a, dims, NULL, 3);
         free(a);
         MLPutString(link,
                     image.type == 0 && image.depth == 1 ? \"Bit\" :
                     image.bytes == 1 ? \"Byte\" : \"Bit16\");
         MLPutFunction(link, \"Rule\", 2);
         MLPutSymbol(link, \"ColorSpace\");
         MLPutString(link,
                     image.type == 0 && image.depth == 1 ? \"Automatic\" :
                     image.samples < 3 ? \"Grayscale\" : \"RGB\");
      }
      else {
         MLPutSymbol(link, \"Null\");
         report(libData, \"loadpng\",
                ret < 0 ? \"memerr\" : \"invalid\",
                errmsg[ret < -1 || ret > 2 ? 1 : ret + 1]);
      }
   }
   else {
      MLClearError(link);
      MLNewPacket(link);
      MLPutSymbol(link, \"Null\");
      report(libData, \"loadpng\", \"arg\",
             \"String expected at position 1 as the only argument\");
   }
   MLClearError(link);
   return LIBRARY_NO_ERROR;
}";

If[Head[lib]===String,LibraryUnload[lib]]
lib=CreateLibrary[loadpng$source,"loadpng","CompileOptions"->"-lpng",
  "ShellOutputFunction"->Print]
loadpng=LibraryFunctionLoad[lib,"loadpng",LinkObject,LinkObject]

Then you can:

loadpng["myimage.png"]

which will return the same image as:

Import["myimage.png"]