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 2 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         int 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 = 1;
  47                 shutter = 1;
  48                 gain = 1;
  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, compile ALE with\n"
  56                 "*** ImageMagick support ('make IMAGEMAGICK=1').\n"
  57                 "*** (To do this, you must have a source distribution\n"
  58                 "*** of ALE, and you must have ImageMagick installed.)\n\n",
  59                 filename);
  60         exit(1);
  61 }
  62
  63 static inline int digest_comment(FILE *f, const char *filename, extended_t *extended) {
  64         int next = '#';
  65         int value;
  66         double fvalue, fvalue2;
  67
  68         while (next != '\n' && next != '\r' && next != EOF) {
  69                 while (next == ' ' || next == '\t' || next == '#') {
  70                         next = fgetc(f);
  71                         if (feof(f))
  72                                 error_ppm(filename);
  73                 }
  74
  75                 if (ungetc(next, f) == EOF) {
  76                         assert(0);
  77                         fprintf(stderr, "Unable to ungetc().");
  78                         exit(1);
  79                 }
  80
  81                 fvalue2 = 1;
  82
  83                 if (extended->is_extended && fscanf(f, "Black-level: %d", &value) == 1)
  84                         extended->black_level = value;
  85                 else if (extended->is_extended && fscanf(f, "ISO: %lf", &fvalue) == 1)
  86                         extended->gain = fvalue / 100;
  87                 else if (extended->is_extended && fscanf(f, "Gain: %lf", &fvalue) == 1)
  88                         extended->gain = fvalue;
  89                 else if (extended->is_extended && fscanf(f, "Aperture: %lf", &fvalue) == 1)
  90                         extended->aperture = fvalue;
  91                 else if (extended->is_extended && fscanf(f, "Shutter: %lf/%lf", &fvalue, &fvalue2) > 0)
  92                         extended->shutter = fvalue / fvalue2;
  93                 else if (next != '\n' && next != '\r' && next != EOF)
  94                         next = fgetc(f);
  95
  96                 next = fgetc(f);
  97         }
  98
  99         return next;
 100 }
 101
 102 static inline void eat_comments(FILE *f, const char *filename, extended_t *extended) {
 103         int next = ' ';
 104
 105         while (next == ' ' || next == '\n' || next == '\t' || next == '#' || next == '\r') {
 106                 next = fgetc(f);
 107                 if (next == '#')
 108                         next = digest_comment(f, filename, extended);
 109                 if (feof(f))
 110                         error_ppm(filename);
 111         }
 112
 113         if (ungetc(next, f) == EOF) {
 114                 assert(0);
 115                 fprintf(stderr, "Unable to ungetc().");
 116                 exit(1);
 117         }
 118 }
 119
 120 static inline int is_eppm(const char *filename) {
 121         char m1, m2, m3, m4;
 122         int n;
 123         extended_t extended;
 124         FILE *f = fopen(filename, "rb");
 125
 126         if (f == NULL)
 127                 return 0;
 128
 129         /* Magic */
 130
 131         eat_comments(f, filename, &extended);   /* XXX - should we eat comments here? */
 132         n = fscanf(f, "%c%c%c%c", &m1, &m2, &m3, &m4);
 133
 134         fclose(f);
 135
 136         if (n != 4 || m1 != 'P' || (m2 != '6' && m2 != '3') || m3 != '#' || m4 != 'E')
 137                 return 0;
 138
 139         return 1;
 140 }
 141
 142 static inline image *read_ppm(const char *filename, exposure *e, unsigned int bayer, int init_reference_gain = 0) {
 143         unsigned int i, j, k;
 144         image *im;
 145         unsigned char m1, m2, val;
 146         int m3, m4;
 147         int ival;
 148         int w, h, mcv;
 149         int n;
 150         struct extended_t extended;
 151         FILE *f = fopen(filename, "rb");
 152
 153         if (f == NULL) {
 154                 fprintf(stderr, "\n\nUnable to open '%s'.\n\n", filename);
 155                 exit(1);
 156         }
 157
 158         assert(f);
 159
 160         /* Magic */
 161
 162         eat_comments(f, filename, &extended);   /* XXX - should we eat comments here? */
 163         n = fscanf(f, "%c%c", &m1, &m2);
 164
 165         if (n != 2 || m1 != 'P' || (m2 != '6' && m2 != '3'))
 166                 error_ppm(filename);
 167
 168         assert(n == 2 && m1 == 'P' && (m2 == '6' || m2 == '3'));
 169
 170         /* Extended flag */
 171
 172         m3 = fgetc(f);
 173
 174         if (m3 == '#') {
 175                 m4 = fgetc(f);
 176                 if (m4 == 'E')
 177                         extended.is_extended = 1;
 178                 else while (m4 != EOF && m4 != '\n' && m4 != '\r')
 179                         m4 = fgetc(f);
 180         } else if (ungetc(m3, f) == EOF) {
 181                 assert(0);
 182                 fprintf(stderr, "Unable to ungetc().");
 183                 exit(1);
 184         }
 185
 186         /* Width */
 187
 188         eat_comments(f, filename, &extended);
 189         n = fscanf(f, " %d", &w);
 190         assert(n == 1);
 191
 192         if (n != 1)
 193                 error_ppm(filename);
 194
 195         /* Height */
 196
 197         eat_comments(f, filename, &extended);
 198         n = fscanf(f, "%d", &h);
 199         assert(n == 1);
 200
 201         if (n != 1)
 202                 error_ppm(filename);
 203
 204         /* Maximum component value */
 205
 206         eat_comments(f, filename, &extended);
 207         n = fscanf(f, "%d", &mcv);
 208         assert(n == 1);
 209         assert(mcv <= 65535 || m2 == '3');
 210
 211         if (n != 1 || (mcv > 65535 && m2 == '6'))
 212                 error_ppm(filename);
 213
 214         /* Make a new image */
 215
 216         if (bayer == IMAGE_BAYER_NONE)
 217                 im = new image_ale_real(h, w, 3, "file", e);
 218         else
 219                 im = new image_bayer_ale_real(h, w, 3, bayer, "file", e);
 220
 221         assert (im);
 222
 223         /* Trailing whitespace */
 224
 225         if (fgetc(f) == EOF) {
 226                 assert(0);
 227                 error_ppm(filename);
 228         }
 229
 230         /* Pixels */
 231
 232         for (i = 0; i < im->height(); i++)
 233         for (j = 0; j < im->width();  j++) {
 234                 pixel p;
 235                 for (k = 0; k < im->depth();  k++) {
 236
 237                         if (m2 == '6') {
 238
 239                                 /* Binary data */
 240
 241                                 n = fscanf(f, "%c", &val);
 242                                 assert (n == 1);
 243
 244                                 if (n != 1)
 245                                         error_ppm(filename);
 246
 247                                 ival = val;
 248
 249                                 if (mcv > 255) {
 250                                         n = fscanf(f, "%c", &val);
 251                                         assert(n == 1);
 252
 253                                         if (n != 1)
 254                                                 error_ppm(filename);
 255
 256                                         ival = (ival << 8) | val;
 257                                 }
 258
 259                         } else {
 260
 261                                 /* ASCII data */
 262
 263                                 eat_comments(f, filename, &extended);
 264
 265                                 n = fscanf(f, "%d", &ival);
 266
 267                                 assert (n == 1);
 268                                 if (n != 1)
 269                                         error_ppm(filename);
 270                         }
 271
 272                         p[k] = (ale_real) (ival - (extended.is_extended ? extended.black_level : 0))
 273                              / (ale_real) (mcv  - (extended.is_extended ? extended.black_level : 0));
 274                 }
 275
 276                 im->set_pixel(i, j, e->linearize(p));
 277         }
 278
 279         /* Handle exposure and gain */
 280
 281         if (extended.is_extended) {
 282                 double combined_gain = (1 / pow(extended.aperture, 2))
 283                                      * extended.shutter
 284                                      * extended.gain;
 285                 if (init_reference_gain)
 286                         exposure::set_gain_reference(combined_gain);
 287                 else
 288                         e->set_gain_multiplier(exposure::get_gain_reference()
 289                                              / combined_gain);
 290         }
 291
 292         /* Done */
 293
 294         fclose(f);
 295
 296         return im;
 297 }
 298
 299 static inline void write_ppm(const char *filename, const image *im, exposure *e,
 300                 unsigned int mcv, int plain, int rezero, int exposure_scale, double nn_defined_radius) {
 301         unsigned int i, j, k;
 302         FILE *f = fopen(filename, "wb");
 303
 304         if (f == NULL) {
 305                 fprintf(stderr, "\n\nUnable to open '%s'.\n\n", filename);
 306                 exit(1);
 307         }
 308
 309         assert(f);
 310
 311         /*
 312          * Output a plain (ASCII) or raw (binary) PPM file
 313          */
 314
 315         /* Magic */
 316
 317         if (plain)
 318                 fprintf(f, "P3 ");
 319         else
 320                 fprintf(f, "P6 ");
 321
 322         /* Width */
 323
 324         fprintf(f, "%d ", im->width());
 325
 326         /* Height */
 327
 328         fprintf(f, "%d ", im->height());
 329
 330         /* Maximum component value */
 331
 332         fprintf(f, "%d\n", mcv);
 333
 334         /* Automatic exposure adjustment information */
 335
 336         ale_real maxval = 1;
 337         ale_real minval = (rezero ? im->minval() : 0);
 338         if (minval > 0)
 339                 minval = 0;
 340         pixel minval_pixel(minval, minval, minval);
 341
 342         if (exposure_scale) {
 343                 ale_real new_maxval = im->maxval();
 344
 345                 if (new_maxval > maxval)
 346                         maxval = new_maxval;
 347         }
 348
 349         /* Pixels */
 350
 351         for (i = 0; i < im->height(); i++)
 352         for (j = 0; j < im->width();  j++) {
 353                 pixel value = im->get_pixel(i, j);
 354
 355                 /*
 356                  * Get nearest-neighbor defined values.
 357                  *
 358                  * XXX: While this implementation is correct, it is inefficient
 359                  * for large radii.  A better implementation would search
 360                  * perimeters of squares of ever-increasing radius, tracking
 361                  * the best-so-far data until the square perimeter exceeded the
 362                  * best-so-far radius.
 363                  */
 364
 365                 for (k = 0; k < 3; k++)
 366                 if (isnan(value[k]))
 367                 for (int radius = 1; radius <= nn_defined_radius; radius++) {
 368                         double nearest_radius_squared = (radius + 1) * (radius + 1);
 369                         for (int ii = -radius; ii <= radius; ii++)
 370                         for (int jj = -radius; jj <= radius; jj++) {
 371                                 if (!im->in_bounds(point(i + ii, j + jj)))
 372                                         continue;
 373                                 if (ii * ii + jj * jj < nearest_radius_squared
 374                                  && finite(im->get_pixel(i + ii, j + jj)[k])) {
 375                                         value[k] = im->get_pixel(i + ii, j + jj)[k];
 376                                         nearest_radius_squared = ii * ii + jj * jj;
 377                                 }
 378                         }
 379                         if (nearest_radius_squared < (radius + 1) * (radius + 1))
 380                                 break;
 381                 }
 382
 383                 pixel exposure_adjust = (value - minval_pixel)
 384                                       / (maxval - minval);
 385                 pixel unlinearized = (e->unlinearize(exposure_adjust)).clamp();
 386
 387                 pixel rescaled = unlinearized * (ale_real) mcv;
 388
 389                 for (k = 0; k < im->depth();  k++) {
 390
 391                         uint16_t output_value = (uint16_t) round(rescaled[k]);
 392
 393                         if (plain) {
 394                                 fprintf(f, "%d ", output_value);
 395                         } else {
 396                                 if (mcv > 255)
 397                                         fprintf(f, "%c", output_value >> 8);
 398                                 fprintf(f, "%c", 0xff & output_value);
 399                         }
 400                 }
 401
 402                 if (plain)
 403                         fprintf(f, "\n");
 404
 405         }
 406
 407         /* Done */
 408
 409         fclose(f);
 410 }
 411
 412 #endif