channel separation, color decorrelation, and prediction filtering as preprocessing for lossless RGB image compression

image_predict.c

#if 0  // self-compiling code
gcc -std=c99 -Wall -Wextra -pedantic -g -O3 -march=native $0 -o image_predict || exit 1
exec ./image_predict $*
#endif

// Example code for pre-filtering operations that transform RGB images into
// representations that are better suited for lossless compression
// (color channel separation and decorrelation, pixel prediction).

// Takes a .ppm file as an input and generates multiple files in the form
//   image_predict_<sep>_<pred>_data.pgm
// with various separator/predictor choices.

// This program does *not* perform the actual compression; this has to be
// done externally. It also does not decode the _data.pgm files, as it's a
// pure proof-of-concept implementation, not a proper codec; it verifies
// that the transformation and back-transformation is lossless though.

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>
#include <unistd.h>

////////////////////////////////////////////////////////////////////////////////

typedef struct _pnm {
    int width, height, ncomp;
    uint8_t *data;
    uint8_t rawdata[1];
} pnm_image_t;

pnm_image_t*pnm_load(const char*filename){FILE*f;pnm_image_t*i=NULL;size_t s;
char*p,t;size_t n,h[4]={0,0,0,0};if(!filename)return NULL;f=fopen(filename,"rb"
);if(!f)return NULL;fseek(f,0,SEEK_END);s=(size_t)ftell(f);if(!s||(s>99999999))
goto _;i=(pnm_image_t*)malloc(sizeof(pnm_image_t)+s-1);if(!i)goto _;fseek(f,0,
SEEK_SET);p=(char*)i->rawdata;if((fread((void*)p,1,s,f)!=s)||(*p++!='P'))goto _
;for(n=0,t=0;(n<4)&&--s;p++){if(t==2)t=(*p=='\n')?0:2;else if isdigit(*p){h[n]=
(10*h[n])+*p-'0';t=1;}else if(isspace(*p)||(*p=='#')){if(t)n++;t=(*p=='#')?2:0;
}else break;}i->width=(int)h[1];i->height=(int)h[2];i->ncomp=(h[0]==6)?3:1;i->
data=(unsigned char*)p;if((n==4)&&((h[0]==5)||(h[0]==6))&&h[1]&&(h[1]<65536)&&h
[2]&&(h[2]<65536)&&h[3]&&(h[3]<256)&&(s>=(h[1]*h[2]*i->ncomp))){fclose(f);
return i;}_:fclose(f);free((void*)i);return NULL;}

////////////////////////////////////////////////////////////////////////////////

typedef void separate_func(const uint8_t* src, uint8_t* p0, uint8_t* p1, uint8_t* p2, int npix);
typedef void interleave_func(const uint8_t* p0, const uint8_t* p1, const uint8_t* p2, uint8_t* dest, int npix);

void separate_rgb(const uint8_t* src, uint8_t* r, uint8_t* g, uint8_t* b, int npix) {
    while (npix--) {
        *r++ = *src++;
        *g++ = *src++;
        *b++ = *src++;
    }
}

void interleave_rgb(const uint8_t* r, const uint8_t* g, const uint8_t* b, uint8_t* dest, int npix) {
    while (npix--) {
        *dest++ = *r++;
        *dest++ = *g++;
        *dest++ = *b++;
    }
}

void separate_bdgdr(const uint8_t* src, uint8_t* p_b, uint8_t* p_dg, uint8_t* p_dr, int npix) {
    while (npix--) {
        uint8_t r = *src++;
        uint8_t g = *src++;
        uint8_t b = *src++;
        *p_b++ = b;
        *p_dg++ = g - b + 128;
        *p_dr++ = r - b + 128;
    }
}

void interleave_bdgdr(const uint8_t* p_b, const uint8_t* p_dg, const uint8_t* p_dr, uint8_t* dest, int npix) {
    while (npix--) {
        uint8_t b = *p_b++;
        uint8_t dg = *p_dg++;
        uint8_t dr = *p_dr++;
        *dest++ = b + dr - 128;
        *dest++ = b + dg - 128;
        *dest++ = b;
    }
}

void separate_gdrdb(const uint8_t* src, uint8_t* p_g, uint8_t* p_dr, uint8_t* p_db, int npix) {
    while (npix--) {
        uint8_t r = *src++;
        uint8_t g = *src++;
        uint8_t b = *src++;
        *p_g++ = g;
        *p_dr++ = r - g + 128;
        *p_db++ = b - g + 128;
    }
}

void interleave_gdrdb(const uint8_t* p_g, const uint8_t* p_dr, const uint8_t* p_db, uint8_t* dest, int npix) {
    while (npix--) {
        uint8_t g = *p_g++;
        uint8_t dr = *p_dr++;
        uint8_t db = *p_db++;
        *dest++ = g + dr - 128;
        *dest++ = g;
        *dest++ = g + db - 128;
    }
}

void separate_gab(const uint8_t* src, uint8_t* p_g, uint8_t* p_a, uint8_t* p_b, int npix) {
    while (npix--) {
        uint8_t r = *src++;
        uint8_t g = *src++;
        uint8_t b = *src++;
        *p_g++ = g;
        *p_a++ = r - g + 128;
        *p_b++ = b - ((r + g) >> 1) + 128;
    }
}

void interleave_gab(const uint8_t* p_g, const uint8_t* p_a, const uint8_t* p_b, uint8_t* dest, int npix) {
    while (npix--) {
        uint8_t g = *p_g++;
        uint8_t a = *p_a++;
        uint8_t b = *p_b++;
        uint8_t r = g + a - 128;
        *dest++ = r;
        *dest++ = g;
        *dest++ = ((r + g) >> 1) + b - 128;
    }
}

static const struct separator {
    const char *name;
    separate_func* sep;
    interleave_func* unsep;
    const char* desc;
} separators[] = {
    { "rgb",   separate_rgb,   interleave_rgb,   "RGB decomposition, no decorrelation" },
    { "bdgdr", separate_bdgdr, interleave_bdgdr, "B (G-B) (R-B) decomposition" },
    { "gdrdb", separate_gdrdb, interleave_gdrdb, "G (R-G) (B-G) decomposition" },
    { "gab",   separate_gab,   interleave_gab,   "G (R-G) (B-avg(R,G)) decomposition" },
    { NULL, }
};

////////////////////////////////////////////////////////////////////////////////

typedef void predict_func(uint8_t* buffer, int width, int height);

void null_predictor(uint8_t* buffer, int width, int height) {
    (void)buffer, (void)width, (void)height;
}

void predict_left(uint8_t* p, int width, int height) {
    p += width * (height - 1) + width - 1;
    for (int y = height - 1;  y >= 0;  --y, --p) {
        for (int x = width - 1;  x;  --x, --p) {
            p[0] -= p[-1] - 128;  // center: predict from left pixel
        }
        if (y) {
            p[0] -= p[-width] - 128;  // left edge: predict from top pixel
        }  // at x=0, y=0, do nothing
    }
}

void recon_left(uint8_t* p, int width, int height) {
    for (int y = 0;  y < height;  ++y) {
        if (y) {
            p[0] += p[-width] - 128;  // left edge: predict from top pixel
        }
        ++p;
        for (int x = 1;  x < width;  ++x, ++p) {
            p[0] += p[-1] - 128;  // center: predict from left pixel
        }
    }
}

void predict_grad(uint8_t* p, int width, int height) {
    p += width * (height - 1) + width - 1;
    for (int y = height - 1;  y;  --y, --p) {
        for (int x = width - 1;  x;  --x, --p) {
            p[0] -= p[-1] + p[-width] - p[-width-1] - 128;  // center: predict N + E - NE
        }
        p[0] -= p[-width] - 128;  // left edge: predict from top pixel
    }
    for (int x = width - 1;  x;  --x, --p) {
        p[0] -= p[-1] - 128;  // top edge: predict from left pixel
    }
}

void recon_grad(uint8_t* p, int width, int height) {
    ++p;
    for (int x = 1;  x < width;  ++x, ++p) {
        p[0] += p[-1] - 128;  // top edge: predict from left pixel
    }
    for (int y = 1;  y < height;  ++y) {
        if (y) {
            p[0] += p[-width] - 128;  // left edge: predict from top pixel
        }
        ++p;
        for (int x = 1;  x < width;  ++x, ++p) {
            p[0] += p[-1] + p[-width] - p[-width-1] - 128;  // center: predict N + E - NE
        }
    }
}

inline uint8_t clamp_pred(int16_t a, int16_t b, int16_t c) {
    int16_t max = (a > b) ? ((a > c) ? a : c) : ((b > c) ? b : c);
    int16_t min = (a < b) ? ((a < c) ? a : c) : ((b < c) ? b : c);
    int16_t pred = a + b - c;
    return (uint8_t)((pred < min) ? min : (pred > max) ? max : pred);
}

void predict_clamp(uint8_t* p, int width, int height) {
    p += width * (height - 1) + width - 1;
    for (int y = height - 1;  y;  --y, --p) {
        for (int x = width - 1;  x;  --x, --p) {
            p[0] -= clamp_pred(p[-1], p[-width], p[-width-1]) - 128;  // center
        }
        p[0] -= p[-width] - 128;  // left edge: predict from top pixel
    }
    for (int x = width - 1;  x;  --x, --p) {
        p[0] -= p[-1] - 128;  // top edge: predict from left pixel
    }
}

void recon_clamp(uint8_t* p, int width, int height) {
    ++p;
    for (int x = 1;  x < width;  ++x, ++p) {
        p[0] += p[-1] - 128;  // top edge: predict from left pixel
    }
    for (int y = 1;  y < height;  ++y) {
        if (y) {
            p[0] += p[-width] - 128;  // left edge: predict from top pixel
        }
        ++p;
        for (int x = 1;  x < width;  ++x, ++p) {
            p[0] += clamp_pred(p[-1], p[-width], p[-width-1]) - 128;  // center
        }
    }
}

static const struct predictor {
    const char *name;
    predict_func* pred;
    predict_func* unpred;
    const char *desc;
} predictors[] = {
    { "raw",   null_predictor, null_predictor, "no prediction" },
    { "left",  predict_left,   recon_left,     "predict from left pixel" },
    { "grad",  predict_grad,   recon_grad,     "predict gradient (L+T-LT)" },
    { "clamp", predict_clamp,  recon_clamp,    "predict gradient with clamping" },
    { NULL, }
};

////////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[]) {
    if (argc != 2) {
        printf("Usage: %s <input.ppm>\n", argv[0]);
        return 2;
    }

    pnm_image_t* img = pnm_load(argv[1]);
    if (!img || (img->ncomp != 3)) {
        fprintf(stderr, "Error: input must be a valid 8-bit RGB PPM file\n");
        return 1;
    }

    void *outbuf = malloc(img->width * img->height * 3);
    void *checkbuf = malloc(img->width * img->height * 3);
    assert(outbuf != NULL);
    uint8_t *planes[3];
    planes[0] = (uint8_t*) outbuf;
    planes[1] = planes[0] + (img->width * img->height);
    planes[2] = planes[1] + (img->width * img->height);

    printf("Separators:\n");
    for (const struct separator* sep = separators;  sep->name;  ++sep) {
        printf("  - %-10s -> %s\n", sep->name, sep->desc);
    }
    printf("Prediction filters:\n");
    for (const struct predictor* pred = predictors;  pred->name;  ++pred) {
        printf("  - %-10s -> %s\n", pred->name, pred->desc);
    }

    for (const struct separator* sep = separators;  sep->name;  ++sep) {
        for (const struct predictor* pred = predictors;  pred->name;  ++pred) {
            char variant[32], filename[64];
            sprintf(variant, "%s_%s", sep->name, pred->name);
            
            sep->sep(img->data, planes[0], planes[1], planes[2], img->width * img->height);
            pred->pred(planes[0], img->width, img->height);
            pred->pred(planes[1], img->width, img->height);
            pred->pred(planes[2], img->width, img->height);

            sprintf(filename, "image_predict_%s_data.pgm", variant);
            FILE *f = fopen(filename, "wb");
            if (f) {
                fprintf(f, "P5\n%d %d\n255\n", img->width, img->height * 3);
                fwrite(outbuf, img->width * img->height * 3, 1, f);
                fclose(f);
            }

            pred->unpred(planes[0], img->width, img->height);
            pred->unpred(planes[1], img->width, img->height);
            pred->unpred(planes[2], img->width, img->height);
            sep->unsep(planes[0], planes[1], planes[2], (uint8_t*)checkbuf, img->width * img->height);

            sprintf(filename, "image_predict_%s_recon.ppm", variant);
            if (!memcmp((const void*)img->data, checkbuf, img->width * img->height * 3)) {
                printf("%s: OK\n", variant);
                unlink(filename);
            } else {
                fprintf(stderr, "WARNING: %s reconstruction doesn't match original!\n", variant);
                f = fopen(filename, "wb");
                if (f) {
                    fprintf(f, "P6\n%d %d\n255\n", img->width, img->height);
                    fwrite(checkbuf, img->width * img->height * 3, 1, f);
                    fclose(f);
                }
            } // comparison
        } // predictor loop
    } // separator loop

    free(checkbuf);
    free(outbuf);
    free(img);
    return 0;
}