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         class scale_by_half_threaded : public thread::decompose_domain {
 529                 image *is;
 530                 const image *iu;
 531         protected:
 532                 void subdomain_algorithm(unsigned int thread,
 533                                 int i_min, int i_max, int j_min, int j_max) {
 534
 535                         unsigned int ui_min = (unsigned int) i_min;
 536                         unsigned int ui_max = (unsigned int) i_max;
 537                         unsigned int uj_min = (unsigned int) j_min;
 538                         unsigned int uj_max = (unsigned int) j_max;
 539
 540                         for (unsigned int i = ui_min; i < ui_max; i++)
 541                         for (unsigned int j = uj_min; j < uj_max; j++) {
 542                                 is->set_pixel(i, j,
 543
 544                                       ( ( ((i > 0 && j > 0)
 545                                             ? iu->get_pixel(2 * i - 1, 2 * j - 1) * (ale_real) 0.0625
 546                                             : pixel(0, 0, 0))
 547                                         + ((i > 0)
 548                                             ? iu->get_pixel(2 * i - 1, 2 * j) * 0.125
 549                                             : pixel(0, 0, 0))
 550                                         + ((i > 0 && j < is->width() - 1)
 551                                             ? iu->get_pixel(2 * i - 1, 2 * j + 1) * 0.0625
 552                                             : pixel(0, 0, 0))
 553                                         + ((j > 0)
 554                                             ? iu->get_pixel(2 * i, 2 * j - 1) * 0.125
 555                                             : pixel(0, 0, 0))
 556                                         + iu->get_pixel(2 * i, 2 * j) * 0.25
 557                                         + ((j < is->width() - 1)
 558                                             ? iu->get_pixel(2 * i, 2 * j + 1) * 0.125
 559                                             : pixel(0, 0, 0))
 560                                         + ((i < is->height() - 1 && j > 0)
 561                                             ? iu->get_pixel(2 * i + 1, 2 * j - 1) * 0.0625
 562                                             : pixel(0, 0, 0))
 563                                         + ((i < is->height() - 1)
 564                                             ? iu->get_pixel(2 * i + 1, 2 * j) * 0.125
 565                                             : pixel(0, 0, 0))
 566                                         + ((i < is->height() && j < is->width() - 1)
 567                                             ? iu->get_pixel(2 * i + 1, 2 * j + 1) * 0.0625
 568                                             : pixel(0, 0, 0)))
 569
 570                                      /
 571
 572                                         ( ((i > 0 && j > 0)
 573                                             ? 0.0625
 574                                             : 0)
 575                                         + ((i > 0)
 576                                             ? 0.125
 577                                             : 0)
 578                                         + ((i > 0 && j < is->width() - 1)
 579                                             ? 0.0625
 580                                             : 0)
 581                                         + ((j > 0)
 582                                             ? 0.125
 583                                             : 0)
 584                                         + 0.25
 585                                         + ((j < is->width() - 1)
 586                                             ? 0.125
 587                                             : 0)
 588                                         + ((i < is->height() - 1 && j > 0)
 589                                             ? 0.0625
 590                                             : 0)
 591                                         + ((i < is->height() - 1)
 592                                             ? 0.125
 593                                             : 0)
 594                                         + ((i < is->height() && j < is->width() - 1)
 595                                             ? 0.0625
 596                                             : 0) ) ) );
 597                         }
 598                 }
 599
 600         public:
 601                 scale_by_half_threaded(image *_is, const image *_iu)
 602                                         : decompose_domain(0, _is->height(),
 603                                                            0, _is->width()) {
 604                         is = _is;
 605                         iu = _iu;
 606                 }
 607         };
 608
 609         image *scale_by_half(char *name) const {
 610                 ale_pos f = 0.5;
 611
 612                 image *is = scale_generator(
 613                         (int) floor(height() * f),
 614                         (int) floor(width()  * f), depth(), name);
 615
 616                 assert(is);
 617
 618                 scale_by_half_threaded sbht(is, this);
 619                 sbht.run();
 620
 621                 is->_offset = point(_offset[0] * f, _offset[1] * f);
 622
 623                 return is;
 624         }
 625
 626         /*
 627          * Scale by half.  This function uses externally-provided weights,
 628          * multiplied by the following filter:
 629          *
 630          *      1/16    1/8     1/16
 631          *      1/8     1/4     1/8
 632          *      1/16    1/8     1/16
 633          *
 634          * Values are normalized so that the sum of the weights of contributing
 635          * pixels is 1.
 636          */
 637         image *scale_by_half(const image *weights, char *name) const {
 638
 639                 if (weights == NULL)
 640                         return scale_by_half(name);
 641
 642                 ale_pos f = 0.5;
 643
 644                 image *is = scale_generator(
 645                         (int) floor(height() * f),
 646                         (int) floor(width()  * f), depth(), name);
 647
 648                 assert(is);
 649
 650                 for (unsigned int i = 0; i < is->height(); i++)
 651                 for (unsigned int j = 0; j < is->width(); j++) {
 652
 653                         pixel value = pixel
 654
 655                               ( ( ((i > 0 && j > 0)
 656                                     ? get_pixel(2 * i - 1, 2 * j - 1)
 657                                     * weights->get_pixel(2 * i - 1, 2 * j - 1)
 658                                     * (ale_real) 0.0625
 659                                     : pixel(0, 0, 0))
 660                                 + ((i > 0)
 661                                     ? get_pixel(2 * i - 1, 2 * j)
 662                                     * weights->get_pixel(2 * i - 1, 2 * j)
 663                                     * 0.125
 664                                     : pixel(0, 0, 0))
 665                                 + ((i > 0 && j < is->width() - 1)
 666                                     ? get_pixel(2 * i - 1, 2 * j + 1)
 667                                     * weights->get_pixel(2 * i - 1, 2 * j + 1)
 668                                     * 0.0625
 669                                     : pixel(0, 0, 0))
 670                                 + ((j > 0)
 671                                     ? get_pixel(2 * i, 2 * j - 1)
 672                                     * weights->get_pixel(2 * i, 2 * j - 1)
 673                                     * 0.125
 674                                     : pixel(0, 0, 0))
 675                                 + get_pixel(2 * i, 2 * j)
 676                                     * weights->get_pixel(2 * i, 2 * j)
 677                                     * 0.25
 678                                 + ((j < is->width() - 1)
 679                                     ? get_pixel(2 * i, 2 * j + 1)
 680                                     * weights->get_pixel(2 * i, 2 * j + 1)
 681                                     * 0.125
 682                                     : pixel(0, 0, 0))
 683                                 + ((i < is->height() - 1 && j > 0)
 684                                     ? get_pixel(2 * i + 1, 2 * j - 1)
 685                                     * weights->get_pixel(2 * i + 1, 2 * j - 1)
 686                                     * 0.0625
 687                                     : pixel(0, 0, 0))
 688                                 + ((i < is->height() - 1)
 689                                     ? get_pixel(2 * i + 1, 2 * j)
 690                                     * weights->get_pixel(2 * i + 1, 2 * j)
 691                                     * 0.125
 692                                     : pixel(0, 0, 0))
 693                                 + ((i < is->height() && j < is->width() - 1)
 694                                     ? get_pixel(2 * i + 1, 2 * j + 1)
 695                                     * weights->get_pixel(2 * i + 1, 2 * j + 1)
 696                                     * 0.0625
 697                                     : pixel(0, 0, 0)))
 698
 699                              /
 700
 701                                 ( ((i > 0 && j > 0)
 702                                     ? weights->get_pixel(2 * i - 1, 2 * j - 1)
 703                                     * 0.0625
 704                                     : pixel(0, 0, 0))
 705                                 + ((i > 0)
 706                                     ? weights->get_pixel(2 * i - 1, 2 * j)
 707                                     * 0.125
 708                                     : pixel(0, 0, 0))
 709                                 + ((i > 0 && j < is->width() - 1)
 710                                     ? weights->get_pixel(2 * i - 1, 2 * j + 1)
 711                                     * 0.0625
 712                                     : pixel(0, 0, 0))
 713                                 + ((j > 0)
 714                                     ? weights->get_pixel(2 * i, 2 * j - 1)
 715                                     * 0.125
 716                                     : pixel(0, 0, 0))
 717                                 + weights->get_pixel(2 * i, 2 * j)
 718                                 * 0.25
 719                                 + ((j < is->width() - 1)
 720                                     ? weights->get_pixel(2 * i, 2 * j + 1)
 721                                     * 0.125
 722                                     : pixel(0, 0, 0))
 723                                 + ((i < is->height() - 1 && j > 0)
 724                                     ? weights->get_pixel(2 * i + 1, 2 * j - 1)
 725                                     * 0.0625
 726                                     : pixel(0, 0, 0))
 727                                 + ((i < is->height() - 1)
 728                                     ? weights->get_pixel(2 * i + 1, 2 * j)
 729                                     * 0.125
 730                                     : pixel(0, 0, 0))
 731                                 + ((i < is->height() && j < is->width() - 1)
 732                                     ? weights->get_pixel(2 * i + 1, 2 * j + 1)
 733                                     * 0.0625
 734                                     : pixel(0, 0, 0)) ) );
 735
 736                         for (int k = 0; k < 3; k++)
 737                                 if (!finite(value[k]))
 738                                         value[k] = 0;
 739
 740                         is->set_pixel(i, j, value);
 741                 }
 742
 743                 is->_offset = point(_offset[0] * f, _offset[1] * f);
 744
 745                 return is;
 746         }
 747
 748         /*
 749          * Scale an image definition array by 1/2.
 750          *
 751          * ALE considers an image definition array as a special kind of image
 752          * weight array (typedefs of which should appear below the definition
 753          * of this class).  ALE uses nonzero pixel values to mean 'defined' and
 754          * zero values to mean 'undefined'.  Through this interpretation, the
 755          * image weight array implementation that ALE uses allows image weight
 756          * arrays to also serve as image definition arrays.
 757          *
 758          * Whereas scaling of image weight arrays is not generally obvious in
 759          * either purpose or method, ALE requires that image definition arrays
 760          * be scalable.  (Note that in the special case where weight is treated
 761          * as certainty, using a geometric mean is probably correct.)
 762          *
 763          * We currently use a geometric mean to implement scaling of
 764          * definition arrays.
 765          */
 766
 767         class defined_scale_by_half_threaded : public thread::decompose_domain {
 768                 image *is;
 769                 const image *iu;
 770         protected:
 771                 void subdomain_algorithm(unsigned int thread,
 772                                 int i_min, int i_max, int j_min, int j_max) {
 773
 774                         unsigned int ui_min = (unsigned int) i_min;
 775                         unsigned int ui_max = (unsigned int) i_max;
 776                         unsigned int uj_min = (unsigned int) j_min;
 777                         unsigned int uj_max = (unsigned int) j_max;
 778
 779                         for (unsigned int i = ui_min; i < ui_max; i++)
 780                         for (unsigned int j = uj_min; j < uj_max; j++) {
 781
 782                                 pixel value = pixel
 783
 784                                       ( ( ((i > 0 && j > 0)
 785                                             ? ppow(iu->get_pixel(2 * i - 1, 2 * j - 1), 0.0625)
 786                                             : pixel(0, 0, 0))
 787                                         * ((i > 0)
 788                                             ? ppow(iu->get_pixel(2 * i - 1, 2 * j), 0.125)
 789                                             : pixel(0, 0, 0))
 790                                         * ((i > 0 && j < is->width() - 1)
 791                                             ? ppow(iu->get_pixel(2 * i - 1, 2 * j + 1), 0.0625)
 792                                             : pixel(0, 0, 0))
 793                                         * ((j > 0)
 794                                             ? ppow(iu->get_pixel(2 * i, 2 * j - 1), 0.125)
 795                                             : pixel(0, 0, 0))
 796                                         * ppow(iu->get_pixel(2 * i, 2 * j), 0.25)
 797                                         * ((j < is->width() - 1)
 798                                             ? ppow(iu->get_pixel(2 * i, 2 * j + 1), 0.125)
 799                                             : pixel(0, 0, 0))
 800                                         * ((i < is->height() - 1 && j > 0)
 801                                             ? ppow(iu->get_pixel(2 * i + 1, 2 * j - 1), 0.0625)
 802                                             : pixel(0, 0, 0))
 803                                         * ((i < is->height() - 1)
 804                                             ? ppow(iu->get_pixel(2 * i + 1, 2 * j), 0.125)
 805                                             : pixel(0, 0, 0))
 806                                         * ((i < is->height() && j < is->width() - 1)
 807                                             ? ppow(iu->get_pixel(2 * i + 1, 2 * j + 1), 0.0625)
 808                                             : pixel(0, 0, 0))));
 809
 810
 811                                 for (int k = 0; k < 3; k++)
 812                                         if (!finite(value[k]))
 813                                                 value[k] = 0;
 814
 815                                 is->set_pixel(i, j, value);
 816                         }
 817                 }
 818
 819         public:
 820                 defined_scale_by_half_threaded(image *_is, const image *_iu)
 821                                         : decompose_domain(0, _is->height(),
 822                                                            0, _is->width()) {
 823                         is = _is;
 824                         iu = _iu;
 825                 }
 826         };
 827
 828         image *defined_scale_by_half(char *name) const {
 829                 ale_pos f = 0.5;
 830
 831                 image *is = scale_generator(
 832                         (int) floor(height() * f),
 833                         (int) floor(width()  * f), depth(), name);
 834
 835                 assert(is);
 836
 837                 defined_scale_by_half_threaded dsbht(is, this);
 838                 dsbht.run();
 839
 840                 is->_offset = point(_offset[0] * f, _offset[1] * f);
 841
 842                 return is;
 843         }
 844
 845         /*
 846          * Return an image scaled by some factor != 1.0, using bilinear
 847          * interpolation.
 848          */
 849         image *scale(ale_pos f, char *name, int defined = 0) const {
 850
 851                 /*
 852                  * We probably don't want to scale images by a factor of 1.0,
 853                  * or by non-positive values.
 854                  */
 855                 assert (f != 1.0 && f > 0);
 856
 857                 if (f > 1.0) {
 858                         image *is = scale_generator(
 859                                 (int) floor(height() * f),
 860                                 (int) floor(width()  * f), depth(), name);
 861
 862                         assert(is);
 863
 864                         unsigned int i, j, k;
 865
 866                         for (i = 0; i < is->height(); i++)
 867                         for (j = 0; j < is->width(); j++)
 868                         for (k = 0; k < is->depth(); k++)
 869                                 is->set_pixel(i, j,
 870                                         get_scaled_bl(point(i, j), f, defined));
 871
 872                         is->_offset = point(_offset[0] * f, _offset[1] * f);
 873
 874                         return is;
 875                 } else if (f == 0.5) {
 876                         if (defined == 0)
 877                                 return scale_by_half(name);
 878                         else
 879                                 return defined_scale_by_half(name);
 880                 } else {
 881                         image *is = scale(2*f, name, defined);
 882                         image *result = is->scale(0.5, name, defined);
 883                         delete is;
 884                         return result;
 885                 }
 886
 887         }
 888
 889         /*
 890          * Extend the image area to the top, bottom, left, and right,
 891          * initializing the new image areas with black pixels.  Negative values
 892          * shrink the image.
 893          */
 894         virtual void extend(int top, int bottom, int left, int right) = 0;
 895
 896         /*
 897          * Clone
 898          */
 899         image *clone(char *name) const {
 900                 image *ic = scale_generator(
 901                         height(), width(), depth(), name);
 902
 903                 assert(ic);
 904
 905                 for (unsigned int i = 0; i < height(); i++)
 906                 for (unsigned int j = 0; j < width();  j++)
 907                         ic->set_pixel(i, j,
 908                                 get_pixel(i, j));
 909
 910
 911                 ic->_offset = _offset;
 912
 913                 ic->_apm_memo = _apm_memo;
 914                 ic->_acm_memo = _acm_memo;
 915                 ic->_accm_memo = _accm_memo;
 916                 ic->_apm = _apm;
 917                 ic->_acm = _acm;
 918                 ic->_accm = _accm;
 919
 920                 return ic;
 921         }
 922
 923         /*
 924          * Calculate the average (mean) clamped magnitude of a channel across
 925          * all pixels in an image.  The magnitude is clamped to the range of
 926          * real inputs.
 927          */
 928         ale_real avg_channel_clamped_magnitude(unsigned int k) const {
 929
 930                 /*
 931                  * This is a memoized function
 932                  */
 933
 934                 assert (k < _depth);
 935
 936                 avg_channel_clamped_magnitude_memo();
 937                 return _accm[k];
 938         }
 939
 940         pixel avg_channel_clamped_magnitude() const {
 941                 avg_channel_clamped_magnitude_memo();
 942                 return _accm;
 943         }
 944
 945         /*
 946          * Calculate the average (mean) magnitude of a channel across all
 947          * pixels in an image.
 948          */
 949         ale_real avg_channel_magnitude(unsigned int k) const {
 950
 951                 /*
 952                  * This is a memoized function
 953                  */
 954
 955                 assert (k < _depth);
 956
 957                 avg_channel_magnitude_memo();
 958                 return _acm[k];
 959         }
 960
 961         pixel avg_channel_magnitude() const {
 962                 avg_channel_magnitude_memo();
 963                 return _acm;
 964         }
 965
 966         /*
 967          * Calculate the average (mean) magnitude of a pixel (where magnitude
 968          * is defined as the mean of the channel values).
 969          */
 970         ale_real avg_pixel_magnitude() const {
 971                 unsigned int i, j, k;
 972
 973                 ale_accum accumulator;
 974                 ale_accum divisor = 0;
 975
 976                 if (_apm_memo)
 977                         return _apm;
 978
 979                 _apm_memo = 1;
 980                 accumulator = 0;
 981
 982                 for (i = 0; i < _dimy;  i++)
 983                 for (j = 0; j < _dimx;  j++) {
 984                         pixel value = get_pixel(i, j);
 985
 986                         for (k = 0; k < _depth; k++)
 987                                 if (finite(value[k])) {
 988                                         accumulator += value[k];
 989                                         divisor++;
 990                                 }
 991                 }
 992
 993                 accumulator /= divisor;
 994
 995                 _apm = accumulator;
 996
 997                 return _apm;
 998         }
 999
1000         virtual ~image() {
1001         }
1002 };
1003
1004 #endif