?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
doc/web: Place examples section more visibly.
[Ale.git] / d2 / image.h
blob9ece0be10a105d15752907203f76fbfadcfd2600
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 static double resident;
53 unsigned int _dimx, _dimy, _depth;
54 point _offset;
55 const char *name;
56 mutable exposure *_exp;
57 unsigned int bayer;
58 private:
59 /*
60 * Memoized function variables. We may want to change these even when
61 * *this is constant.
62 */
63 mutable int _apm_memo;
64 mutable ale_real _apm;
65 mutable int _accm_memo;
66 mutable pixel _accm;
67 mutable int _acm_memo;
68 mutable pixel _acm;
69
70 void avg_channel_clamped_magnitude_memo() const {
71 unsigned int i, j, k;
72 pixel_accum accumulator;
73 pixel_accum divisor;
74
75 if (_accm_memo)
76 return;
77
78 _accm_memo = 1;
79
80 accumulator = pixel_accum(0, 0, 0);
81
82 for (i = 0; i < _dimy; i++)
83 for (j = 0; j < _dimx; j++) {
84 pixel value = get_pixel(i, j);
85
86 for (k = 0; k < _depth; k++)
87 if (finite(value[k])) {
88 if (value[k] > 1)
89 value[k] = 1;
90 if (value[k] < 0)
91 value[k] = 0;
92 accumulator[k] += value[k];
93 divisor[k] += 1;
94 }
95 }
96
97 accumulator /= divisor;
98
99 _accm = accumulator;
100 }
101
102 void avg_channel_magnitude_memo() const {
103 unsigned int i, j, k;
104 pixel_accum accumulator;
105 pixel_accum divisor;
106
107 if (_acm_memo)
108 return;
109
110 _acm_memo = 1;
111
112 accumulator = pixel_accum(0, 0, 0);
113
114 for (i = 0; i < _dimy; i++)
115 for (j = 0; j < _dimx; j++) {
116 pixel value = get_pixel(i, j);
117
118 for (k = 0; k < _depth; k++)
119 if (finite(value[k])) {
120 accumulator[k] += value[k];
121 divisor[k] += 1;
122 }
123 }
124
125 accumulator /= divisor;
126
127 _acm = accumulator;
128 }
129
130 protected:
131 void image_updated() {
132 _apm_memo = 0;
133 _acm_memo = 0;
134 _accm_memo = 0;
135 }
136
137 public:
138 static void set_resident(double r) {
139 resident = r;
140 }
141
142 static double get_resident() {
143 return resident;
144 }
145
146 image (unsigned int dimy, unsigned int dimx, unsigned int depth,
147 const char *name = "anonymous", exposure *_exp = NULL,
148 unsigned int bayer = IMAGE_BAYER_NONE) {
149
150 assert (depth == 3);
151 _depth = 3;
152
153 _dimx = dimx;
154 _dimy = dimy;
155 _offset = point(0, 0);
156 _apm_memo = 0;
157 _acm_memo = 0;
158 _accm_memo = 0;
159 this->name = name;
160 this->_exp = _exp;
161 this->bayer = bayer;
162
163 if (_exp != NULL)
164 _exp->add_listener(this, name);
165 }
166
167 unsigned int get_bayer() const {
168 return bayer;
169 }
170
171 virtual char get_channels(int i, int j) const {
172 return 0x7;
173 }
174
175 virtual unsigned int bayer_color(unsigned int i, unsigned int j) const {
176 assert(0);
177 }
178
179 double storage_size() const {
180 if (bayer != IMAGE_BAYER_NONE)
181 return _dimx * _dimy * sizeof(ale_real);
182
183 return 3 * _dimx * _dimy * sizeof(ale_real);
184 }
185
186 exposure &exp() const {
187 return *_exp;
188 }
189
190 point offset() const {
191 return _offset;
192 }
193
194 void set_offset(int i, int j) {
195 _offset[0] = i;
196 _offset[1] = j;
197 }
198
199 void set_offset(point p) {
200 _offset = p;
201 }
202
203 unsigned int width() const {
204 return _dimx;
205 }
206
207 unsigned int height() const {
208 return _dimy;
209 }
210
211 unsigned int depth() const {
212 return _depth;
213 }
214
215 virtual void set_pixel(unsigned int y, unsigned int x, spixel p) = 0;
216
217 virtual spixel get_pixel(unsigned int y, unsigned int x) const = 0;
218
219 virtual spixel get_raw_pixel(unsigned int y, unsigned int x) const {
220 return ((const image *)this)->get_pixel(y, x);
221 }
222
223 virtual void add_pixel(unsigned int y, unsigned int x, pixel p) {
224 assert(0);
225 }
226
227 virtual void mul_pixel(unsigned int y, unsigned int x, pixel p) {
228 assert(0);
229 }
230
231 virtual void div_pixel(unsigned int y, unsigned int x, pixel p) {
232 assert(0);
233 }
234
235 virtual void add_chan(unsigned int y, unsigned int x, unsigned int k, ale_real c) {
236 assert(0);
237 }
238
239 virtual void div_chan(unsigned int y, unsigned int x, unsigned int k, ale_real c) {
240 assert(0);
241 }
242
243 virtual void set_chan(unsigned int y, unsigned int x, unsigned int k, ale_sreal c) = 0;
244
245 virtual ale_sreal get_chan(unsigned int y, unsigned int x, unsigned int k) const = 0;
246
247 ale_real maxval() const {
248 ale_real result = get_pixel(0, 0)[0];
249
250 for (unsigned int i = 0; i < _dimy; i++)
251 for (unsigned int j = 0; j < _dimx; j++) {
252 pixel p = get_pixel(i, j);
253
254 for (unsigned int k = 0; k < _depth; k++)
255 if (p[k] > result || !finite(result))
256 result = p[k];
257 }
258
259 return result;
260 }
261
262 ale_real minval() const {
263 ale_real result = get_pixel(0, 0)[0];
264
265 for (unsigned int i = 0; i < _dimy; i++)
266 for (unsigned int j = 0; j < _dimx; j++) {
267 pixel p = get_pixel(i, j);
268
269 for (unsigned int k = 0; k < _depth; k++)
270 if (p[k] < result || !finite(result))
271 result = p[k];
272 }
273
274 return result;
275 }
276
277 /*
278 * Get the maximum difference among adjacent pixels.
279 */
280
281 pixel get_max_diff(unsigned int i, unsigned int j) const {
282 assert(i <= _dimy - 1);
283 assert(j <= _dimx - 1);
284
285 pixel max = get_pixel(i, j), min = get_pixel(i, j);
286
287 for (int ii = -1; ii <= 1; ii++)
288 for (int jj = -1; jj <= 1; jj++) {
289 int iii = i + ii;
290 int jjj = j + jj;
291
292 if (iii < 0)
293 continue;
294 if (jjj < 0)
295 continue;
296 if ((unsigned int) iii > _dimy - 1)
297 continue;
298 if ((unsigned int) jjj > _dimx - 1)
299 continue;
300
301 pixel p = get_pixel(iii, jjj);
302
303 for (int d = 0; d < 3; d++) {
304 if (p[d] > max[d])
305 max[d] = p[d];
306 if (p[d] < min[d])
307 min[d] = p[d];
308 }
309 }
310
311 return max - min;
312 }
313
314 pixel get_max_diff(point x) const {
315 assert (x[0] >= 0);
316 assert (x[1] >= 0);
317 assert (x[0] <= _dimy - 1);
318 assert (x[1] <= _dimx - 1);
319
320 unsigned int i = (unsigned int) round(x[0]);
321 unsigned int j = (unsigned int) round(x[1]);
322
323 return get_max_diff(i, j);
324 }
325
326 int in_bounds(point x) const {
327 if (x[0] < 0
328 || x[1] < 0
329 || x[0] > height() - 1
330 || x[1] > width() - 1)
331 return 0;
332 if (!x.defined())
333 return 0;
334 return 1;
335 }
336
337 /*
338 * Get a color value at a given position using bilinear interpolation between the
339 * four nearest pixels.
340 */
341 pixel get_bl(point x, int defined = 0) const {
342 // fprintf(stderr, "get_bl x=%f %f\n", (double) x[0], (double) x[1]);
343
344 pixel result;
345
346 assert (x[0] >= 0);
347 assert (x[1] >= 0);
348 assert (x[0] <= _dimy - 1);
349 assert (x[1] <= _dimx - 1);
350
351 int lx = (int) floor(x[1]);
352 int hx = (int) floor(x[1]) + 1;
353 int ly = (int) floor(x[0]);
354 int hy = (int) floor(x[0]) + 1;
355
356 // fprintf(stderr, "get_bl l=%d %d h=%d %d\n", ly, lx, hy, hx);
357
358 pixel neighbor[4];
359 ale_real factor[4];
360
361 neighbor[0] = get_pixel(ly, lx);
362 neighbor[1] = get_pixel(hy % _dimy, lx);
363 neighbor[2] = get_pixel(hy % _dimy, hx % _dimx);
364 neighbor[3] = get_pixel(ly, hx % _dimx);
365
366 // for (int d = 0; d < 4; d++)
367 // fprintf(stderr, "neighbor_%d=%f %f %f\n", d,
368 // (double) neighbor[d][0],
369 // (double) neighbor[d][1],
370 // (double) neighbor[d][2]);
371
372 factor[0] = (ale_real) (hx - x[1]) * (ale_real) (hy - x[0]);
373 factor[1] = (ale_real) (hx - x[1]) * (ale_real) (x[0] - ly);
374 factor[2] = (ale_real) (x[1] - lx) * (ale_real) (x[0] - ly);
375 factor[3] = (ale_real) (x[1] - lx) * (ale_real) (hy - x[0]);
376
377 // for (int d = 0; d < 4; d++)
378 // fprintf(stderr, "factor_%d=%f\n", d,
379 // (double) factor[d]);
380
381 /*
382 * Use bilinear and/or geometric interpolation
383 */
384
385 if (defined == 0) {
386 result = pixel(0, 0, 0);
387
388 for (int n = 0; n < 4; n++)
389 result += factor[n] * neighbor[n];
390 } else {
391 #if 0
392 /*
393 * Calculating the geometric mean may be expensive on
394 * some platforms (e.g., those without floating-point
395 * support.
396 */
397
398 result = pixel(1, 1, 1);
399
400 for (int n = 0; n < 4; n++)
401 result *= ppow(neighbor[n], factor[n]);
402 #else
403 /*
404 * Taking the minimum value may be cheaper than
405 * calculating a geometric mean.
406 */
407
408 result = neighbor[0];
409
410 for (int n = 1; n < 4; n++)
411 for (int k = 0; k < 3; k++)
412 if (neighbor[n][k] < result[k])
413 result[k] = neighbor[n][k];
414 #endif
415 }
416
417 // fprintf(stderr, "result=%f %f %f\n",
418 // (double) result[0],
419 // (double) result[1],
420 // (double) result[2]);
421
422 return result;
423 }
424
425 pixel get_scaled_bl(point x, ale_pos f, int defined = 0) const {
426 point scaled(
427 x[0]/f <= height() - 1
428 ? (x[0]/f)
429 : (ale_pos) (height() - 1),
430 x[1]/f <= width() - 1
431 ? (x[1]/f)
432 : (ale_pos) (width() - 1));
433
434 return get_bl(scaled, defined);
435 }
436
437
438 /*
439 * Make a new image suitable for receiving scaled values.
440 */
441 virtual image *scale_generator(int height, int width, int depth, const char *name) const = 0;
442
443 /*
444 * Generate an image of medians within a given radius
445 */
446
447 image *medians(int radius) const {
448
449 assert (radius >= 0);
450
451 image *is = scale_generator(height(), width(), depth(), "median");
452 assert(is);
453
454 for (unsigned int i = 0; i < height(); i++)
455 for (unsigned int j = 0; j < width(); j++) {
456
457 std::vector<ale_real> p[3];
458
459 for (int ii = -radius; ii <= radius; ii++)
460 for (int jj = -radius; jj <= radius; jj++) {
461 int iii = i + ii;
462 int jjj = j + jj;
463
464 if (in_bounds(point(iii, jjj)))
465 for (int k = 0; k < 3; k++)
466 if (finite(get_pixel(iii, jjj)[k]))
467 p[k].push_back(get_pixel(iii, jjj)[k]);
468 }
469
470 is->set_pixel(i, j, d2::pixel::undefined());
471
472 for (int k = 0; k < 3; k++) {
473 std::sort(p[k].begin(), p[k].end());
474
475 unsigned int pkc = p[k].size();
476
477 if (pkc == 0)
478 continue;
479
480 if (pkc % 2 == 0)
481 is->set_chan(i, j, k,
482 (p[k][pkc / 2] + p[k][pkc / 2 - 1]) / 2);
483 else
484 is->set_chan(i, j, k,
485 p[k][pkc / 2]);
486 }
487 }
488
489 return is;
490 }
491
492 /*
493 * Generate an image of differences of the first channel. The first
494 * coordinate differences are stored in the first channel, second in the
495 * second channel.
496 */
497
498 image *fcdiffs() const {
499 image *is = scale_generator(height(), width(), depth(), "diff");
500
501 assert(is);
502
503 for (unsigned int i = 0; i < height(); i++)
504 for (unsigned int j = 0; j < width(); j++) {
505
506 if (i + 1 < height()
507 && i > 0
508 && !finite(get_chan(i, j, 0))) {
509
510 is->set_chan(i, j, 0, (get_chan(i + 1, j, 0)
511 - get_chan(i - 1, j, 0)) / 2);
512
513 } else if (i + 1 < height()
514 && i > 0
515 && finite(get_chan(i + 1, j, 0))
516 && finite(get_chan(i - 1, j, 0))) {
517
518 is->set_chan(i, j, 0, ((get_chan(i, j, 0) - get_chan(i - 1, j, 0))
519 + (get_chan(i + 1, j, 0) - get_chan(i, j, 0))) / 2);
520
521 } else if (i + 1 < height()
522 && finite(get_chan(i + 1, j, 0))) {
523
524 is->set_chan(i, j, 0, get_chan(i + 1, j, 0) - get_chan(i, j, 0));
525
526 } else if (i > 0
527 && finite(get_chan(i - 1, j, 0))) {
528
529 is->set_chan(i, j, 0, get_chan(i, j, 0) - get_chan(i - 1, j, 0));
530
531 } else {
532 is->set_chan(i, j, 0, 0);
533 }
534
535 if (j + 1 < width()
536 && j > 0
537 && !finite(get_chan(i, j, 0))) {
538
539 is->set_chan(i, j, 1, (get_chan(i, j + 1, 0) - get_chan(i, j - 1, 0)) / 2);
540
541 } else if (j + 1 < width()
542 && j > 0
543 && finite(get_chan(i, j + 1, 0))
544 && finite(get_chan(i, j - 1, 0))) {
545
546 is->set_chan(i, j, 1, ((get_chan(i, j, 0) - get_chan(i, j - 1, 0))
547 + (get_chan(i, j + 1, 0) - get_chan(i, j, 0))) / 2);
548
549 } else if (j + 1 < width() && finite(get_chan(i, j + 1, 0))) {
550
551 is->set_chan(i, j, 1, get_chan(i, j + 1, 0) - get_chan(i, j, 0));
552
553 } else if (j > 0 && finite(get_chan(i, j - 1, 0))) {
554
555 is->set_chan(i, j, 1, get_chan(i, j, 0) - get_chan(i, j - 1, 0));
556
557 } else {
558 is->set_chan(i, j, 1, 0);
559 }
560 }
561
562 return is;
563 }
564
565 /*
566 * Generate an image of median (within a given radius) difference of the
567 * first channel.
568 */
569
570 image *fcdiff_median(int radius) const {
571 image *diff = fcdiffs();
572
573 assert(diff);
574
575 image *median = diff->medians(radius);
576
577 assert(median);
578
579 delete diff;
580
581 return median;
582 }
583
584 /*
585 * Scale by half. We use the following filter:
586 *
587 * 1/16 1/8 1/16
588 * 1/8 1/4 1/8
589 * 1/16 1/8 1/16
590 *
591 * At the edges, these values are normalized so that the sum of the
592 * weights of contributing pixels is 1.
593 */
594 class scale_by_half_threaded : public thread::decompose_domain {
595 image *is;
596 const image *iu;
597 protected:
598 void subdomain_algorithm(unsigned int thread,
599 int i_min, int i_max, int j_min, int j_max) {
600
601 ale_real _0625 = (ale_real) 0.0625;
602 ale_real _125 = (ale_real) 0.125;
603 ale_real _25 = (ale_real) 0.25;
604 ale_real _0 = (ale_real) 0;
605
606 unsigned int ui_min = (unsigned int) i_min;
607 unsigned int ui_max = (unsigned int) i_max;
608 unsigned int uj_min = (unsigned int) j_min;
609 unsigned int uj_max = (unsigned int) j_max;
610
611 for (unsigned int i = ui_min; i < ui_max; i++)
612 for (unsigned int j = uj_min; j < uj_max; j++) {
613 is->set_pixel(i, j,
614
615 ( ( ((i > 0 && j > 0)
616 ? iu->get_pixel(2 * i - 1, 2 * j - 1) * _0625
617 : pixel(0, 0, 0))
618 + ((i > 0)
619 ? iu->get_pixel(2 * i - 1, 2 * j) * _125
620 : pixel(0, 0, 0))
621 + ((i > 0 && j < is->width() - 1)
622 ? iu->get_pixel(2 * i - 1, 2 * j + 1) * _0625
623 : pixel(0, 0, 0))
624 + ((j > 0)
625 ? iu->get_pixel(2 * i, 2 * j - 1) * _125
626 : pixel(0, 0, 0))
627 + iu->get_pixel(2 * i, 2 * j) * _25
628 + ((j < is->width() - 1)
629 ? iu->get_pixel(2 * i, 2 * j + 1) * _125
630 : pixel(0, 0, 0))
631 + ((i < is->height() - 1 && j > 0)
632 ? iu->get_pixel(2 * i + 1, 2 * j - 1) * _0625
633 : pixel(0, 0, 0))
634 + ((i < is->height() - 1)
635 ? iu->get_pixel(2 * i + 1, 2 * j) * _125
636 : pixel(0, 0, 0))
637 + ((i < is->height() && j < is->width() - 1)
638 ? iu->get_pixel(2 * i + 1, 2 * j + 1) * _0625
639 : pixel(0, 0, 0)))
640
641 /
642
643 ( ((i > 0 && j > 0)
644 ? _0625
645 : _0)
646 + ((i > 0)
647 ? _125
648 : _0)
649 + ((i > 0 && j < is->width() - 1)
650 ? _0625
651 : _0)
652 + ((j > 0)
653 ? _125
654 : _0)
655 + _25
656 + ((j < is->width() - 1)
657 ? _125
658 : _0)
659 + ((i < is->height() - 1 && j > 0)
660 ? _0625
661 : _0)
662 + ((i < is->height() - 1)
663 ? _125
664 : _0)
665 + ((i < is->height() && j < is->width() - 1)
666 ? _0625
667 : _0) ) ) );
668 }
669 }
670
671 public:
672 scale_by_half_threaded(image *_is, const image *_iu)
673 : decompose_domain(0, _is->height(),
674 0, _is->width()) {
675 is = _is;
676 iu = _iu;
677 }
678 };
679
680 image *scale_by_half(const char *name) const {
681 ale_pos f = 0.5;
682
683 image *is = scale_generator(
684 (int) floor(height() * (double) f),
685 (int) floor(width() * (double) f), depth(), name);
686
687 assert(is);
688
689 scale_by_half_threaded sbht(is, this);
690 sbht.run();
691
692 is->_offset = point(_offset[0] * f, _offset[1] * f);
693
694 return is;
695 }
696
697 /*
698 * Scale by half. This function uses externally-provided weights,
699 * multiplied by the following filter:
700 *
701 * 1/16 1/8 1/16
702 * 1/8 1/4 1/8
703 * 1/16 1/8 1/16
704 *
705 * Values are normalized so that the sum of the weights of contributing
706 * pixels is 1.
707 */
708 image *scale_by_half(const image *weights, const char *name) const {
709
710 if (weights == NULL)
711 return scale_by_half(name);
712
713 ale_pos f = 0.5;
714
715 image *is = scale_generator(
716 (int) floor(height() * (double) f),
717 (int) floor(width() * (double) f), depth(), name);
718
719 assert(is);
720
721 for (unsigned int i = 0; i < is->height(); i++)
722 for (unsigned int j = 0; j < is->width(); j++) {
723
724 pixel value = pixel
725
726 ( ( ((i > 0 && j > 0)
727 ? (pixel) get_pixel(2 * i - 1, 2 * j - 1)
728 * (pixel) weights->get_pixel(2 * i - 1, 2 * j - 1)
729 * (ale_real) 0.0625
730 : pixel(0, 0, 0))
731 + ((i > 0)
732 ? (pixel) get_pixel(2 * i - 1, 2 * j)
733 * (pixel) weights->get_pixel(2 * i - 1, 2 * j)
734 * 0.125
735 : pixel(0, 0, 0))
736 + ((i > 0 && j < is->width() - 1)
737 ? (pixel) get_pixel(2 * i - 1, 2 * j + 1)
738 * (pixel) weights->get_pixel(2 * i - 1, 2 * j + 1)
739 * 0.0625
740 : pixel(0, 0, 0))
741 + ((j > 0)
742 ? (pixel) get_pixel(2 * i, 2 * j - 1)
743 * (pixel) weights->get_pixel(2 * i, 2 * j - 1)
744 * 0.125
745 : pixel(0, 0, 0))
746 + get_pixel(2 * i, 2 * j)
747 * (pixel) weights->get_pixel(2 * i, 2 * j)
748 * 0.25
749 + ((j < is->width() - 1)
750 ? (pixel) get_pixel(2 * i, 2 * j + 1)
751 * (pixel) weights->get_pixel(2 * i, 2 * j + 1)
752 * 0.125
753 : pixel(0, 0, 0))
754 + ((i < is->height() - 1 && j > 0)
755 ? (pixel) get_pixel(2 * i + 1, 2 * j - 1)
756 * (pixel) weights->get_pixel(2 * i + 1, 2 * j - 1)
757 * 0.0625
758 : pixel(0, 0, 0))
759 + ((i < is->height() - 1)
760 ? (pixel) get_pixel(2 * i + 1, 2 * j)
761 * (pixel) weights->get_pixel(2 * i + 1, 2 * j)
762 * 0.125
763 : pixel(0, 0, 0))
764 + ((i < is->height() && j < is->width() - 1)
765 ? (pixel) get_pixel(2 * i + 1, 2 * j + 1)
766 * (pixel) weights->get_pixel(2 * i + 1, 2 * j + 1)
767 * 0.0625
768 : pixel(0, 0, 0)))
769
770 /
771
772 ( ((i > 0 && j > 0)
773 ? weights->get_pixel(2 * i - 1, 2 * j - 1)
774 * 0.0625
775 : pixel(0, 0, 0))
776 + ((i > 0)
777 ? weights->get_pixel(2 * i - 1, 2 * j)
778 * 0.125
779 : pixel(0, 0, 0))
780 + ((i > 0 && j < is->width() - 1)
781 ? weights->get_pixel(2 * i - 1, 2 * j + 1)
782 * 0.0625
783 : pixel(0, 0, 0))
784 + ((j > 0)
785 ? weights->get_pixel(2 * i, 2 * j - 1)
786 * 0.125
787 : pixel(0, 0, 0))
788 + weights->get_pixel(2 * i, 2 * j)
789 * 0.25
790 + ((j < is->width() - 1)
791 ? weights->get_pixel(2 * i, 2 * j + 1)
792 * 0.125
793 : pixel(0, 0, 0))
794 + ((i < is->height() - 1 && j > 0)
795 ? weights->get_pixel(2 * i + 1, 2 * j - 1)
796 * 0.0625
797 : pixel(0, 0, 0))
798 + ((i < is->height() - 1)
799 ? weights->get_pixel(2 * i + 1, 2 * j)
800 * 0.125
801 : pixel(0, 0, 0))
802 + ((i < is->height() && j < is->width() - 1)
803 ? weights->get_pixel(2 * i + 1, 2 * j + 1)
804 * 0.0625
805 : pixel(0, 0, 0)) ) );
806
807 for (int k = 0; k < 3; k++)
808 if (!finite(value[k]))
809 value[k] = 0;
810
811 is->set_pixel(i, j, value);
812 }
813
814 is->_offset = point(_offset[0] * f, _offset[1] * f);
815
816 return is;
817 }
818
819 /*
820 * Scale an image definition array by 1/2.
821 *
822 * ALE considers an image definition array as a special kind of image
823 * weight array (typedefs of which should appear below the definition
824 * of this class). ALE uses nonzero pixel values to mean 'defined' and
825 * zero values to mean 'undefined'. Through this interpretation, the
826 * image weight array implementation that ALE uses allows image weight
827 * arrays to also serve as image definition arrays.
828 *
829 * Whereas scaling of image weight arrays is not generally obvious in
830 * either purpose or method, ALE requires that image definition arrays
831 * be scalable. (Note that in the special case where weight is treated
832 * as certainty, using a geometric mean is probably correct.)
833 *
834 * We currently use a geometric mean to implement scaling of
835 * definition arrays.
836 */
837
838 class defined_scale_by_half_threaded : public thread::decompose_domain {
839 image *is;
840 const image *iu;
841 protected:
842 void subdomain_algorithm(unsigned int thread,
843 int i_min, int i_max, int j_min, int j_max) {
844
845 #if 0
846 ale_real _0625 = (ale_real) 0.0625;
847 ale_real _125 = (ale_real) 0.125;
848 ale_real _25 = (ale_real) 0.25;
849 #endif
850
851 int ui_min = (int) i_min;
852 int ui_max = (int) i_max;
853 int uj_min = (int) j_min;
854 int uj_max = (int) j_max;
855
856 for (int i = ui_min; i < ui_max; i++)
857 for (int j = uj_min; j < uj_max; j++) {
858
859 #if 0
860
861 /*
862 * Calculate a geometric mean; this approach
863 * may be expensive on some platforms (e.g.,
864 * those without floating-point support in
865 * hardware).
866 */
867
868 pixel value = pixel
869
870 ( ( ((i > 0 && j > 0)
871 ? ppow(iu->get_pixel(2 * i - 1, 2 * j - 1), _0625)
872 : pixel(0, 0, 0))
873 * ((i > 0)
874 ? ppow(iu->get_pixel(2 * i - 1, 2 * j), _125)
875 : pixel(0, 0, 0))
876 * ((i > 0 && j < is->width() - 1)
877 ? ppow(iu->get_pixel(2 * i - 1, 2 * j + 1), _0625)
878 : pixel(0, 0, 0))
879 * ((j > 0)
880 ? ppow(iu->get_pixel(2 * i, 2 * j - 1), _125)
881 : pixel(0, 0, 0))
882 * ppow(iu->get_pixel(2 * i, 2 * j), _25)
883 * ((j < is->width() - 1)
884 ? ppow(iu->get_pixel(2 * i, 2 * j + 1), _125)
885 : pixel(0, 0, 0))
886 * ((i < is->height() - 1 && j > 0)
887 ? ppow(iu->get_pixel(2 * i + 1, 2 * j - 1), _0625)
888 : pixel(0, 0, 0))
889 * ((i < is->height() - 1)
890 ? ppow(iu->get_pixel(2 * i + 1, 2 * j), _125)
891 : pixel(0, 0, 0))
892 * ((i < is->height() && j < is->width() - 1)
893 ? ppow(iu->get_pixel(2 * i + 1, 2 * j + 1), _0625)
894 : pixel(0, 0, 0))));
895 #else
896
897 pixel value = iu->get_pixel(2 * i, 2 * j);
898
899 for (int ii = 2 * i - 1; ii <= 2 * i + 1; ii++)
900 for (int jj = 2 * j - 1; jj <= 2 * j + 1; jj++) {
901 if (ii < 0
902 || jj < 0
903 || ii > (int) iu->height() - 1
904 || jj > (int) iu->height() - 1)
905 continue;
906
907 pixel value2 = iu->get_pixel(ii, jj);
908
909 for (int k = 0; k < 3; k++)
910 if (value2[k] < value[k]
911 || !finite(value2[k])) /* propagate non-finites */
912 value[k] = value2[k];
913 }
914
915 #endif
916
917
918 for (int k = 0; k < 3; k++)
919 if (!finite(value[k]))
920 value[k] = 0;
921
922 is->set_pixel(i, j, value);
923 }
924 }
925
926 public:
927 defined_scale_by_half_threaded(image *_is, const image *_iu)
928 : decompose_domain(0, _is->height(),
929 0, _is->width()) {
930 is = _is;
931 iu = _iu;
932 }
933 };
934
935 image *defined_scale_by_half(const char *name) const {
936 ale_pos f = 0.5;
937
938 image *is = scale_generator(
939 (int) floor(height() * (double) f),
940 (int) floor(width() * (double) f), depth(), name);
941
942 assert(is);
943
944 defined_scale_by_half_threaded dsbht(is, this);
945 dsbht.run();
946
947 is->_offset = point(_offset[0] * f, _offset[1] * f);
948
949 return is;
950 }
951
952 /*
953 * Return an image scaled by some factor != 1.0, using bilinear
954 * interpolation.
955 */
956 image *scale(ale_pos f, const char *name, int defined = 0) const {
957
958 /*
959 * We probably don't want to scale images by a factor of 1.0,
960 * or by non-positive values.
961 */
962 assert (f != 1.0 && f > 0);
963
964 if (f > 1.0) {
965 image *is = scale_generator(
966 (int) floor(height() * (double) f),
967 (int) floor(width() * (double) f), depth(), name);
968
969 assert(is);
970
971 unsigned int i, j, k;
972
973 for (i = 0; i < is->height(); i++)
974 for (j = 0; j < is->width(); j++)
975 for (k = 0; k < is->depth(); k++)
976 is->set_pixel(i, j,
977 get_scaled_bl(point(i, j), f, defined));
978
979 is->_offset = point(_offset[0] * f, _offset[1] * f);
980
981 return is;
982 } else if (f == 0.5) {
983 if (defined == 0)
984 return scale_by_half(name);
985 else
986 return defined_scale_by_half(name);
987 } else {
988 image *is = scale(2*f, name, defined);
989 image *result = is->scale(0.5, name, defined);
990 delete is;
991 return result;
992 }
993
994 }
995
996 /*
997 * Extend the image area to the top, bottom, left, and right,
998 * initializing the new image areas with black pixels. Negative values
999 * shrink the image.
1000 */
1001 virtual image *_extend(int top, int bottom, int left, int right) = 0;
1002
1003 static void extend(image **i, int top, int bottom, int left, int right) {
1004 image *is = (*i)->_extend(top, bottom, left, right);
1005
1006 if (is != NULL) {
1007 delete (*i);
1008 *i = is;
1009 }
1010 }
1011
1012 /*
1013 * Clone
1014 */
1015 image *clone(const char *name) const {
1016 image *ic = scale_generator(
1017 height(), width(), depth(), name);
1018
1019 assert(ic);
1020
1021 for (unsigned int i = 0; i < height(); i++)
1022 for (unsigned int j = 0; j < width(); j++)
1023 ic->set_pixel(i, j,
1024 get_pixel(i, j));
1025
1026
1027 ic->_offset = _offset;
1028
1029 ic->_apm_memo = _apm_memo;
1030 ic->_acm_memo = _acm_memo;
1031 ic->_accm_memo = _accm_memo;
1032 ic->_apm = _apm;
1033 ic->_acm = _acm;
1034 ic->_accm = _accm;
1035
1036 return ic;
1037 }
1038
1039 /*
1040 * Calculate the average (mean) clamped magnitude of a channel across
1041 * all pixels in an image. The magnitude is clamped to the range of
1042 * real inputs.
1043 */
1044 ale_real avg_channel_clamped_magnitude(unsigned int k) const {
1045
1046 /*
1047 * This is a memoized function
1048 */
1049
1050 assert (k < _depth);
1051
1052 avg_channel_clamped_magnitude_memo();
1053 return _accm[k];
1054 }
1055
1056 pixel avg_channel_clamped_magnitude() const {
1057 avg_channel_clamped_magnitude_memo();
1058 return _accm;
1059 }
1060
1061 /*
1062 * Calculate the average (mean) magnitude of a channel across all
1063 * pixels in an image.
1064 */
1065 ale_real avg_channel_magnitude(unsigned int k) const {
1066
1067 /*
1068 * This is a memoized function
1069 */
1070
1071 assert (k < _depth);
1072
1073 avg_channel_magnitude_memo();
1074 return _acm[k];
1075 }
1076
1077 pixel avg_channel_magnitude() const {
1078 avg_channel_magnitude_memo();
1079 return _acm;
1080 }
1081
1082 /*
1083 * Calculate the average (mean) magnitude of a pixel (where magnitude
1084 * is defined as the mean of the channel values).
1085 */
1086 ale_real avg_pixel_magnitude() const {
1087 unsigned int i, j, k;
1088
1089 ale_accum accumulator;
1090 ale_accum divisor = 0;
1091
1092 if (_apm_memo)
1093 return _apm;
1094
1095 _apm_memo = 1;
1096 accumulator = 0;
1097
1098 for (i = 0; i < _dimy; i++)
1099 for (j = 0; j < _dimx; j++) {
1100 pixel value = get_pixel(i, j);
1101
1102 for (k = 0; k < _depth; k++)
1103 if (finite(value[k])) {
1104 accumulator += value[k];
1105 divisor += 1;
1106 }
1107 }
1108
1109 accumulator /= divisor;
1110
1111 _apm = accumulator;
1112
1113 return _apm;
1114 }
1115
1116 virtual ~image() {
1117 }
1118 };
1119
1120 #endif
synthetic capture engine and renderer