d2/ppm.h

   1 // Copyright 2002 David Hilvert <dhilvert@auricle.dyndns.org>,
   2 //                              <dhilvert@ugcs.caltech.edu>
   3
   4 /*  This file is part of the Anti-Lamenessing Engine.
   5
   6     The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 3 of the License, or
   9     (at your option) any later version.
  10
  11     The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with the Anti-Lamenessing Engine; if not, write to the Free Software
  18     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19 */
  20
  21 /*
  22  * ppm.h: Read and write PPM files.
  23  */
  24
  25 #ifndef __ppm_h__
  26 #define __ppm_h__
  27
  28 #include "image_ale_real.h"
  29 #include "image_bayer_ale_real.h"
  30 #include "exposure/exposure.h"
  31
  32 /*
  33  * Extended attributes
  34  */
  35
  36 struct extended_t {
  37         int is_extended;
  38         ale_real black_level;
  39         ale_real aperture; /* 1 == f/1.0, 1.4 == f/1.4, etc. */
  40         ale_real shutter;  /* 1 == 1 sec, 0.5 == 1/2 sec, etc. */
  41         ale_real gain;     /* 1 == ISO 100, 2 == ISO 200, etc. */
  42
  43         extended_t() {
  44                 is_extended = 0;
  45                 black_level = 0;
  46                 aperture = 0;
  47                 shutter = 0;
  48                 gain = 0;
  49         }
  50 };
  51
  52 static inline void error_ppm(const char *filename) {
  53         fprintf(stderr,
  54                 "\n\n*** '%s' doesn't look like a PPM file.\n"
  55                 "\n*** (To handle other file types, use a version of ALE with\n"
  56                 "*** ImageMagick support enabled.)\n\n",
  57                 filename);
  58         exit(1);
  59 }
  60
  61 static inline int digest_comment(FILE *f, const char *filename, extended_t *extended) {
  62         int next = '#';
  63         int value;
  64         double fvalue, fvalue2;
  65
  66         while (next != '\n' && next != '\r' && next != EOF) {
  67                 while (next == ' ' || next == '\t' || next == '#') {
  68                         next = fgetc(f);
  69                         if (feof(f))
  70                                 error_ppm(filename);
  71                 }
  72
  73                 if (ungetc(next, f) == EOF) {
  74                         assert(0);
  75                         fprintf(stderr, "Unable to ungetc().");
  76                         exit(1);
  77                 }
  78
  79                 fvalue2 = 1;
  80
  81                 if (extended->is_extended && fscanf(f, "Black-level: %d", &value) == 1)
  82                         extended->black_level = value;
  83                 else if (extended->is_extended && fscanf(f, "ISO: %lf", &fvalue) == 1)
  84                         extended->gain = fvalue / 100;
  85                 else if (extended->is_extended && fscanf(f, "Gain: %lf", &fvalue) == 1)
  86                         extended->gain = fvalue;
  87                 else if (extended->is_extended && fscanf(f, "Aperture: %lf", &fvalue) == 1)
  88                         extended->aperture = fvalue;
  89                 else if (extended->is_extended && fscanf(f, "Shutter: %lf/%lf", &fvalue, &fvalue2) > 0)
  90                         extended->shutter = fvalue / fvalue2;
  91                 else if (next != '\n' && next != '\r' && next != EOF)
  92                         next = fgetc(f);
  93
  94                 next = fgetc(f);
  95         }
  96
  97         return next;
  98 }
  99
 100 static inline void eat_comments(FILE *f, const char *filename, extended_t *extended) {
 101         int next = ' ';
 102
 103         while (next == ' ' || next == '\n' || next == '\t' || next == '#' || next == '\r') {
 104                 next = fgetc(f);
 105                 if (next == '#')
 106                         next = digest_comment(f, filename, extended);
 107                 if (feof(f))
 108                         error_ppm(filename);
 109         }
 110
 111         if (ungetc(next, f) == EOF) {
 112                 assert(0);
 113                 fprintf(stderr, "Unable to ungetc().");
 114                 exit(1);
 115         }
 116 }
 117
 118 static inline int is_eppm(const char *filename) {
 119         char m1, m2, m3, m4;
 120         int n;
 121         extended_t extended;
 122         FILE *f = fopen(filename, "rb");
 123
 124         if (f == NULL)
 125                 return 0;
 126
 127         /* Magic */
 128
 129         eat_comments(f, filename, &extended);   /* XXX - should we eat comments here? */
 130         n = fscanf(f, "%c%c%c%c", &m1, &m2, &m3, &m4);
 131
 132         fclose(f);
 133
 134         if (n != 4 || m1 != 'P' || (m2 != '6' && m2 != '3') || m3 != '#' || m4 != 'E')
 135                 return 0;
 136
 137         return 1;
 138 }
 139
 140 static inline image *read_ppm(const char *filename, exposure *e, unsigned int bayer, int init_reference_gain = 0) {
 141         unsigned int i, j, k;
 142         image *im;
 143         unsigned char m1, m2, val;
 144         int m3, m4;
 145         int ival;
 146         int w, h, mcv;
 147         int n;
 148         struct extended_t extended;
 149         FILE *f = fopen(filename, "rb");
 150
 151         if (f == NULL) {
 152                 fprintf(stderr, "\n\nUnable to open '%s'.\n\n", filename);
 153                 exit(1);
 154         }
 155
 156         assert(f);
 157
 158         /* Magic */
 159
 160         eat_comments(f, filename, &extended);   /* XXX - should we eat comments here? */
 161         n = fscanf(f, "%c%c", &m1, &m2);
 162
 163         if (n != 2 || m1 != 'P' || (m2 != '6' && m2 != '3'))
 164                 error_ppm(filename);
 165
 166         assert(n == 2 && m1 == 'P' && (m2 == '6' || m2 == '3'));
 167
 168         /* Extended flag */
 169
 170         m3 = fgetc(f);
 171
 172         if (m3 == '#') {
 173                 m4 = fgetc(f);
 174                 if (m4 == 'E')
 175                         extended.is_extended = 1;
 176                 else while (m4 != EOF && m4 != '\n' && m4 != '\r')
 177                         m4 = fgetc(f);
 178         } else if (ungetc(m3, f) == EOF) {
 179                 assert(0);
 180                 fprintf(stderr, "Unable to ungetc().");
 181                 exit(1);
 182         }
 183
 184         /* Width */
 185
 186         eat_comments(f, filename, &extended);
 187         n = fscanf(f, " %d", &w);
 188         assert(n == 1);
 189
 190         if (n != 1)
 191                 error_ppm(filename);
 192
 193         /* Height */
 194
 195         eat_comments(f, filename, &extended);
 196         n = fscanf(f, "%d", &h);
 197         assert(n == 1);
 198
 199         if (n != 1)
 200                 error_ppm(filename);
 201
 202         /* Maximum component value */
 203
 204         eat_comments(f, filename, &extended);
 205         n = fscanf(f, "%d", &mcv);
 206         assert(n == 1);
 207         assert(mcv <= 65535 || m2 == '3');
 208
 209         if (extended.black_level == 0) {
 210                 extended.black_level = e->get_black_level();
 211         } else {
 212                 extended.black_level /= mcv;
 213         }
 214
 215         if (n != 1 || (mcv > 65535 && m2 == '6'))
 216                 error_ppm(filename);
 217
 218         /* Make a new image */
 219
 220         if (bayer == IMAGE_BAYER_NONE)
 221                 im = new image_ale_real(h, w, 3, "file", e);
 222         else
 223                 im = new image_bayer_ale_real(h, w, 3, bayer, "file", e);
 224
 225         assert (im);
 226
 227         /* Trailing whitespace */
 228
 229         if (fgetc(f) == EOF) {
 230                 assert(0);
 231                 error_ppm(filename);
 232         }
 233
 234         /* Pixels */
 235
 236         for (i = 0; i < im->height(); i++)
 237         for (j = 0; j < im->width();  j++) {
 238                 pixel p;
 239                 for (k = 0; k < im->depth();  k++) {
 240
 241                         if (m2 == '6') {
 242
 243                                 /* Binary data */
 244
 245                                 n = fscanf(f, "%c", &val);
 246                                 assert (n == 1);
 247
 248                                 if (n != 1)
 249                                         error_ppm(filename);
 250
 251                                 ival = val;
 252
 253                                 if (mcv > 255) {
 254                                         n = fscanf(f, "%c", &val);
 255                                         assert(n == 1);
 256
 257                                         if (n != 1)
 258                                                 error_ppm(filename);
 259
 260                                         ival = (ival << 8) | val;
 261                                 }
 262
 263                         } else {
 264
 265                                 /* ASCII data */
 266
 267                                 eat_comments(f, filename, &extended);
 268
 269                                 n = fscanf(f, "%d", &ival);
 270
 271                                 assert (n == 1);
 272                                 if (n != 1)
 273                                         error_ppm(filename);
 274                         }
 275
 276                         p[k] = ((ale_real) ival) / ((ale_real) mcv);
 277
 278                 }
 279
 280                 pixel p_linear = (e->linearize(p) - e->get_multiplier() * extended.black_level)
 281                                / (1 - extended.black_level);
 282
 283                 im->set_pixel(i, j, p_linear);
 284         }
 285
 286         /* Handle exposure and gain */
 287
 288         if (extended.is_extended) {
 289                 if (extended.aperture != 0
 290                  || extended.shutter != 0
 291                  || extended.gain != 0) {
 292
 293                         if (extended.aperture == 0)
 294                                 extended.aperture = 1;
 295                         if (extended.shutter == 0)
 296                                 extended.shutter = 1;
 297                         if (extended.gain == 0)
 298                                 extended.gain = 1;
 299
 300                         double combined_gain = (1 / pow(extended.aperture, 2))
 301                                              * extended.shutter
 302                                              * extended.gain;
 303
 304                         if (init_reference_gain)
 305                                 exposure::set_gain_reference(combined_gain);
 306                         else
 307                                 e->set_gain_multiplier(exposure::get_gain_reference()
 308                                                      / combined_gain);
 309                 }
 310         }
 311
 312         /* Done */
 313
 314         fclose(f);
 315
 316         return im;
 317 }
 318
 319 static inline void write_ppm(const char *filename, const image *im, exposure *e,
 320                 unsigned int mcv, int plain, int rezero, int exposure_scale, double nn_defined_radius) {
 321         unsigned int i, j, k;
 322         FILE *f = fopen(filename, "wb");
 323
 324         if (f == NULL) {
 325                 fprintf(stderr, "\n\nUnable to open '%s'.\n\n", filename);
 326                 exit(1);
 327         }
 328
 329         assert(f);
 330
 331         /*
 332          * Output a plain (ASCII) or raw (binary) PPM file
 333          */
 334
 335         /* Magic */
 336
 337         if (plain)
 338                 fprintf(f, "P3 ");
 339         else
 340                 fprintf(f, "P6 ");
 341
 342         /* Width */
 343
 344         fprintf(f, "%d ", im->width());
 345
 346         /* Height */
 347
 348         fprintf(f, "%d ", im->height());
 349
 350         /* Maximum component value */
 351
 352         fprintf(f, "%d\n", mcv);
 353
 354         /* Automatic exposure adjustment information */
 355
 356         ale_real maxval = 1;
 357         ale_real minval = (rezero ? im->minval() : 0);
 358         if (minval > 0)
 359                 minval = 0;
 360         pixel minval_pixel(minval, minval, minval);
 361
 362         if (exposure_scale) {
 363                 ale_real new_maxval = im->maxval();
 364
 365                 if (new_maxval > maxval)
 366                         maxval = new_maxval;
 367         }
 368
 369         /* Pixels */
 370
 371         for (i = 0; i < im->height(); i++)
 372         for (j = 0; j < im->width();  j++) {
 373                 pixel value = im->get_pixel(i, j);
 374
 375                 /*
 376                  * Get nearest-neighbor defined values.
 377                  *
 378                  * XXX: While this implementation is correct, it is inefficient
 379                  * for large radii.  A better implementation would search
 380                  * perimeters of squares of ever-increasing radius, tracking
 381                  * the best-so-far data until the square perimeter exceeded the
 382                  * best-so-far radius.
 383                  */
 384
 385                 for (k = 0; k < 3; k++)
 386                 if (isnan(value[k]))
 387                 for (int radius = 1; radius <= nn_defined_radius; radius++) {
 388                         double nearest_radius_squared = (radius + 1) * (radius + 1);
 389                         for (int ii = -radius; ii <= radius; ii++)
 390                         for (int jj = -radius; jj <= radius; jj++) {
 391                                 if (!im->in_bounds(point(i + ii, j + jj)))
 392                                         continue;
 393                                 if (ii * ii + jj * jj < nearest_radius_squared
 394                                  && finite(im->get_pixel(i + ii, j + jj)[k])) {
 395                                         value[k] = im->get_pixel(i + ii, j + jj)[k];
 396                                         nearest_radius_squared = ii * ii + jj * jj;
 397                                 }
 398                         }
 399                         if (nearest_radius_squared < (radius + 1) * (radius + 1))
 400                                 break;
 401                 }
 402
 403                 pixel exposure_adjust = (value - minval_pixel)
 404                                       / (maxval - minval);
 405                 pixel unlinearized = (e->unlinearize(exposure_adjust)).clamp();
 406
 407                 pixel rescaled = unlinearized * (ale_real) mcv;
 408
 409                 for (k = 0; k < im->depth();  k++) {
 410
 411                         uint16_t output_value = (uint16_t) round(rescaled[k]);
 412
 413                         if (plain) {
 414                                 fprintf(f, "%d ", output_value);
 415                         } else {
 416                                 if (mcv > 255)
 417                                         fprintf(f, "%c", output_value >> 8);
 418                                 fprintf(f, "%c", 0xff & output_value);
 419                         }
 420                 }
 421
 422                 if (plain)
 423                         fprintf(f, "\n");
 424
 425         }
 426
 427         /* Done */
 428
 429         fclose(f);
 430 }
 431
 432 #endif