imread.c

   1 /* PHOROS -- Photogrammetric Road Survey */
   2 /* Copyright (C) 2016 Bert Burgemeister */
   3
   4 /* This program is free software; you can redistribute it and/or modify */
   5 /* it under the terms of the GNU General Public License as published by */
   6 /* the Free Software Foundation; either version 2 of the License, or */
   7 /* (at your option) any later version. */
   8
   9 /* This program is distributed in the hope that it will be useful, */
  10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
  11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
  12 /* GNU General Public License for more details. */
  13
  14 /* You should have received a copy of the GNU General Public License along */
  15 /* with this program; if not, write to the Free Software Foundation, Inc., */
  16 /* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
  17
  18 #include <stdlib.h>
  19 #include <stdio.h>
  20 #include <stdint.h>
  21 #include <malloc.h>
  22 #include <math.h>
  23 #include <png.h>
  24
  25 void
  26 cplt_horiz(unsigned char *rawp, png_bytep row, int y, int x, int color, int width)
  27 {
  28         row[3 * x + color] = (rawp[-1] + rawp[1]) / 2;
  29 }
  30
  31 void
  32 cplt_vert(unsigned char *rawp, png_bytep row, int y, int x, int color, int width)
  33 {
  34         row[3 * x + color] = (rawp[-width] + rawp[width]) / 2;
  35 }
  36
  37 void
  38 cplt_squ(unsigned char *rawp, png_bytep row, int y, int x, int color, int width)
  39 {
  40         row[3 * x + color] = (rawp[-1] + rawp[1] + rawp[-width] + rawp[width]) / 4;
  41 }
  42
  43 void
  44 cplt_diag(unsigned char *rawp, png_bytep row, int y, int x, int color, int width)
  45 {
  46         row[3 * x  + color] = (rawp[-width - 1] + rawp[-width + 1] +
  47                                rawp[width - 1] + rawp[width + 1]) / 4;
  48 }
  49
  50 void
  51 cplt_g_on_r(unsigned char *rawp, png_bytep row, int y, int x, int width)
  52 {
  53         cplt_horiz(rawp, row, y, x, 0, width);
  54         cplt_vert(rawp, row, y, x, 2, width);
  55 }
  56
  57 void
  58 cplt_g_on_b(unsigned char *rawp, png_bytep row, int y, int x, int width)
  59 {
  60         cplt_horiz(rawp, row, y, x, 2, width);
  61         cplt_vert(rawp, row, y, x, 0, width);
  62 }
  63
  64 void
  65 cplt_r(unsigned char *rawp, png_bytep row, int y, int x, int width)
  66 {
  67         cplt_squ(rawp, row, y, x, 1, width);
  68         cplt_diag(rawp, row, y, x, 2, width);
  69 }
  70
  71 void
  72 cplt_b(unsigned char *rawp, png_bytep row, int y, int x, int width)
  73 {
  74         cplt_squ(rawp, row, y, x, 1, width);
  75         cplt_diag(rawp, row, y, x, 0, width);
  76 }
  77
  78 void
  79 colrz_r(unsigned char *rawp, png_bytep row, int y, int x)
  80 {
  81         row[3 * x + 0] = rawp[0];
  82 }
  83
  84 void
  85 colrz_g(unsigned char *rawp, png_bytep row, int y, int x)
  86 {
  87         row[3 * x + 1] = rawp[0];
  88 }
  89
  90 void
  91 colrz_b(unsigned char *rawp, png_bytep row, int y, int x)
  92 {
  93         row[3 * x + 2] = rawp[0];
  94 }
  95
  96 void
  97 raise_color(png_bytep row, int x, double *colr_raisr)
  98 {
  99         int i, val;
 100
 101         for (i = 0; i < 3; i++) {
 102                 val = round(row[3 * x + i] * colr_raisr[i]);
 103                 if (val > 255)
 104                         val = 255;
 105                 row[3 * x + i] = val;
 106         }
 107 }
 108
 109 void
 110 raise_noop(png_bytep row, int x, double *colr_raisr)
 111 {
 112 }
 113
 114 /* structure to store in-memory PNG in */
 115 struct mem_encode
 116 {
 117         char *buffer;
 118         size_t size;
 119 };
 120
 121 void
 122 write_png_data(png_structp png_ptr, png_bytep data, png_size_t length)
 123 {
 124         struct mem_encode *p = (struct mem_encode *)png_get_io_ptr(png_ptr);
 125         size_t nsize = p->size + length;
 126         /* allocate or grow buffer */
 127         if (p->buffer)
 128                 p->buffer = realloc(p->buffer, nsize);
 129         else
 130                 p->buffer = malloc(nsize);
 131         if (!p->buffer)
 132                 png_error(png_ptr, "Write Error");
 133         /* copy new bytes to end of buffer */
 134         memcpy(p->buffer + p->size, data, length);
 135         p->size += length;
 136 }
 137
 138 void
 139 flush_png(png_structp png_ptr)
 140 {
 141 }
 142
 143 int
 144 uncompressed2png(struct mem_encode *mem_png, int width, int height, int channels,
 145                  int *bayerpat, double *color_raiser, unsigned char *uncompressed)
 146 {
 147         int retval = 0;
 148         png_structp png_ptr = NULL;
 149         png_infop info_ptr = NULL;
 150         png_bytep row = NULL;
 151         void (*colrz_ev_ev)(unsigned char *, png_bytep, int, int) = NULL;
 152         void (*colrz_ev_od)(unsigned char *, png_bytep, int, int) = NULL;
 153         void (*colrz_od_ev)(unsigned char *, png_bytep, int, int) = NULL;
 154         void (*colrz_od_od)(unsigned char *, png_bytep, int, int) = NULL;
 155         void (*cplt_ev_ev)(unsigned char *, png_bytep, int, int, int) = NULL;
 156         void (*cplt_ev_od)(unsigned char *, png_bytep, int, int, int) = NULL;
 157         void (*cplt_od_ev)(unsigned char *, png_bytep, int, int, int) = NULL;
 158         void (*cplt_od_od)(unsigned char *, png_bytep, int, int, int) = NULL;
 159         void (*raise)(unsigned char *, int, double *) = raise_noop;
 160         int x, y, x_ev, x_od, y_ev, y_od;
 161         int i;
 162
 163         for (i = 0; i < 3; i++)
 164                 if (fabs(color_raiser[i] - 1) > 0.00001)
 165                         raise = raise_color;
 166         if (channels == 3) {
 167                 if (bayerpat[0] == 0x0000ff) { /* red */
 168                         colrz_ev_ev = colrz_r;
 169                         colrz_ev_od = colrz_g;
 170                         colrz_od_ev = colrz_g;
 171                         colrz_od_od = colrz_b;
 172                         cplt_ev_ev = cplt_r;
 173                         cplt_ev_od = cplt_g_on_r;
 174                         cplt_od_ev = cplt_g_on_b;
 175                         cplt_od_od = cplt_b;
 176                 } else if (bayerpat[0] == 0xff0000) { /* blue */
 177                         colrz_ev_ev = colrz_b;
 178                         colrz_ev_od = colrz_g;
 179                         colrz_od_ev = colrz_g;
 180                         colrz_od_od = colrz_r;
 181                         cplt_ev_ev = cplt_b;
 182                         cplt_ev_od = cplt_g_on_b;
 183                         cplt_od_ev = cplt_g_on_r;
 184                         cplt_od_od = cplt_b;
 185                 } else if (bayerpat[0] == 0x00ff00) { /* green */
 186                         if (bayerpat[1] == 0x0000ff) { /* red */
 187                                 colrz_ev_ev = colrz_g;
 188                                 colrz_ev_od = colrz_r;
 189                                 colrz_od_ev = colrz_b;
 190                                 colrz_od_od = colrz_g;
 191                                 cplt_ev_ev = cplt_g_on_r;
 192                                 cplt_ev_od = cplt_r;
 193                                 cplt_od_ev = cplt_b;
 194                                 cplt_od_od = cplt_g_on_b;
 195                         } else if (bayerpat[1] == 0xff0000) { /* blue */
 196                                 colrz_ev_ev = colrz_g;
 197                                 colrz_ev_od = colrz_b;
 198                                 colrz_od_ev = colrz_r;
 199                                 colrz_od_od = colrz_g;
 200                                 cplt_ev_ev = cplt_g_on_b;
 201                                 cplt_ev_od = cplt_b;
 202                                 cplt_od_ev = cplt_r;
 203                                 cplt_od_od = cplt_g_on_r;
 204                         } else {
 205                                 retval = 2; /* all green is not a Bayer pattern */
 206                                 goto finalize;
 207                         }
 208                 } else {
 209                         retval = 3;     /* first byte neither 0x00 nor 0xff */
 210                         goto finalize;
 211                 }
 212         }
 213         png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
 214         if (png_ptr == NULL) {
 215                 retval = 11; /* Could not allocate write struct */
 216                 goto finalize;
 217         }
 218         info_ptr = png_create_info_struct(png_ptr);
 219         if (info_ptr == NULL) {
 220                 retval = 12; /* Could not allocate info struct */
 221                 goto finalize;
 222         }
 223         if (setjmp(png_jmpbuf(png_ptr))) {
 224                 retval = 13; /* Error during png creation */
 225                 goto finalize;
 226         }
 227         mem_png->buffer = NULL;
 228         mem_png->size = 0;
 229         png_set_write_fn(png_ptr, mem_png, write_png_data, flush_png);
 230         /* Write header (8 bit colour depth) */
 231         row = malloc(channels * width * sizeof(png_byte));
 232         if (row == NULL) {
 233                 retval = 21;
 234                 goto finalize;
 235         }
 236         /* Write image data */
 237         if (channels == 3) {
 238                 png_set_IHDR(png_ptr, info_ptr, width, height,
 239                              8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
 240                              PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
 241                 png_write_info(png_ptr, info_ptr);
 242                 for (y_ev = 0, y_od = 1;
 243                      y_ev < height && y_od < height;
 244                      y_ev += 2, y_od += 2) {
 245                         for (x_ev = 0, x_od = 1;
 246                              x_ev < width && x_od < width;
 247                              x_ev += 2, x_od += 2) {
 248                                 colrz_ev_ev(uncompressed + y_ev * width + x_ev,
 249                                             row, y_ev, x_ev);
 250                                 colrz_ev_od(uncompressed + y_ev * width + x_od,
 251                                             row, y_ev, x_od);
 252                                 if (y_ev > 0 && y_ev < height - 1) {
 253                                         cplt_ev_ev(uncompressed + y_ev * width + x_ev,
 254                                                    row, y_ev, x_ev, width);
 255                                         cplt_ev_od(uncompressed + y_ev * width + x_od,
 256                                                    row, y_ev, x_od, width);
 257                                 }
 258                                 raise(row, x_ev, color_raiser);
 259                                 raise(row, x_od, color_raiser);
 260                         }
 261                         png_write_row(png_ptr, row);
 262                         for (x_ev = 0, x_od = 1;
 263                              x_ev < width && x_od < width;
 264                              x_ev += 2, x_od += 2) {
 265                                 colrz_od_ev(uncompressed + y_od * width + x_ev,
 266                                             row, y_od, x_ev);
 267                                 colrz_od_od(uncompressed + y_od * width + x_od,
 268                                             row, y_od, x_od);
 269                                 if (y_od < height - 1) {
 270                                         cplt_od_ev(uncompressed + y_od * width + x_ev,
 271                                                    row, y_od, x_ev, width);
 272                                         cplt_od_od(uncompressed + y_od * width + x_od,
 273                                                    row, y_od, x_od, width);
 274                                 }
 275                                 raise(row, x_ev, color_raiser);
 276                                 raise(row, x_od, color_raiser);
 277                         }
 278                         png_write_row(png_ptr, row);
 279                 }
 280         } else if (channels == 1) {
 281                 png_set_IHDR(png_ptr, info_ptr, width, height,
 282                              8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE,
 283                              PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
 284                 png_write_info(png_ptr, info_ptr);
 285                 for (y = 0; y < height; y++) {
 286                         for (x = 0; x < width; x++) {
 287                                 row[x] = uncompressed[y * width + x];
 288                         }
 289                         png_write_row(png_ptr, row);
 290                 }
 291         } else {
 292                 retval = 6;       /* wrong number of channels */
 293                 goto finalize;
 294         }
 295         png_write_end(png_ptr, NULL);
 296 finalize:
 297         if (info_ptr != NULL) png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
 298         if (png_ptr != NULL) png_destroy_write_struct(&png_ptr, NULL);
 299         if (row != NULL) free(row);
 300         return retval;
 301 }
 302
 303 int
 304 huffdecode(int width, int height, unsigned char *uncompressed,
 305            unsigned char hcode[static 4 * 511], unsigned char hlen[static 511],
 306            unsigned char *compressed)
 307 {
 308         int i, j, m;
 309         int maxlen = 0;
 310         int code;
 311         struct huffval {
 312                 int val;
 313                 int len;
 314         } *htable;
 315         int row, column, cidx;
 316         uint64_t hc_mask;
 317         unsigned int rowhead;
 318
 319         for (i = 0; i < 511; i++)
 320                 if (hlen[i] > maxlen)
 321                         maxlen = hlen[i];
 322         if ((htable = malloc((1 << maxlen) * sizeof(struct huffval))) == NULL)
 323                 return 31;
 324         /* htable may well become too big to fit into memory.  Maybe
 325         we should act a bit smarter and handle longer huffman codes
 326         differently. */
 327         for (i = 0; i < 511; i++) {
 328                 if (hlen[i] > 0) {
 329                         code = hcode[4 * i] << 24 | hcode[4 * i + 1] << 16 |
 330                                 hcode[4 * i + 2] << 8 | hcode[4 * i + 3];
 331                         code <<= maxlen - hlen[i];
 332                         htable[code].val = i - 255;
 333                         htable[code].len = hlen[i];
 334                         for (j = 0; j < (1 << maxlen - hlen[i]); j++) {
 335                                 htable[code | j].val = i - 255;
 336                                 htable[code | j].len = hlen[i];
 337                         }
 338                 }
 339         }
 340         hc_mask = 0;
 341         for (i = 0, m = 1; i < maxlen; i++, m <<= 1)
 342                 hc_mask |= m;
 343         cidx = 0;
 344         for (row = 0; row < height; row++) {
 345                 rowhead = compressed[cidx / 8] << 16;
 346                 cidx += 8;
 347                 rowhead |= compressed[cidx / 8] << 8;
 348                 cidx += 8;
 349                 rowhead |= compressed[cidx / 8];
 350                 rowhead >>= 8 - cidx % 8;
 351                 uncompressed[width * row] = (rowhead >> 8) & 0xff;
 352                 uncompressed[width * row + 1] = rowhead & 0xff;
 353                 for (column = 2; column < width; column++) {
 354                         div_t cidx_d, maxlen_d;
 355                         int rem_bits, r;
 356                         uint64_t hc;
 357
 358                         cidx_d = div(cidx, 8);
 359                         rem_bits = maxlen - (8 - cidx_d.rem);
 360                         cidx += 8 - cidx_d.rem;
 361                         maxlen_d = div(maxlen, 8);
 362                         hc = compressed[cidx_d.quot];
 363                         for (i = rem_bits; i > 0; i -= 8, cidx += 8) {
 364                                 hc <<= 8;
 365                                 hc |= compressed[cidx / 8];
 366                         }
 367                         if ((r = rem_bits % 8) > 0) {
 368                                 hc >>= 8 - r;
 369                                 cidx -= 8 - r;
 370                         }
 371                         hc &= hc_mask;
 372                         uncompressed[width * row + column] =
 373                                 uncompressed[width * row + column - 2] - htable[hc].val;
 374                         cidx -= maxlen - htable[hc].len;
 375                 }
 376         }
 377         free(htable);
 378         return 0;
 379 }
 380
 381 void
 382 reverse(unsigned char *base, size_t size)
 383 {
 384         int i, j;
 385         char tmp;
 386
 387         for (i = 0, j = size - 1; i < j; i++, j--) {
 388                 tmp = base[i];
 389                 base[i] = base[j];
 390                 base[j] = tmp;
 391         }
 392 }
 393
 394 void
 395 brighten(unsigned char *base, size_t size)
 396 {
 397         size_t i;
 398         unsigned char min = UCHAR_MAX, max = 0, range;
 399
 400         for (i = 0; i < size; i++) {
 401                 if (base[i] < min)
 402                         min = base[i];
 403                 if (base[i] > max)
 404                         max = base[i];
 405         }
 406         range = max - min;
 407         if (max < 200)
 408                 for (i = 0; i < size; i++)
 409                         base[i] = (base[i] - min) * UCHAR_MAX / range;
 410 }
 411
 412 int
 413 png2mem(char *path, int start, int len, int width, int height,
 414         int channels, int *bayer_pattern, int compr_mode,
 415         unsigned char *uncompressed, unsigned char *compressed,
 416         struct mem_encode *mem_png,
 417         int reversep, int brightenp, double *color_raiser)
 418 {
 419         FILE *in;
 420         unsigned char hlen[511], hcode[511 * 4];
 421         int retval;
 422
 423         if ((in = fopen(path, "r")) == NULL)
 424                 return 1;
 425         fseek(in, start, SEEK_CUR);
 426         if (compr_mode == 1 || compr_mode == 2) {
 427                 fread(hcode, sizeof(hcode[0]), 511 * 4, in);
 428                 fread(hlen, sizeof(hlen[0]), 511, in);
 429                 fread(compressed, sizeof(compressed[0]), len, in);
 430                 fclose(in);
 431                 if ((retval = huffdecode(width, height, uncompressed, hcode,
 432                                          hlen, compressed)) != 0)
 433                         return retval;
 434                 if (reversep)
 435                         reverse(uncompressed, width * height);
 436                 if (brightenp)
 437                         brighten(uncompressed, width * height);
 438                 if ((retval = uncompressed2png(mem_png, width, height, channels,
 439                                                bayer_pattern, color_raiser,
 440                                                uncompressed)) != 0)
 441                         return retval;
 442                 return 0;
 443         } else if (compr_mode == 0) { /* untested */
 444                 fread(uncompressed, sizeof(uncompressed[0]), width * height, in);
 445                 fclose(in);
 446                 if (reversep)
 447                         reverse(uncompressed, width * height);
 448                 if (brightenp)
 449                         brighten(uncompressed, width * height);
 450                 if ((retval = uncompressed2png(mem_png, width, height, channels,
 451                                                bayer_pattern, color_raiser,
 452                                                uncompressed)) != 0)
 453                         return retval;
 454                 return 0;
 455         } else {
 456                 fclose(in);
 457                 return 5;
 458         }
 459 }
 460
 461 int
 462 ping(int n) {
 463         return n;
 464 }
 465
 466 int
 467 main(int argc, char *argv[])
 468 {
 469         FILE *fp;
 470         unsigned char compressed[1390353];
 471         unsigned char uncompressed[1700 * 1500];
 472         int width = 1700, height = 1500;
 473         struct mem_encode mp;
 474
 475         setvbuf(stdout, NULL, _IONBF, 0);
 476
 477         png2mem("mitsa005_CCD_Front_PULNIX_13.pictures", 2247 + 17, 1390353, width, height, 3, (int[2]){0x00ff00, 0x0000ff}, 2, uncompressed, compressed, &mp, 0, 0, (double[3]){1, 1, 1});
 478         fp = fopen("o0.png", "wb");
 479         if (fp == NULL) {
 480                 fprintf(stderr, "Could not open file %s for writing\n", "o0.png");
 481         }
 482         fwrite(mp.buffer, mp.size, sizeof(mp.buffer[0]), fp);
 483
 484         png2mem("mitsa005_CCD_Front_PULNIX_13.pictures", 2247 + 17, 1390353, width, height, 3, (int[2]){0x00ff00, 0x0000ff}, 2, uncompressed, compressed, &mp, 0, 1, (double[3]){1, 2.5, 1.5});
 485         fp = fopen("o1.png", "wb");
 486         if (fp == NULL) {
 487                 fprintf(stderr, "Could not open file %s for writing\n", "o1.png");
 488         }
 489         fwrite(mp.buffer, mp.size, sizeof(mp.buffer[0]), fp);
 490         if(mp.buffer)
 491                 free(mp.buffer);
 492 }