d2/image.h

   1 // Copyright 2002, 2003, 2004 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  * image.h: Abstract base class for the internal representations of images used
  23  * by ALE.
  24  */
  25
  26 #ifndef __image_h__
  27 #define __image_h__
  28
  29 #include "point.h"
  30 #include "pixel.h"
  31 #include "exposure/exposure.h"
  32
  33 #define IMAGE_BAYER_NONE 0
  34
  35 /*
  36  * This constant indicates that some other default value should be filled in.
  37  */
  38
  39 #define IMAGE_BAYER_DEFAULT 0x8
  40
  41 /*
  42  * Do not change these values without inspecting
  43  * image_bayer_ale_real::r_*_offset().
  44  */
  45 #define IMAGE_BAYER_RGBG 0x4 /* binary 100 */
  46 #define IMAGE_BAYER_GBGR 0x5 /* binary 101 */
  47 #define IMAGE_BAYER_GRGB 0x6 /* binary 110 */
  48 #define IMAGE_BAYER_BGRG 0x7 /* binary 111 */
  49
  50 class image : protected exposure::listener {
  51 protected:
  52         unsigned int _dimx, _dimy, _depth;
  53         point _offset;
  54         char *name;
  55         mutable exposure *_exp;
  56         unsigned int bayer;
  57 private:
  58         /*
  59          * Memoized function variables.  We may want to change these even when
  60          * *this is constant.
  61          */
  62         mutable int _apm_memo;
  63         mutable ale_real _apm;
  64         mutable int _accm_memo;
  65         mutable pixel _accm;
  66         mutable int _acm_memo;
  67         mutable pixel _acm;
  68
  69         void avg_channel_clamped_magnitude_memo() const {
  70                 unsigned int i, j, k;
  71                 pixel_accum accumulator;
  72                 pixel_accum divisor;
  73
  74                 if (_accm_memo)
  75                         return;
  76
  77                 _accm_memo = 1;
  78
  79                 accumulator = pixel_accum(0, 0, 0);
  80
  81                 for (i = 0; i < _dimy;  i++)
  82                 for (j = 0; j < _dimx;  j++) {
  83                         pixel value = get_pixel(i, j);
  84
  85                         for (k = 0; k < _depth; k++)
  86                                 if (finite(value[k])) {
  87                                         if (value[k] > 1)
  88                                                 value[k] = 1;
  89                                         if (value[k] < 0)
  90                                                 value[k] = 0;
  91                                         accumulator[k] += value[k];
  92                                         divisor[k] += 1;
  93                                 }
  94                 }
  95
  96                 accumulator /= divisor;
  97
  98                 _accm = accumulator;
  99         }
 100
 101         void avg_channel_magnitude_memo() const {
 102                 unsigned int i, j, k;
 103                 pixel_accum accumulator;
 104                 pixel_accum divisor;
 105
 106                 if (_acm_memo)
 107                         return;
 108
 109                 _acm_memo = 1;
 110
 111                 accumulator = pixel_accum(0, 0, 0);
 112
 113                 for (i = 0; i < _dimy;  i++)
 114                 for (j = 0; j < _dimx;  j++) {
 115                         pixel value = get_pixel(i, j);
 116
 117                         for (k = 0; k < _depth; k++)
 118                                 if (finite(value[k])) {
 119                                         accumulator[k] += value[k];
 120                                         divisor[k] += 1;
 121                                 }
 122                 }
 123
 124                 accumulator /= divisor;
 125
 126                 _acm = accumulator;
 127         }
 128
 129 protected:
 130         void image_updated() {
 131                 _apm_memo = 0;
 132                 _acm_memo = 0;
 133                 _accm_memo = 0;
 134         }
 135
 136 public:
 137         image (unsigned int dimy, unsigned int dimx, unsigned int depth,
 138                         char *name = "anonymous", exposure *_exp = NULL,
 139                         unsigned int bayer = IMAGE_BAYER_NONE) {
 140
 141                 assert (depth == 3);
 142                 _depth = 3;
 143
 144                 _dimx = dimx;
 145                 _dimy = dimy;
 146                 _offset = point(0, 0);
 147                 _apm_memo = 0;
 148                 _acm_memo = 0;
 149                 _accm_memo = 0;
 150                 this->name = name;
 151                 this->_exp = _exp;
 152                 this->bayer = bayer;
 153
 154                 if (_exp != NULL)
 155                         _exp->add_listener(this, name);
 156         }
 157
 158         unsigned int get_bayer() const {
 159                 return bayer;
 160         }
 161
 162         virtual char get_channels(int i, int j) const {
 163                 return 0x7;
 164         }
 165
 166         double storage_size() const {
 167                 if (bayer != IMAGE_BAYER_NONE)
 168                         return _dimx * _dimy * sizeof(ale_real);
 169
 170                 return 3 * _dimx * _dimy * sizeof(ale_real);
 171         }
 172
 173         exposure &exp() const {
 174                 return *_exp;
 175         }
 176
 177         point offset() const {
 178                 return _offset;
 179         }
 180
 181         void set_offset(int i, int j) {
 182                 _offset[0] = i;
 183                 _offset[1] = j;
 184         }
 185
 186         void set_offset(point p) {
 187                 _offset = p;
 188         }
 189
 190         unsigned int width() const {
 191                 return _dimx;
 192         }
 193
 194         unsigned int height() const {
 195                 return _dimy;
 196         }
 197
 198         unsigned int depth() const {
 199                 return _depth;
 200         }
 201
 202         virtual spixel &pix(unsigned int y, unsigned int x) = 0;
 203
 204         virtual const spixel &pix(unsigned int y, unsigned int x) const = 0;
 205
 206         virtual ale_real &chan(unsigned int y, unsigned int x, unsigned int k) = 0;
 207
 208         virtual const ale_real &chan(unsigned int y, unsigned int x, unsigned int k) const = 0;
 209
 210         virtual void set_pixel(unsigned int y, unsigned int x, spixel p) {
 211                 pix(y, x) = p;
 212         }
 213
 214         virtual pixel get_pixel(unsigned int y, unsigned int x) const {
 215                 return ((const image *)this)->pix(y, x);
 216         }
 217
 218         virtual pixel get_raw_pixel(unsigned int y, unsigned int x) const {
 219                 return ((const image *)this)->get_pixel(y, x);
 220         }
 221
 222         ale_real maxval() const {
 223                 ale_real result = get_pixel(0, 0)[0];
 224
 225                 for (unsigned int i = 0; i < _dimy; i++)
 226                 for (unsigned int j = 0; j < _dimx; j++) {
 227                         pixel p = get_pixel(i, j);
 228
 229                         for (unsigned int k = 0; k < _depth; k++)
 230                                 if (p[k] > result || !finite(result))
 231                                         result = p[k];
 232                 }
 233
 234                 return result;
 235         }
 236
 237         ale_real minval() const {
 238                 ale_real result = get_pixel(0, 0)[0];
 239
 240                 for (unsigned int i = 0; i < _dimy; i++)
 241                 for (unsigned int j = 0; j < _dimx; j++) {
 242                         pixel p = get_pixel(i, j);
 243
 244                         for (unsigned int k = 0; k < _depth; k++)
 245                                 if (p[k] < result || !finite(result))
 246                                         result = p[k];
 247                 }
 248
 249                 return result;
 250         }
 251
 252         /*
 253          * Get the maximum difference among adjacent pixels.
 254          */
 255
 256         pixel get_max_diff(unsigned int i, unsigned int j) const {
 257                 assert(i <= _dimy - 1);
 258                 assert(j <= _dimx - 1);
 259
 260                 pixel max = get_pixel(i, j), min = get_pixel(i, j);
 261
 262                 for (int ii = -1; ii <= 1; ii++)
 263                 for (int jj = -1; jj <= 1; jj++) {
 264                         int iii = i + ii;
 265                         int jjj = j + jj;
 266
 267                         if (iii < 0)
 268                                 continue;
 269                         if (jjj < 0)
 270                                 continue;
 271                         if ((unsigned int) iii > _dimy - 1)
 272                                 continue;
 273                         if ((unsigned int) jjj > _dimx - 1)
 274                                 continue;
 275
 276                         pixel p = get_pixel(iii, jjj);
 277
 278                         for (int d = 0; d < 3; d++) {
 279                                 if (p[d] > max[d])
 280                                         max[d] = p[d];
 281                                 if (p[d] < min[d])
 282                                         min[d] = p[d];
 283                         }
 284                 }
 285
 286                 return max - min;
 287         }
 288
 289         pixel get_max_diff(point x) const {
 290                 assert (x[0] >= 0);
 291                 assert (x[1] >= 0);
 292                 assert (x[0] <= _dimy - 1);
 293                 assert (x[1] <= _dimx - 1);
 294
 295                 unsigned int i = (unsigned int) round(x[0]);
 296                 unsigned int j = (unsigned int) round(x[1]);
 297
 298                 return get_max_diff(i, j);
 299         }
 300
 301         int in_bounds(point x) const {
 302                 if (x[0] < 0
 303                  || x[1] < 0
 304                  || x[0] > height() - 1
 305                  || x[1] > width() - 1)
 306                         return 0;
 307                 if (!x.defined())
 308                         return 0;
 309                 return 1;
 310         }
 311
 312         /*
 313          * Get a color value at a given position using bilinear interpolation between the
 314          * four nearest pixels.
 315          */
 316         pixel get_bl(point x, int defined = 0) const {
 317
 318                 pixel result;
 319
 320                 assert (x[0] >= 0);
 321                 assert (x[1] >= 0);
 322                 assert (x[0] <= _dimy - 1);
 323                 assert (x[1] <= _dimx - 1);
 324
 325                 int lx = (int) floor(x[1]);
 326                 int hx = (int) floor(x[1]) + 1;
 327                 int ly = (int) floor(x[0]);
 328                 int hy = (int) floor(x[0]) + 1;
 329
 330                 pixel neighbor[4];
 331                 ale_pos factor[4];
 332
 333                 neighbor[0] = get_pixel(ly, lx);
 334                 neighbor[1] = get_pixel(hy % _dimy, lx);
 335                 neighbor[2] = get_pixel(hy % _dimy, hx % _dimx);
 336                 neighbor[3] = get_pixel(ly, hx % _dimx);
 337
 338                 factor[0] = (hx - x[1]) * (hy - x[0]);
 339                 factor[1] = (hx - x[1]) * (x[0] - ly);
 340                 factor[2] = (x[1] - lx) * (x[0] - ly);
 341                 factor[3] = (x[1] - lx) * (hy - x[0]);
 342
 343                 /*
 344                  * Use bilinear and/or geometric interpolation
 345                  */
 346
 347                 if (defined == 0) {
 348                         result = pixel(0, 0, 0);
 349
 350                         for (int n = 0; n < 4; n++)
 351                                 result += factor[n] * neighbor[n];
 352                 } else {
 353                         result = pixel(1, 1, 1);
 354
 355                         for (int n = 0; n < 4; n++)
 356                                 result *= ppow(neighbor[n], factor[n]);
 357                 }
 358
 359                 return result;
 360         }
 361
 362         pixel get_scaled_bl(point x, ale_pos f, int defined = 0) const {
 363                 point scaled(
 364                                 x[0]/f <= height() - 1
 365                                        ? (x[0]/f)
 366                                        : (height() - 1),
 367                                 x[1]/f <= width() - 1
 368                                        ? (x[1]/f)
 369                                        : (width() - 1));
 370
 371                 return get_bl(scaled, defined);
 372         }
 373
 374
 375         /*
 376          * Make a new image suitable for receiving scaled values.
 377          */
 378         virtual image *scale_generator(int height, int width, int depth, char *name) const = 0;
 379
 380         /*
 381          * Generate an image of medians within a given radius
 382          */
 383
 384         image *medians(int radius) const {
 385
 386                 assert (radius >= 0);
 387
 388                 image *is = scale_generator(height(), width(), depth(), "median");
 389                 assert(is);
 390
 391                 for (unsigned int i = 0; i < height(); i++)
 392                 for (unsigned int j = 0; j < width(); j++) {
 393
 394                         std::vector<ale_real> p[3];
 395
 396                         for (int ii = -radius; ii <= radius; ii++)
 397                         for (int jj = -radius; jj <= radius; jj++) {
 398                                 int iii = i + ii;
 399                                 int jjj = j + jj;
 400
 401                                 if (in_bounds(point(iii, jjj)))
 402                                         for (int k = 0; k < 3; k++)
 403                                                 if (finite(get_pixel(iii, jjj)[k]))
 404                                                         p[k].push_back(get_pixel(iii, jjj)[k]);
 405                         }
 406
 407                         is->pix(i, j) = d2::pixel::undefined();
 408
 409                         for (int k = 0; k < 3; k++) {
 410                                 std::sort(p[k].begin(), p[k].end());
 411
 412                                 unsigned int pkc = p[k].size();
 413
 414                                 if (pkc == 0)
 415                                         continue;
 416
 417                                 if (pkc % 2 == 0)
 418                                         is->chan(i, j, k) = (p[k][pkc / 2] + p[k][pkc / 2 - 1]) / 2;
 419                                 else
 420                                         is->chan(i, j, k) = p[k][pkc / 2];
 421                         }
 422                 }
 423
 424                 return is;
 425         }
 426
 427         /*
 428          * Generate an image of differences of the first channel.  The first
 429          * coordinate differences are stored in the first channel, second in the
 430          * second channel.
 431          */
 432
 433         image *fcdiffs() const {
 434                 image *is = scale_generator(height(), width(), depth(), "diff");
 435
 436                 assert(is);
 437
 438                 for (unsigned int i = 0; i < height(); i++)
 439                 for (unsigned int j = 0; j < width(); j++) {
 440
 441                         if (i + 1 < height()
 442                          && i > 0
 443                          && !finite(chan(i, j, 0))) {
 444
 445                                 is->chan(i, j, 0) = (chan(i + 1, j, 0) - chan(i - 1, j, 0)) / 2;
 446
 447                         } else if (i + 1 < height()
 448                                 && i > 0
 449                                 && finite(chan(i + 1, j, 0))
 450                                 && finite(chan(i - 1, j, 0))) {
 451
 452                                 is->chan(i, j, 0) = ((chan(i, j, 0) - chan(i - 1, j, 0))
 453                                                    + (chan(i + 1, j, 0) - chan(i, j, 0))) / 2;
 454
 455                         } else if (i + 1 < height()
 456                                 && finite(chan(i + 1, j, 0))) {
 457
 458                                 is->chan(i, j, 0) = chan(i + 1, j, 0) - chan(i, j, 0);
 459
 460                         } else if (i > 0
 461                                 && finite(chan(i - 1, j, 0))) {
 462
 463                                 is->chan(i, j, 0) = chan(i, j, 0) - chan(i - 1, j, 0);
 464
 465                         } else {
 466                                 is->chan(i, j, 0) = 0;
 467                         }
 468
 469                         if (j + 1 < width()
 470                          && j > 0
 471                          && !finite(chan(i, j, 0))) {
 472
 473                                 is->chan(i, j, 1) = (chan(i, j + 1, 0) - chan(i, j - 1, 0)) / 2;
 474
 475                         } else if (j + 1 < width()
 476                                 && j > 0
 477                                 && finite(chan(i, j + 1, 0))
 478                                 && finite(chan(i, j - 1, 0))) {
 479
 480                                 is->chan(i, j, 1) = ((chan(i, j, 0) - chan(i, j - 1, 0))
 481                                                    + (chan(i, j + 1, 0) - chan(i, j, 0))) / 2;
 482
 483                         } else if (j + 1 < width() && finite(chan(i, j + 1, 0))) {
 484
 485                                 is->chan(i, j, 1) = chan(i, j + 1, 0) - chan(i, j, 0);
 486
 487                         } else if (j > 0 && finite(chan(i, j - 1, 0))) {
 488
 489                                 is->chan(i, j, 1) = chan(i, j, 0) - chan(i, j - 1, 0);
 490
 491                         } else {
 492                                 is->chan(i, j, 1) = 0;
 493                         }
 494                 }
 495
 496                 return is;
 497         }
 498
 499         /*
 500          * Generate an image of median (within a given radius) difference of the
 501          * first channel.
 502          */
 503
 504         image *fcdiff_median(int radius) const {
 505                 image *diff = fcdiffs();
 506
 507                 assert(diff);
 508
 509                 image *median = diff->medians(radius);
 510
 511                 assert(median);
 512
 513                 delete diff;
 514
 515                 return median;
 516         }
 517
 518         /*
 519          * Scale by half.  We use the following filter:
 520          *
 521          *      1/16    1/8     1/16
 522          *      1/8     1/4     1/8
 523          *      1/16    1/8     1/16
 524          *
 525          * At the edges, these values are normalized so that the sum of the
 526          * weights of contributing pixels is 1.
 527          */
 528         image *scale_by_half(char *name) const {
 529                 ale_pos f = 0.5;
 530
 531                 image *is = scale_generator(
 532                         (int) floor(height() * f),
 533                         (int) floor(width()  * f), depth(), name);
 534
 535                 assert(is);
 536
 537                 for (unsigned int i = 0; i < is->height(); i++)
 538                 for (unsigned int j = 0; j < is->width(); j++)
 539
 540                         is->set_pixel(i, j,
 541
 542                               ( ( ((i > 0 && j > 0)
 543                                     ? get_pixel(2 * i - 1, 2 * j - 1) * (ale_real) 0.0625
 544                                     : pixel(0, 0, 0))
 545                                 + ((i > 0)
 546                                     ? get_pixel(2 * i - 1, 2 * j) * 0.125
 547                                     : pixel(0, 0, 0))
 548                                 + ((i > 0 && j < is->width() - 1)
 549                                     ? get_pixel(2 * i - 1, 2 * j + 1) * 0.0625
 550                                     : pixel(0, 0, 0))
 551                                 + ((j > 0)
 552                                     ? get_pixel(2 * i, 2 * j - 1) * 0.125
 553                                     : pixel(0, 0, 0))
 554                                 + get_pixel(2 * i, 2 * j) * 0.25
 555                                 + ((j < is->width() - 1)
 556                                     ? get_pixel(2 * i, 2 * j + 1) * 0.125
 557                                     : pixel(0, 0, 0))
 558                                 + ((i < is->height() - 1 && j > 0)
 559                                     ? get_pixel(2 * i + 1, 2 * j - 1) * 0.0625
 560                                     : pixel(0, 0, 0))
 561                                 + ((i < is->height() - 1)
 562                                     ? get_pixel(2 * i + 1, 2 * j) * 0.125
 563                                     : pixel(0, 0, 0))
 564                                 + ((i < is->height() && j < is->width() - 1)
 565                                     ? get_pixel(2 * i + 1, 2 * j + 1) * 0.0625
 566                                     : pixel(0, 0, 0)))
 567
 568                              /
 569
 570                                 ( ((i > 0 && j > 0)
 571                                     ? 0.0625
 572                                     : 0)
 573                                 + ((i > 0)
 574                                     ? 0.125
 575                                     : 0)
 576                                 + ((i > 0 && j < is->width() - 1)
 577                                     ? 0.0625
 578                                     : 0)
 579                                 + ((j > 0)
 580                                     ? 0.125
 581                                     : 0)
 582                                 + 0.25
 583                                 + ((j < is->width() - 1)
 584                                     ? 0.125
 585                                     : 0)
 586                                 + ((i < is->height() - 1 && j > 0)
 587                                     ? 0.0625
 588                                     : 0)
 589                                 + ((i < is->height() - 1)
 590                                     ? 0.125
 591                                     : 0)
 592                                 + ((i < is->height() && j < is->width() - 1)
 593                                     ? 0.0625
 594                                     : 0) ) ) );
 595
 596                 is->_offset = point(_offset[0] * f, _offset[1] * f);
 597
 598                 return is;
 599         }
 600
 601         /*
 602          * Scale by half.  This function uses externally-provided weights,
 603          * multiplied by the following filter:
 604          *
 605          *      1/16    1/8     1/16
 606          *      1/8     1/4     1/8
 607          *      1/16    1/8     1/16
 608          *
 609          * Values are normalized so that the sum of the weights of contributing
 610          * pixels is 1.
 611          */
 612         image *scale_by_half(const image *weights, char *name) const {
 613
 614                 if (weights == NULL)
 615                         return scale_by_half(name);
 616
 617                 ale_pos f = 0.5;
 618
 619                 image *is = scale_generator(
 620                         (int) floor(height() * f),
 621                         (int) floor(width()  * f), depth(), name);
 622
 623                 assert(is);
 624
 625                 for (unsigned int i = 0; i < is->height(); i++)
 626                 for (unsigned int j = 0; j < is->width(); j++) {
 627
 628                         pixel value = pixel
 629
 630                               ( ( ((i > 0 && j > 0)
 631                                     ? get_pixel(2 * i - 1, 2 * j - 1)
 632                                     * weights->get_pixel(2 * i - 1, 2 * j - 1)
 633                                     * (ale_real) 0.0625
 634                                     : pixel(0, 0, 0))
 635                                 + ((i > 0)
 636                                     ? get_pixel(2 * i - 1, 2 * j)
 637                                     * weights->get_pixel(2 * i - 1, 2 * j)
 638                                     * 0.125
 639                                     : pixel(0, 0, 0))
 640                                 + ((i > 0 && j < is->width() - 1)
 641                                     ? get_pixel(2 * i - 1, 2 * j + 1)
 642                                     * weights->get_pixel(2 * i - 1, 2 * j + 1)
 643                                     * 0.0625
 644                                     : pixel(0, 0, 0))
 645                                 + ((j > 0)
 646                                     ? get_pixel(2 * i, 2 * j - 1)
 647                                     * weights->get_pixel(2 * i, 2 * j - 1)
 648                                     * 0.125
 649                                     : pixel(0, 0, 0))
 650                                 + get_pixel(2 * i, 2 * j)
 651                                     * weights->get_pixel(2 * i, 2 * j)
 652                                     * 0.25
 653                                 + ((j < is->width() - 1)
 654                                     ? get_pixel(2 * i, 2 * j + 1)
 655                                     * weights->get_pixel(2 * i, 2 * j + 1)
 656                                     * 0.125
 657                                     : pixel(0, 0, 0))
 658                                 + ((i < is->height() - 1 && j > 0)
 659                                     ? get_pixel(2 * i + 1, 2 * j - 1)
 660                                     * weights->get_pixel(2 * i + 1, 2 * j - 1)
 661                                     * 0.0625
 662                                     : pixel(0, 0, 0))
 663                                 + ((i < is->height() - 1)
 664                                     ? get_pixel(2 * i + 1, 2 * j)
 665                                     * weights->get_pixel(2 * i + 1, 2 * j)
 666                                     * 0.125
 667                                     : pixel(0, 0, 0))
 668                                 + ((i < is->height() && j < is->width() - 1)
 669                                     ? get_pixel(2 * i + 1, 2 * j + 1)
 670                                     * weights->get_pixel(2 * i + 1, 2 * j + 1)
 671                                     * 0.0625
 672                                     : pixel(0, 0, 0)))
 673
 674                              /
 675
 676                                 ( ((i > 0 && j > 0)
 677                                     ? weights->get_pixel(2 * i - 1, 2 * j - 1)
 678                                     * 0.0625
 679                                     : pixel(0, 0, 0))
 680                                 + ((i > 0)
 681                                     ? weights->get_pixel(2 * i - 1, 2 * j)
 682                                     * 0.125
 683                                     : pixel(0, 0, 0))
 684                                 + ((i > 0 && j < is->width() - 1)
 685                                     ? weights->get_pixel(2 * i - 1, 2 * j + 1)
 686                                     * 0.0625
 687                                     : pixel(0, 0, 0))
 688                                 + ((j > 0)
 689                                     ? weights->get_pixel(2 * i, 2 * j - 1)
 690                                     * 0.125
 691                                     : pixel(0, 0, 0))
 692                                 + weights->get_pixel(2 * i, 2 * j)
 693                                 * 0.25
 694                                 + ((j < is->width() - 1)
 695                                     ? weights->get_pixel(2 * i, 2 * j + 1)
 696                                     * 0.125
 697                                     : pixel(0, 0, 0))
 698                                 + ((i < is->height() - 1 && j > 0)
 699                                     ? weights->get_pixel(2 * i + 1, 2 * j - 1)
 700                                     * 0.0625
 701                                     : pixel(0, 0, 0))
 702                                 + ((i < is->height() - 1)
 703                                     ? weights->get_pixel(2 * i + 1, 2 * j)
 704                                     * 0.125
 705                                     : pixel(0, 0, 0))
 706                                 + ((i < is->height() && j < is->width() - 1)
 707                                     ? weights->get_pixel(2 * i + 1, 2 * j + 1)
 708                                     * 0.0625
 709                                     : pixel(0, 0, 0)) ) );
 710
 711                         for (int k = 0; k < 3; k++)
 712                                 if (!finite(value[k]))
 713                                         value[k] = 0;
 714
 715                         is->set_pixel(i, j, value);
 716                 }
 717
 718                 is->_offset = point(_offset[0] * f, _offset[1] * f);
 719
 720                 return is;
 721         }
 722
 723         /*
 724          * Return an image scaled by some factor != 1.0, using bilinear
 725          * interpolation.
 726          */
 727         image *scale(ale_pos f, char *name, int defined = 0) const {
 728
 729                 /*
 730                  * We probably don't want to scale images by a factor of 1.0,
 731                  * or by non-positive values.
 732                  */
 733                 assert (f != 1.0 && f > 0);
 734
 735                 if (f > 1.0) {
 736                         image *is = scale_generator(
 737                                 (int) floor(height() * f),
 738                                 (int) floor(width()  * f), depth(), name);
 739
 740                         assert(is);
 741
 742                         unsigned int i, j, k;
 743
 744                         for (i = 0; i < is->height(); i++)
 745                         for (j = 0; j < is->width(); j++)
 746                         for (k = 0; k < is->depth(); k++)
 747                                 is->set_pixel(i, j,
 748                                         get_scaled_bl(point(i, j), f, defined));
 749
 750                         is->_offset = point(_offset[0] * f, _offset[1] * f);
 751
 752                         return is;
 753                 } else if (f == 0.5) {
 754                         return scale_by_half(name);
 755                 } else {
 756                         image *is = scale(2*f, name);
 757                         image *result = is->scale(0.5, name);
 758                         delete is;
 759                         return result;
 760                 }
 761
 762         }
 763
 764         /*
 765          * Scale an image definition array by 1/2.
 766          *
 767          * ALE considers an image definition array as a special kind of image
 768          * weight array (typedefs of which should appear below the definition
 769          * of this class).  ALE uses nonzero pixel values to mean 'defined' and
 770          * zero values to mean 'undefined'.  Through this interpretation, the
 771          * image weight array implementation that ALE uses allows image weight
 772          * arrays to also serve as image definition arrays.
 773          *
 774          * Whereas scaling of image weight arrays is not generally obvious in
 775          * either purpose or method, ALE requires that image definition arrays
 776          * be scalable.  (Note that in the special case where weight is treated
 777          * as certainty, using a geometric mean is probably correct.)
 778          *
 779          * We currently use a geometric mean to implement scaling of
 780          * definition arrays.
 781          */
 782
 783         image *defined_scale_by_half(char *name) const {
 784                 ale_pos f = 0.5;
 785
 786                 image *is = scale_generator(
 787                         (int) floor(height() * f),
 788                         (int) floor(width()  * f), depth(), name);
 789
 790                 assert(is);
 791
 792                 for (unsigned int i = 0; i < is->height(); i++)
 793                 for (unsigned int j = 0; j < is->width(); j++) {
 794
 795                         pixel value = pixel
 796
 797                               ( ( ((i > 0 && j > 0)
 798                                     ? ppow(get_pixel(2 * i - 1, 2 * j - 1), 0.0625)
 799                                     : pixel(0, 0, 0))
 800                                 * ((i > 0)
 801                                     ? ppow(get_pixel(2 * i - 1, 2 * j), 0.125)
 802                                     : pixel(0, 0, 0))
 803                                 * ((i > 0 && j < is->width() - 1)
 804                                     ? ppow(get_pixel(2 * i - 1, 2 * j + 1), 0.0625)
 805                                     : pixel(0, 0, 0))
 806                                 * ((j > 0)
 807                                     ? ppow(get_pixel(2 * i, 2 * j - 1), 0.125)
 808                                     : pixel(0, 0, 0))
 809                                 * ppow(get_pixel(2 * i, 2 * j), 0.25)
 810                                 * ((j < is->width() - 1)
 811                                     ? ppow(get_pixel(2 * i, 2 * j + 1), 0.125)
 812                                     : pixel(0, 0, 0))
 813                                 * ((i < is->height() - 1 && j > 0)
 814                                     ? ppow(get_pixel(2 * i + 1, 2 * j - 1), 0.0625)
 815                                     : pixel(0, 0, 0))
 816                                 * ((i < is->height() - 1)
 817                                     ? ppow(get_pixel(2 * i + 1, 2 * j), 0.125)
 818                                     : pixel(0, 0, 0))
 819                                 * ((i < is->height() && j < is->width() - 1)
 820                                     ? ppow(get_pixel(2 * i + 1, 2 * j + 1), 0.0625)
 821                                     : pixel(0, 0, 0))));
 822
 823
 824                         for (int k = 0; k < 3; k++)
 825                                 if (!finite(value[k]))
 826                                         value[k] = 0;
 827
 828                         is->set_pixel(i, j, value);
 829                 }
 830
 831                 is->_offset = point(_offset[0] * f, _offset[1] * f);
 832
 833                 return is;
 834         }
 835
 836         /*
 837          * Extend the image area to the top, bottom, left, and right,
 838          * initializing the new image areas with black pixels.  Negative values
 839          * shrink the image.
 840          */
 841         virtual void extend(int top, int bottom, int left, int right) = 0;
 842
 843         /*
 844          * Clone
 845          */
 846         image *clone(char *name) const {
 847                 image *ic = scale_generator(
 848                         height(), width(), depth(), name);
 849
 850                 assert(ic);
 851
 852                 for (unsigned int i = 0; i < height(); i++)
 853                 for (unsigned int j = 0; j < width();  j++)
 854                         ic->set_pixel(i, j,
 855                                 get_pixel(i, j));
 856
 857
 858                 ic->_offset = _offset;
 859
 860                 ic->_apm_memo = _apm_memo;
 861                 ic->_acm_memo = _acm_memo;
 862                 ic->_accm_memo = _accm_memo;
 863                 ic->_apm = _apm;
 864                 ic->_acm = _acm;
 865                 ic->_accm = _accm;
 866
 867                 return ic;
 868         }
 869
 870         /*
 871          * Calculate the average (mean) clamped magnitude of a channel across
 872          * all pixels in an image.  The magnitude is clamped to the range of
 873          * real inputs.
 874          */
 875         ale_real avg_channel_clamped_magnitude(unsigned int k) const {
 876
 877                 /*
 878                  * This is a memoized function
 879                  */
 880
 881                 assert (k < _depth);
 882
 883                 avg_channel_clamped_magnitude_memo();
 884                 return _accm[k];
 885         }
 886
 887         pixel avg_channel_clamped_magnitude() const {
 888                 avg_channel_clamped_magnitude_memo();
 889                 return _accm;
 890         }
 891
 892         /*
 893          * Calculate the average (mean) magnitude of a channel across all
 894          * pixels in an image.
 895          */
 896         ale_real avg_channel_magnitude(unsigned int k) const {
 897
 898                 /*
 899                  * This is a memoized function
 900                  */
 901
 902                 assert (k < _depth);
 903
 904                 avg_channel_magnitude_memo();
 905                 return _acm[k];
 906         }
 907
 908         pixel avg_channel_magnitude() const {
 909                 avg_channel_magnitude_memo();
 910                 return _acm;
 911         }
 912
 913         /*
 914          * Calculate the average (mean) magnitude of a pixel (where magnitude
 915          * is defined as the mean of the channel values).
 916          */
 917         ale_real avg_pixel_magnitude() const {
 918                 unsigned int i, j, k;
 919
 920                 ale_accum accumulator;
 921                 ale_accum divisor = 0;
 922
 923                 if (_apm_memo)
 924                         return _apm;
 925
 926                 _apm_memo = 1;
 927                 accumulator = 0;
 928
 929                 for (i = 0; i < _dimy;  i++)
 930                 for (j = 0; j < _dimx;  j++) {
 931                         pixel value = get_pixel(i, j);
 932
 933                         for (k = 0; k < _depth; k++)
 934                                 if (finite(value[k])) {
 935                                         accumulator += value[k];
 936                                         divisor++;
 937                                 }
 938                 }
 939
 940                 accumulator /= divisor;
 941
 942                 _apm = accumulator;
 943
 944                 return _apm;
 945         }
 946
 947         virtual ~image() {
 948         }
 949 };
 950
 951 #endif