MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
image_tools.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#ifndef MVE_IMAGE_TOOLS_HEADER
11#define MVE_IMAGE_TOOLS_HEADER
12
13#include <iostream>
14#include <limits>
15#include <complex>
16#include <type_traits>
17
18#include "util/exception.h"
19#include "math/accum.h"
20#include "math/functions.h"
21#include "mve/defines.h"
22#include "mve/image.h"
23
26
27/*
28 * ----------------------- Image conversions ---------------------
29 */
30
35FloatImage::Ptr
36byte_to_float_image (ByteImage::ConstPtr image);
37
41DoubleImage::Ptr
42byte_to_double_image (ByteImage::ConstPtr image);
43
48ByteImage::Ptr
49float_to_byte_image (FloatImage::ConstPtr image,
50 float vmin = 0.0f, float vmax = 1.0f);
51
56ByteImage::Ptr
57double_to_byte_image (DoubleImage::ConstPtr image,
58 double vmin = 0.0, double vmax = 1.0);
59
64ByteImage::Ptr
65int_to_byte_image (IntImage::ConstPtr image);
66
71ByteImage::Ptr
72raw_to_byte_image (RawImage::ConstPtr image,
73 uint16_t vmin = 0, uint16_t vmax = 65535);
74
79FloatImage::Ptr
80raw_to_float_image(RawImage::ConstPtr image);
81
86template <typename SRC, typename DST>
87typename Image<DST>::Ptr
88type_to_type_image (typename Image<SRC>::ConstPtr image);
89
93template <typename T>
94void
95find_min_max_value (typename mve::Image<T>::ConstPtr image, T* vmin, T* vmax);
96
102void
103float_image_normalize (FloatImage::Ptr image);
104
105/*
106 * ----------------------- Image undistortion ---------------------
107 */
108
115template <typename T>
116typename Image<T>::Ptr
117image_undistort_msps (typename Image<T>::ConstPtr img, double k0, double k1);
118
126template <typename T>
127typename Image<T>::Ptr
128image_undistort_k2k4 (typename Image<T>::ConstPtr img,
129 double focal_length, double k2, double k4);
130
137template <typename T>
138typename Image<T>::Ptr
139image_undistort_vsfm (typename Image<T>::ConstPtr img,
140 double focal_length, double k1);
141
142/*
143 * ------------------- Image scaling and cropping -----------------
144 */
145
151template <typename T>
152typename Image<T>::Ptr
153crop (typename Image<T>::ConstPtr image, int64_t width, int64_t height,
154 int64_t left, int64_t top, T const* fill_color);
155
158{
162 //RESCALE_CUBIC // Not implemented
164
171template <typename T>
172typename Image<T>::Ptr
173rescale (typename Image<T>::ConstPtr image, RescaleInterpolation interp,
174 int64_t width, int64_t height);
175
181template <typename T>
182typename Image<T>::Ptr
183rescale_half_size (typename Image<T>::ConstPtr image);
184
191template <typename T>
192typename Image<T>::Ptr
193rescale_half_size_gaussian (typename Image<T>::ConstPtr image,
194 float sigma = 0.866025403784439f);
195
200template <typename T>
201typename Image<T>::Ptr
202rescale_half_size_subsample (typename Image<T>::ConstPtr image);
203
210template <typename T>
211typename Image<T>::Ptr
212rescale_double_size (typename Image<T>::ConstPtr img);
213
221template <typename T>
222typename Image<T>::Ptr
223rescale_double_size_supersample (typename Image<T>::ConstPtr img);
224
231template <typename T>
232void
233rescale_nearest (typename Image<T>::ConstPtr in, typename Image<T>::Ptr out);
234
241template <typename T>
242void
243rescale_linear (typename Image<T>::ConstPtr in, typename Image<T>::Ptr out);
244
252template <typename T>
253void
254rescale_gaussian (typename Image<T>::ConstPtr in,
255 typename Image<T>::Ptr out, float sigma_factor = 1.0f);
256
257/*
258 * ------------------------- Image blurring --------------------------
259 */
260
265template <typename T>
266typename Image<T>::Ptr
267blur_gaussian (typename Image<T>::ConstPtr in, float sigma);
268
274template <typename T>
275typename Image<T>::Ptr
276blur_boxfilter (typename Image<T>::ConstPtr in, int ks);
277
278/*
279 * ------------------- Image rotation and flipping -------------------
280 */
281
290
299template <typename T>
300typename Image<T>::Ptr
301rotate (typename Image<T>::ConstPtr image, RotateType type);
302
309template <typename T>
310typename Image<T>::Ptr
311rotate (typename Image<T>::ConstPtr image, float angle, T const* fill_color);
312
321
325template <typename T>
326void
327flip (typename Image<T>::Ptr image, FlipType type);
328
329/*
330 * ----------------------- Image desaturation ------------------------
331 */
332
342
355template <typename T>
356typename Image<T>::Ptr
357desaturate (typename Image<T>::ConstPtr image, DesaturateType type);
358
362template <typename T>
363typename Image<T>::Ptr
364expand_grayscale (typename Image<T>::ConstPtr image);
365
367template <typename T>
368void
369reduce_alpha (typename mve::Image<T>::Ptr img); // TODO Blend color?
370
371/* ------------------------- Edge detection ----------------------- */
372
379template <typename T>
380typename mve::Image<T>::Ptr
381sobel_edge (typename mve::Image<T>::ConstPtr img);
382
383/* ------------------------- Miscellaneous ------------------------ */
384
389template <typename T>
390typename Image<T>::Ptr
391subtract (typename Image<T>::ConstPtr i1, typename Image<T>::ConstPtr i2);
392
398template <typename T>
399typename Image<T>::Ptr
400difference (typename Image<T>::ConstPtr i1, typename Image<T>::ConstPtr i2);
401
407template <typename T>
408void
409gamma_correct (typename Image<T>::Ptr image, T const& power);
410
415void
416gamma_correct (ByteImage::Ptr image, float power);
417
428template <typename T>
429void
430gamma_correct_srgb (typename Image<T>::Ptr image);
431
442template <typename T>
443void
444gamma_correct_inv_srgb (typename Image<T>::Ptr image);
445
451template <typename T_IN, typename T_OUT>
452typename Image<T_OUT>::Ptr
453integral_image (typename Image<T_IN>::ConstPtr image);
454
461template <typename T>
462T
463integral_image_area (typename Image<T>::ConstPtr sat,
464 int64_t x1, int64_t y1, int64_t x2, int64_t y2, int64_t cc = 0);
465
470template <typename T>
471typename Image<T>::Ptr
472create_thumbnail (typename Image<T>::ConstPtr image,
473 int64_t thumb_width, int64_t thumb_height);
474
477
478/* ------------------------- Implementation ----------------------- */
479
480#include <cmath>
481#include <stdexcept>
482
485
486template <typename SRC, typename DST>
487typename Image<DST>::Ptr
489{
490 typename Image<DST>::Ptr out = Image<DST>::create();
491 out->allocate(image->width(), image->height(), image->channels());
492 int64_t size = image->get_value_amount();
493 SRC const* src_buf = image->get_data_pointer();
494 DST* dst_buf = out->get_data_pointer();
495 std::copy_n(src_buf, size, dst_buf);
496 return out;
497}
498
499/* ---------------------------------------------------------------- */
500
501template <typename T>
502void
503find_min_max_value (typename mve::Image<T>::ConstPtr image, T* vmin, T* vmax)
504{
505 *vmin = std::numeric_limits<T>::max();
506 *vmax = std::numeric_limits<T>::is_signed
507 ? -std::numeric_limits<T>::max()
508 : std::numeric_limits<T>::min();
509
510 for (T const* ptr = image->begin(); ptr != image->end(); ++ptr)
511 {
512 if (*ptr < *vmin)
513 *vmin = *ptr;
514 if (*ptr > *vmax)
515 *vmax = *ptr;
516 }
517}
518
519/* ---------------------------------------------------------------- */
520
521template <typename T>
522typename Image<T>::Ptr
524 int64_t width, int64_t height)
525{
526 if (width == 0 && height == 0)
527 throw std::invalid_argument("Invalid size request");
528
529 /* Duplicate input image if width and height match image size. */
530 if (width == img->width() && height == img->height())
531 return Image<T>::create(*img);
532
533 /* Keep aspect ratio if one of width or height is 0. */
534 if (width == 0)
535 width = height * img->width() / img->height();
536 else if (height == 0)
537 height = width * img->height() / img->width();
538
539 /* Scale down to an appropriate mipmap level for resizing. */
540 if (interp == RESCALE_NEAREST || interp == RESCALE_LINEAR)
541 {
542 while (2 * width <= img->width() && 2 * height <= img->height())
543 img = rescale_half_size<T>(img);
544 }
545
546 typename Image<T>::Ptr out(Image<T>::create());
547 out->allocate(width, height, img->channels());
548
549 switch (interp)
550 {
551 case RESCALE_NEAREST:
552 rescale_nearest<T>(img, out);
553 break;
554 case RESCALE_LINEAR:
555 rescale_linear<T>(img, out);
556 break;
557#if 0
558 case RESCALE_CUBIC:
559 rescale_cubic(img, out);
560 break;
561#endif
562 case RESCALE_GAUSSIAN:
563 rescale_gaussian<T>(img, out);
564 break;
565
566 default:
567 throw std::invalid_argument("Invalid interpolation type");
568 }
569
570 return out;
571}
572
573/* ---------------------------------------------------------------- */
574
575template <typename T>
576typename Image<T>::Ptr
578{
579 int64_t const iw = img->width();
580 int64_t const ih = img->height();
581 int64_t const ic = img->channels();
582 int64_t const ow = (iw + 1) >> 1;
583 int64_t const oh = (ih + 1) >> 1;
584
585 if (iw < 2 || ih < 2)
586 throw std::invalid_argument("Input image too small for half-sizing");
587
588 typename Image<T>::Ptr out(Image<T>::create());
589 out->allocate(ow, oh, ic);
590
591 int64_t outpos = 0;
592 int64_t rowstride = iw * ic;
593 for (int64_t y = 0; y < oh; ++y)
594 {
595 int64_t irow1 = y * 2 * rowstride;
596 int64_t irow2 = irow1 + rowstride * (y * 2 + 1 < ih);
597
598 for (int64_t x = 0; x < ow; ++x)
599 {
600 int64_t ipix1 = irow1 + x * 2 * ic;
601 int64_t ipix2 = irow2 + x * 2 * ic;
602 int64_t hasnext = (x * 2 + 1 < iw);
603
604 for (int64_t c = 0; c < ic; ++c)
605 out->at(outpos++) = math::interpolate<T>(
606 img->at(ipix1 + c), img->at(ipix1 + ic * hasnext + c),
607 img->at(ipix2 + c), img->at(ipix2 + ic * hasnext + c),
608 0.25f, 0.25f, 0.25f, 0.25f);
609 }
610 }
611
612 return out;
613}
614
615/* ---------------------------------------------------------------- */
616
617template <typename T>
618typename Image<T>::Ptr
620{
621 int64_t const iw = img->width();
622 int64_t const ih = img->height();
623 int64_t const ic = img->channels();
624 int64_t const ow = (iw + 1) >> 1;
625 int64_t const oh = (ih + 1) >> 1;
626
627 if (iw < 2 || ih < 2)
628 throw std::invalid_argument("Invalid input image");
629
630 typename Image<T>::Ptr out(Image<T>::create());
631 out->allocate(ow, oh, ic);
632
633 /*
634 * Weights w1 (4 center px), w2 (8 skewed px) and w3 (4 corner px).
635 * Weights can be normalized by dividing with (4*w1 + 8*w2 + 4*w3).
636 * Since the accumulator is used, normalization is not explicit.
637 */
638 float const w1 = std::exp(-0.5f / (2.0f * MATH_POW2(sigma)));
639 float const w2 = std::exp(-2.5f / (2.0f * MATH_POW2(sigma)));
640 float const w3 = std::exp(-4.5f / (2.0f * MATH_POW2(sigma)));
641
642 int64_t outpos = 0;
643 int64_t const rowstride = iw * ic;
644 for (int64_t y = 0; y < oh; ++y)
645 {
646 /* Init the four row pointers. */
647 int64_t const y2 = y * 2;
648 T const* row[4];
649 row[0] = &img->at(std::max<int64_t>(0, y2 - 1) * rowstride);
650 row[1] = &img->at(y2 * rowstride);
651 row[2] = &img->at(std::min(ih - 1, y2 + 1) * rowstride);
652 row[3] = &img->at(std::min(ih - 1, y2 + 2) * rowstride);
653
654 for (int64_t x = 0; x < ow; ++x)
655 {
656 /* Init four pixel positions for each row. */
657 int64_t x2 = x * 2;
658 int64_t xi[4];
659 xi[0] = std::max<int64_t>(0, x2 - 1) * ic;
660 xi[1] = x2 * ic;
661 xi[2] = std::min(iw - 1, x2 + 1) * ic;
662 xi[3] = std::min(iw - 1, x2 + 2) * ic;
663
664 /* Accumulate 16 values in each channel. */
665 for (int64_t c = 0; c < ic; ++c)
666 {
667 math::Accum<T> accum(T(0));
668 accum.add(row[0][xi[0] + c], w3);
669 accum.add(row[0][xi[1] + c], w2);
670 accum.add(row[0][xi[2] + c], w2);
671 accum.add(row[0][xi[3] + c], w3);
672
673 accum.add(row[1][xi[0] + c], w2);
674 accum.add(row[1][xi[1] + c], w1);
675 accum.add(row[1][xi[2] + c], w1);
676 accum.add(row[1][xi[3] + c], w2);
677
678 accum.add(row[2][xi[0] + c], w2);
679 accum.add(row[2][xi[1] + c], w1);
680 accum.add(row[2][xi[2] + c], w1);
681 accum.add(row[2][xi[3] + c], w2);
682
683 accum.add(row[3][xi[0] + c], w3);
684 accum.add(row[3][xi[1] + c], w2);
685 accum.add(row[3][xi[2] + c], w2);
686 accum.add(row[3][xi[3] + c], w3);
687
688 out->at(outpos++) = accum.normalized();
689 }
690 }
691 }
692
693 return out;
694}
695
696/* ---------------------------------------------------------------- */
697
698template <typename T>
699typename Image<T>::Ptr
701{
702 int64_t const iw = img->width();
703 int64_t const ih = img->height();
704 int64_t const ic = img->channels();
705 int64_t const ow = (iw + 1) >> 1;
706 int64_t const oh = (ih + 1) >> 1;
707 int64_t const irs = iw * ic; // input image row stride
708
709 typename Image<T>::Ptr out(Image<T>::create());
710 out->allocate(ow, oh, ic);
711
712 int64_t iter = 0; // Output image iterator
713 for (int64_t iy = 0; iy < ih; iy += 2)
714 {
715 int64_t rowoff = iy * irs;
716 int64_t pixoff = rowoff;
717 for (int64_t ix = 0; ix < iw; ix += 2)
718 {
719 T const* iptr = &img->at(pixoff);
720 T* optr = &out->at(iter);
721 std::copy(iptr, iptr + ic, optr);
722 pixoff += ic << 1;
723 iter += ic;
724 }
725 }
726
727 return out;
728}
729
730/* ---------------------------------------------------------------- */
731
732template <typename T>
733typename Image<T>::Ptr
735{
736 int64_t const iw = img->width();
737 int64_t const ih = img->height();
738 int64_t const ic = img->channels();
739 int64_t const ow = iw << 1;
740 int64_t const oh = ih << 1;
741 int64_t const irs = iw * ic; // input image row stride
742
743 typename Image<T>::Ptr out(Image<T>::create());
744 out->allocate(ow, oh, ic);
745
746 float w[4] = { 0.75f*0.75f, 0.25f*0.75f, 0.75f*0.25f, 0.25f*0.25f };
747
748 T const* row1 = img->begin();
749 T const* row2 = img->begin();
750
751 int64_t i = 0;
752 for (int64_t y = 0; y < oh; ++y)
753 {
754 /* Uneven row -> advance, even row -> swap. */
755 if (y % 2)
756 row2 = row1 + (y < oh - 1 ? irs : 0);
757 else
758 std::swap(row1, row2);
759
760 T const* px[4] = { row1, row1, row2, row2 };
761 for (int64_t x = 0; x < ow; ++x)
762 {
763 /* Uneven pixel -> advance, even pixel -> swap. */
764 if (x % 2)
765 {
766 int64_t off = (x < ow - 1 ? ic : 0);
767 px[1] = px[0] + off;
768 px[3] = px[2] + off;
769 }
770 else
771 {
772 std::swap(px[0], px[1]);
773 std::swap(px[2], px[3]);
774 }
775
776 for (int64_t c = 0; c < ic; ++c, ++i)
777 out->at(i) = math::interpolate
778 (px[0][c], px[1][c], px[2][c], px[3][c],
779 w[0], w[1], w[2], w[3]);
780 }
781 }
782
783 return out;
784}
785
786/* ---------------------------------------------------------------- */
787
788template <typename T>
789typename Image<T>::Ptr
791{
792 int64_t const iw = img->width();
793 int64_t const ih = img->height();
794 int64_t const ic = img->channels();
795 int64_t const ow = iw << 1;
796 int64_t const oh = ih << 1;
797
798 typename Image<T>::Ptr out(Image<T>::create());
799 out->allocate(ow, oh, ic);
800
801 int64_t witer = 0;
802 for (int64_t y = 0; y < oh; ++y)
803 {
804 bool nexty = (y + 1 < oh);
805 int64_t yoff[2] = { iw * (y >> 1), iw * ((y + nexty) >> 1) };
806 for (int64_t x = 0; x < ow; ++x)
807 {
808 bool nextx = (x + 1 < ow);
809 int64_t xoff[2] = { x >> 1, (x + nextx) >> 1 };
810 T const* val[4] =
811 {
812 &img->at(yoff[0] + xoff[0], 0),
813 &img->at(yoff[0] + xoff[1], 0),
814 &img->at(yoff[1] + xoff[0], 0),
815 &img->at(yoff[1] + xoff[1], 0)
816 };
817
818 for (int64_t c = 0; c < ic; ++c, ++witer)
819 out->at(witer) = math::interpolate
820 (val[0][c], val[1][c], val[2][c], val[3][c],
821 0.25f, 0.25f, 0.25f, 0.25f);
822 }
823 }
824
825 return out;
826}
827
828/* ---------------------------------------------------------------- */
829
830template <typename T>
831void
833{
834 if (img->channels() != out->channels())
835 throw std::invalid_argument("Image channel mismatch");
836
837 int64_t const iw = img->width();
838 int64_t const ih = img->height();
839 int64_t const ic = img->channels();
840 int64_t const ow = out->width();
841 int64_t const oh = out->height();
842
843 int64_t outpos = 0;
844 for (int64_t y = 0; y < oh; ++y)
845 {
846 float ly = ((float)y + 0.5f) * (float)ih / (float)oh;
847 int64_t iy = static_cast<int64_t>(ly);
848 for (int64_t x = 0; x < ow; ++x)
849 {
850 float lx = ((float)x + 0.5f) * (float)iw / (float)ow;
851 int64_t ix = static_cast<int64_t>(lx);
852 for (int64_t c = 0; c < ic; ++c)
853 out->at(outpos++) = img->at(ix,iy,c);
854 }
855 }
856}
857
858/* ---------------------------------------------------------------- */
859
860template <typename T>
861void
863{
864 if (img->channels() != out->channels())
865 throw std::invalid_argument("Image channel mismatch");
866
867 int64_t const iw = img->width();
868 int64_t const ih = img->height();
869 int64_t const ic = img->channels();
870 int64_t const ow = out->width();
871 int64_t const oh = out->height();
872
873 T* out_ptr = out->get_data_pointer();
874 int64_t outpos = 0;
875 for (int64_t y = 0; y < oh; ++y)
876 {
877 float fy = ((float)y + 0.5f) * (float)ih / (float)oh;
878 for (int64_t x = 0; x < ow; ++x, outpos += ic)
879 {
880 float fx = ((float)x + 0.5f) * (float)iw / (float)ow;
881 img->linear_at(fx - 0.5f, fy - 0.5f, out_ptr + outpos);
882 }
883 }
884}
885
886/* ---------------------------------------------------------------- */
887
888template <typename T>
889T
891 float x, float y, int64_t c, float sigma)
892{
893 int64_t const width = img->width();
894 int64_t const height = img->height();
895
896 /*
897 * Calculate kernel size for geometric gaussian
898 * Kernel is cut off at y=1/N, x = sigma * sqrt(2 * ln N).
899 *
900 * For N=256: x = sigma * 3.33.
901 * For N=128: x = sigma * 3.12.
902 * For N=64: x = sigma * 2.884.
903 * For N=32: x = sigma * 2.63.
904 * For N=16: x = sigma * 2.355.
905 * For N=8: x = sigma * 2.04.
906 * For N=4: x = sigma * 1.67.
907 */
908 float ks = sigma * 2.884f;
909
910 /* Calculate min/max kernel position. */
911 float kx_min = std::floor(x - ks);
912 float kx_max = std::ceil(x + ks - 1.0f);
913 float ky_min = std::floor(y - ks);
914 float ky_max = std::ceil(y + ks - 1.0f);
915
916 int64_t kxi_min = static_cast<int64_t>(std::max(0.0f, kx_min));
917 int64_t kxi_max =
918 static_cast<int64_t>(std::min((float)width - 1.0f, kx_max));
919 int64_t kyi_min = static_cast<int64_t>(std::max(0.0f, ky_min));
920 int64_t kyi_max =
921 static_cast<int64_t>(std::min((float)height - 1.0f, ky_max));
922
923 /* Determine pixel weight for kernel bounaries. */
924 float wx_start = kx_min > 0.0f ? kx_min + 1.0f + ks - x : 1.0f;
925 float wx_end = kx_max < (float)width - 1.0f ? ks + x - kx_max : 1.0f;
926 float wy_start = ky_min > 0.0f ? ky_min + 1.0f + ks - y : 1.0f;
927 float wy_end = ky_max < (float)height - 1.0f ? ks + y - ky_max : 1.0f;
928
929 /* Apply kernel. */
930 math::Accum<T> accum(0);
931 for (int64_t yi = kyi_min; yi <= kyi_max; ++yi)
932 for (int64_t xi = kxi_min; xi <= kxi_max; ++xi)
933 {
934 float weight = 1.0f;
935 weight *= (xi == kxi_min ? wx_start : 1.0f);
936 weight *= (xi == kxi_max ? wx_end : 1.0f);
937 weight *= (yi == kyi_min ? wy_start : 1.0f);
938 weight *= (yi == kyi_max ? wy_end : 1.0f);
939 float dx = static_cast<float>(xi) + 0.5f - x;
940 float dy = static_cast<float>(yi) + 0.5f - y;
941 weight *= math::gaussian_xx(dx * dx + dy * dy, sigma);
942
943 accum.add(img->at(xi, yi, c), weight);
944 }
945
946 return accum.normalized();
947}
948
949/* ---------------------------------------------------------------- */
950
951template <typename T>
952void
954 typename Image<T>::Ptr out, float sigma_factor)
955{
956 if (img->channels() != out->channels())
957 throw std::invalid_argument("Image channels mismatch");
958
959 int64_t const ow = out->width();
960 int64_t const oh = out->height();
961 int64_t const oc = out->channels();
962
963 /* Choose gaussian sigma parameter according to scale factor. */
964 float const scale_x = (float)img->width() / (float)ow;
965 float const scale_y = (float)img->height() / (float)oh;
966 float const sigma = sigma_factor * std::max(scale_x, scale_y) / 2.0f;
967
968 /* Iterate pixels of dest image and convolute with gaussians on input. */
969 int64_t i = 0;
970 for (int64_t y = 0; y < oh; ++y)
971 for (int64_t x = 0; x < ow; ++x, ++i)
972 {
973 float xf = ((float)x + 0.5f) * scale_x;
974 float yf = ((float)y + 0.5f) * scale_y;
975 for (int64_t c = 0; c < oc; ++c)
976 out->at(i, c) = gaussian_kernel<T>(img, xf, yf, c, sigma);
977 }
978}
979
980/* ---------------------------------------------------------------- */
981
982template <typename T>
983typename Image<T>::Ptr
984crop (typename Image<T>::ConstPtr image, int64_t width, int64_t height,
985 int64_t left, int64_t top, T const* fill_color)
986{
987 if (width < 0 || height < 0 || !image.get())
988 throw std::invalid_argument("Invalid width/height or null image given");
989
990 typename Image<T>::Ptr out(Image<T>::create());
991 out->allocate(width, height, image->channels());
992
993 int64_t const iw = image->width();
994 int64_t const ih = image->height();
995 int64_t const ic = image->channels();
996
997 /* Output image with fill color if new pixels are revealed. */
998 if (left < 0 || top < 0 || left + width > iw || top + height > ih)
999 out->fill_color(fill_color);
1000
1001 /* Check if input and output have no overlap. */
1002 if (left >= iw || left <= -width || top >= ih || top <= -height)
1003 return out;
1004
1005 /* Copy horizontal overlap for each overlapping row. */
1006 int64_t const overlap =
1007 ic * (std::min(iw, left + width) - std::max<int64_t>(0, left));
1008 for (int64_t y = std::max<int64_t>(0, -top);
1009 y < std::min(height, ih - top); ++y)
1010 {
1011 int64_t lookup_y = top + y;
1012 if (lookup_y >= ih)
1013 break;
1014 T* out_ptr = &out->at(left < 0 ? -left : 0, y, 0);
1015 T const* in_ptr = &image->at(left > 0 ? left : 0, lookup_y, 0);
1016 std::copy(in_ptr, in_ptr + overlap, out_ptr);
1017 }
1018
1019 return out;
1020}
1021
1022/* ---------------------------------------------------------------- */
1023
1024template <typename T>
1025typename Image<T>::Ptr
1026blur_gaussian (typename Image<T>::ConstPtr in, float sigma)
1027{
1028 if (in == nullptr)
1029 throw std::invalid_argument("Null image given");
1030
1031 /* Small sigmas result in literally no change. */
1032 if (MATH_EPSILON_EQ(sigma, 0.0f, 0.1f))
1033 return in->duplicate();
1034
1035 int64_t const w = in->width();
1036 int64_t const h = in->height();
1037 int64_t const c = in->channels();
1038 int64_t const ks = std::ceil(sigma * 2.884f); // Cap kernel at 1/128
1039 std::vector<float> kernel(ks + 1);
1040
1041 /* Fill kernel values. */
1042 for (int64_t i = 0; i < ks + 1; ++i)
1043 kernel[i] = math::gaussian((float)i, sigma);
1044
1045#if 1 // Separated kernel implementation
1046
1047 /* Convolve the image in x direction. */
1048 typename Image<T>::Ptr sep(Image<T>::create(w, h, c));
1049 int64_t px = 0;
1050 for (int64_t y = 0; y < h; ++y)
1051 for (int64_t x = 0; x < w; ++x, ++px)
1052 for (int64_t cc = 0; cc < c; ++cc)
1053 {
1054 math::Accum<T> accum(T(0));
1055 for (int64_t i = -ks; i <= ks; ++i)
1056 {
1057 int64_t idx = math::clamp<int64_t>(x + i, 0, w - 1);
1058 accum.add(in->at(y * w + idx, cc), kernel[std::abs(i)]);
1059 }
1060 sep->at(px, cc) = accum.normalized();
1061 }
1062
1063 /* Convolve the image in y direction. */
1064 typename Image<T>::Ptr out(Image<T>::create(w, h, c));
1065 px = 0;
1066 for (int64_t y = 0; y < h; ++y)
1067 for (int64_t x = 0; x < w; ++x, ++px)
1068 for (int64_t cc = 0; cc < c; ++cc)
1069 {
1070 math::Accum<T> accum(T(0));
1071 for (int64_t i = -ks; i <= ks; ++i)
1072 {
1073 int64_t idx = math::clamp<int64_t>(y + i, 0, h - 1);
1074 accum.add(sep->at(idx * w + x, cc), kernel[std::abs(i)]);
1075 }
1076 out->at(px, cc) = accum.normalized();
1077 }
1078
1079#else // Non-separated kernel implementation
1080
1081 typename Image<T>::Ptr out(Image<T>::create(w, h, c));
1082 int64_t px = 0;
1083 for (int64_t y = 0; y < h; ++y)
1084 for (int64_t x = 0; x < w; ++x, ++px)
1085 for (int64_t cc = 0; cc < c; ++cc)
1086 {
1087 math::Accum<T> accum(T(0));
1088
1089 for (int64_t ky = -ks; ky <= ks; ++ky)
1090 for (int64_t kx = -ks; kx <= ks; ++kx)
1091 {
1092 int64_t idx_x = math::clamp(x + kx, 0, w - 1);
1093 int64_t idx_y = math::clamp(y + ky, 0, h - 1);
1094 accum.add(in->at(idx_y * w + idx_x, cc),
1095 kernel[std::abs(kx)] * kernel[std::abs(ky)]);
1096 }
1097
1098 out->at(px, cc) = accum.normalized();
1099 }
1100
1101#endif
1102
1103 return out;
1104}
1105
1106/* ---------------------------------------------------------------- */
1107
1108template <typename T>
1109typename Image<T>::Ptr
1111{
1112 if (in == nullptr)
1113 throw std::invalid_argument("Null image given");
1114
1115 int64_t w(in->width());
1116 int64_t h(in->height());
1117 int64_t c(in->channels());
1118 int64_t wc = w * c;
1119
1120#if 1
1121 /* Super-fast separated kernel implementation. */
1122 typename Image<T>::Ptr sep(Image<T>::create(w, h, c));
1123 math::Accum<T>* accums = new math::Accum<T>[c];
1124 T const* row = &in->at(0);
1125 T* outrow = &sep->at(0);
1126 for (int64_t y = 0; y < h; ++y, row += w * c, outrow += w * c)
1127 {
1128 for (int64_t cc = 0; cc < c; ++cc) // Reset accumulators
1129 accums[cc] = math::Accum<T>(T(0));
1130
1131 for (int i = 0; i < ks; ++i)
1132 for (int64_t cc = 0; cc < c; ++cc)
1133 accums[cc].add(row[i * c + cc], 1.0f);
1134
1135 for (int64_t x = 0; x < w; ++x)
1136 {
1137 if (x + ks < w - 1)
1138 for (int64_t cc = 0; cc < c; ++cc)
1139 accums[cc].add(row[(x + ks) * c + cc], 1.0f);
1140 if (x > ks)
1141 for (int64_t cc = 0; cc < c; ++cc)
1142 accums[cc].sub(row[(x - ks - 1) * c + cc], 1.0f);
1143 for (int64_t cc = 0; cc < c; ++cc)
1144 outrow[x * c + cc] = accums[cc].normalized();
1145 }
1146 }
1147
1148 /* Second filtering pass with kernel in y-direction. */
1149 typename Image<T>::Ptr out(Image<T>::create(w, h, c));
1150 T const* col = &sep->at(0);
1151 T* outcol = &out->at(0);
1152 for (int64_t x = 0; x < w; ++x, col += c, outcol += c)
1153 {
1154 for (int64_t cc = 0; cc < c; ++cc)
1155 accums[cc] = math::Accum<T>(T(0));
1156
1157 for (int i = 0; i < ks; ++i)
1158 for (int64_t cc = 0; cc < c; ++cc)
1159 accums[cc].add(col[i * wc + cc], 1.0f);
1160
1161 for (int64_t y = 0; y < h; ++y)
1162 {
1163 if (y + ks < h - 1)
1164 for (int64_t cc = 0; cc < c; ++cc)
1165 accums[cc].add(col[(y + ks) * wc + cc], 1.0f);
1166 if (y > ks)
1167 for (int64_t cc = 0; cc < c; ++cc)
1168 accums[cc].sub(col[(y - ks - 1) * wc + cc], 1.0f);
1169 for (int64_t cc = 0; cc < c; ++cc)
1170 outcol[y * wc + cc] = accums[cc].normalized();
1171 }
1172 }
1173
1174 delete [] accums;
1175#endif
1176
1177
1178#if 0 // Slow and naive, but simple implementation
1179 typename Image<T>::Ptr out(Image<T>::create(w, h, c));
1180 for (int64_t y = 0; y < h; ++y)
1181 for (int64_t x = 0; x < w; ++x)
1182 for (int64_t cc = 0; cc < c; ++cc)
1183 {
1184 math::Accum<T> accum(T(0));
1185 for (int ky = -ks; ky <= ks; ++ky)
1186 for (int kx = -ks; kx <= ks; ++kx)
1187 {
1188 int64_t idx_x = math::clamp(x + kx, 0, w - 1);
1189 int64_t idx_y = math::clamp(y + ky, 0, h - 1);
1190 accum.add(in->at(idx_x, idx_y, cc), 1.0f);
1191 }
1192 out->at(x, y, cc) = accum.normalized();
1193 }
1194#endif
1195
1196 return out;
1197}
1198
1199/* ---------------------------------------------------------------- */
1200
1201template <typename T>
1202typename Image<T>::Ptr
1204{
1205 if (image == nullptr)
1206 throw std::invalid_argument("Null image given");
1207
1208 int64_t const iw = image->width();
1209 int64_t const ih = image->height();
1210 int64_t const ic = image->channels();
1211 int64_t const ow = type == ROTATE_180 ? iw : ih;
1212 int64_t const oh = type == ROTATE_180 ? ih : iw;
1213
1214 typename Image<T>::Ptr ret(Image<T>::create());
1215 ret->allocate(ow, oh, ic);
1216
1217 int64_t idx = 0;
1218 for (int64_t y = 0; y < ih; ++y)
1219 for (int64_t x = 0; x < iw; ++x, idx += ic)
1220 {
1221 int64_t dx = x;
1222 int64_t dy = y;
1223 switch (type)
1224 {
1225 case ROTATE_180: dx = iw - x - 1; dy = ih - y - 1; break;
1226 case ROTATE_CW: dx = ih - y - 1; dy = x; break;
1227 case ROTATE_CCW: dx = y; dy = iw - x - 1; break;
1228 case ROTATE_SWAP: dx = y; dy = x; break;
1229 default: break;
1230 }
1231
1232 T const* in_pixel = &image->at(idx);
1233 T* out_pixel = &ret->at(dx, dy, 0);
1234 std::copy(in_pixel, in_pixel + ic, out_pixel);
1235 }
1236
1237 return ret;
1238}
1239
1240/* ---------------------------------------------------------------- */
1241
1242template <typename T>
1243typename Image<T>::Ptr
1244rotate (typename Image<T>::ConstPtr image, float angle, T const* fill_color)
1245{
1246 if (image == nullptr)
1247 throw std::invalid_argument("Null image given");
1248
1249 int64_t const w = image->width();
1250 int64_t const h = image->height();
1251 int64_t const c = image->channels();
1252 float const w2 = static_cast<float>(w - 1) / 2.0f;
1253 float const h2 = static_cast<float>(h - 1) / 2.0f;
1254 typename Image<T>::Ptr ret = Image<T>::create(w, h, c);
1255
1256 float const sin_angle = std::sin(-angle);
1257 float const cos_angle = std::cos(-angle);
1258 T* ret_ptr = ret->begin();
1259 for (int64_t y = 0; y < h; ++y)
1260 for (int64_t x = 0; x < w; ++x, ret_ptr += c)
1261 {
1262 float const sample_x =
1263 cos_angle * (x - w2) - sin_angle * (y - h2) + w2;
1264 float const sample_y =
1265 sin_angle * (x - w2) + cos_angle * (y - h2) + h2;
1266 if (sample_x < -0.5f || sample_x > w - 0.5f
1267 || sample_y < -0.5f || sample_y > h - 0.5f)
1268 std::copy(fill_color, fill_color + c, ret_ptr);
1269 else
1270 image->linear_at(sample_x, sample_y, ret_ptr);
1271 }
1272 return ret;
1273}
1274
1275/* ---------------------------------------------------------------- */
1276
1277template <typename T>
1278void
1279flip (typename Image<T>::Ptr image, FlipType type)
1280{
1281 if (type == FLIP_NONE)
1282 return;
1283 if (image == nullptr)
1284 throw std::invalid_argument("Null image given");
1285
1286 bool const fh = type & FLIP_HORIZONTAL;
1287 bool const fv = type & FLIP_VERTICAL;
1288 int64_t const iw = image->width();
1289 int64_t const ih = image->height();
1290 int64_t const ic = image->channels();
1291 int64_t const max_x = fh ? iw / 2 : iw;
1292 int64_t const max_y = fv && !fh ? ih / 2 : ih;
1293
1294 for (int64_t y = 0; y < max_y; ++y)
1295 for (int64_t x = 0; x < max_x; ++x)
1296 {
1297 int64_t dx = (type & FLIP_HORIZONTAL ? iw - x - 1 : x);
1298 int64_t dy = (type & FLIP_VERTICAL ? ih - y - 1 : y);
1299 T* src = &image->at(x, y, 0);
1300 T* dst = &image->at(dx, dy, 0);
1301 std::swap_ranges(src, src + ic, dst);
1302 }
1303}
1304
1305/* ---------------------------------------------------------------- */
1306
1307template <typename T>
1308inline T
1310{
1311 return *std::max_element(v, v + 3);
1312}
1313
1314template <typename T>
1315inline T
1317{
1318 T const* max = std::max_element(v, v + 3);
1319 T const* min = std::min_element(v, v + 3);
1320 return math::interpolate(*max, *min, 0.5f, 0.5f);
1321}
1322
1323template <typename T>
1324inline T
1326{
1327 return math::interpolate(v[0], v[1], v[2], 0.21f, 0.72f, 0.07f);
1328}
1329
1330template <typename T>
1331inline T
1333{
1334 return math::interpolate(v[0], v[1], v[2], 0.30f, 0.59f, 0.11f);
1335}
1336
1337template <typename T>
1338inline T
1340{
1341 float third(1.0f / 3.0f);
1342 return math::interpolate(v[0], v[1], v[2], third, third, third);
1343}
1344
1345/* ---------------------------------------------------------------- */
1346
1347template <typename T>
1348typename Image<T>::Ptr
1350{
1351 if (img == nullptr)
1352 throw std::invalid_argument("Null image given");
1353
1354 int64_t ic = img->channels();
1355 if (ic != 3 && ic != 4)
1356 throw std::invalid_argument("Image must be RGB or RGBA");
1357
1358 bool has_alpha = (ic == 4);
1359
1360 typename Image<T>::Ptr out(Image<T>::create());
1361 out->allocate(img->width(), img->height(), 1 + has_alpha);
1362
1363 typedef T(*DesaturateFunc)(T const*);
1364 DesaturateFunc func;
1365 switch (type)
1366 {
1367 case DESATURATE_MAXIMUM: func = desaturate_maximum<T>; break;
1368 case DESATURATE_LIGHTNESS: func = desaturate_lightness<T>; break;
1369 case DESATURATE_LUMINOSITY: func = desaturate_luminosity<T>; break;
1370 case DESATURATE_LUMINANCE: func = desaturate_luminance<T>; break;
1371 case DESATURATE_AVERAGE: func = desaturate_average<T>; break;
1372 default: throw std::invalid_argument("Invalid desaturate type");
1373 }
1374
1375 int64_t outpos = 0;
1376 int64_t inpos = 0;
1377 int64_t pixels = img->get_pixel_amount();
1378 for (int64_t i = 0; i < pixels; ++i)
1379 {
1380 T const* v = &img->at(inpos);
1381 out->at(outpos) = func(v);
1382
1383 if (has_alpha)
1384 out->at(outpos + 1) = img->at(inpos + 3);
1385
1386 outpos += 1 + has_alpha;
1387 inpos += 3 + has_alpha;
1388 }
1389
1390 return out;
1391}
1392
1393/* ---------------------------------------------------------------- */
1394
1395template <typename T>
1396typename Image<T>::Ptr
1398{
1399 if (image == nullptr)
1400 throw std::invalid_argument("Null image given");
1401
1402 int64_t const ic = image->channels();
1403 if (ic != 1 && ic != 2)
1404 throw std::invalid_argument("Image must be in G or GA");
1405
1406 bool const has_alpha = (ic == 2);
1407
1408 typename Image<T>::Ptr out(Image<T>::create());
1409 out->allocate(image->width(), image->height(), 3 + has_alpha);
1410
1411 int64_t pixels = image->get_pixel_amount();
1412 for (int64_t i = 0; i < pixels; ++i)
1413 {
1414 out->at(i, 0) = image->at(i, 0);
1415 out->at(i, 1) = image->at(i, 0);
1416 out->at(i, 2) = image->at(i, 0);
1417 if (has_alpha)
1418 out->at(i, 3) = image->at(i, 1);
1419 }
1420
1421 return out;
1422}
1423
1424/* ---------------------------------------------------------------- */
1425
1426template <typename T>
1427void
1429{
1430 int64_t const channels = img->channels();
1431 if (channels != 2 && channels != 4)
1432 throw std::invalid_argument("Image must be in GA or RGBA");
1433 img->delete_channel(channels - 1);
1434}
1435
1436/* ---------------------------------------------------------------- */
1437
1438template <typename T>
1439typename mve::Image<T>::Ptr
1441{
1442 int64_t const width = img->width();
1443 int64_t const height = img->height();
1444 int64_t const chans = img->channels();
1445 int64_t const row_stride = width * chans;
1446
1447 double const max_value = static_cast<double>(std::numeric_limits<T>::max());
1449 (width, height, chans);
1450 T* out_ptr = out->get_data_pointer();
1451
1452 int64_t pos = 0;
1453 for (int64_t y = 0; y < height; ++y)
1454 for (int64_t x = 0; x < width; ++x, pos += chans)
1455 {
1456 if (y == 0 || y == height - 1 || x == 0 || x == width - 1)
1457 {
1458 std::fill(out_ptr + pos, out_ptr + pos + chans, T(0));
1459 continue;
1460 }
1461
1462 for (int64_t cc = 0; cc < chans; ++cc)
1463 {
1464 int64_t const i = pos + cc;
1465 double gx = 1.0 * (double)img->at(i + chans - row_stride)
1466 - 1.0 * (double)img->at(i - chans - row_stride)
1467 + 2.0 * (double)img->at(i + chans)
1468 - 2.0 * (double)img->at(i - chans)
1469 + 1.0 * (double)img->at(i + chans + row_stride)
1470 - 1.0 * (double)img->at(i - chans + row_stride);
1471 double gy = 1.0 * (double)img->at(i + row_stride - chans)
1472 - 1.0 * (double)img->at(i - row_stride - chans)
1473 + 2.0 * (double)img->at(i + row_stride)
1474 - 2.0 * (double)img->at(i - row_stride)
1475 + 1.0 * (double)img->at(i + row_stride + chans)
1476 - 1.0 * (double)img->at(i - row_stride + chans);
1477 double g = std::sqrt(gx * gx + gy * gy);
1478 out_ptr[i] = static_cast<T>(std::min(max_value, g));
1479 }
1480 }
1481
1482 return out;
1483}
1484
1485/* ---------------------------------------------------------------- */
1486
1487template <typename T>
1488typename Image<T>::Ptr
1490{
1491 int64_t const iw = i1->width();
1492 int64_t const ih = i1->height();
1493 int64_t const ic = i1->channels();
1494
1495 if (i1 == nullptr || i2 == nullptr)
1496 throw std::invalid_argument("Null image given");
1497
1498 if (iw != i2->width() || ih != i2->height() || ic != i2->channels())
1499 throw std::invalid_argument("Image dimensions do not match");
1500
1501 typename Image<T>::Ptr out(Image<T>::create());
1502 out->allocate(iw, ih, ic);
1503
1504 // FIXME: That needs speedup! STL algo with pointers.
1505 for (int64_t i = 0; i < i1->get_value_amount(); ++i)
1506 out->at(i) = i1->at(i) - i2->at(i);
1507
1508 return out;
1509}
1510
1511/* ---------------------------------------------------------------- */
1512
1513template <typename T>
1514typename Image<T>::Ptr
1516{
1517 int64_t const iw = i1->width();
1518 int64_t const ih = i1->height();
1519 int64_t const ic = i1->channels();
1520
1521 if (i1 == nullptr || i2 == nullptr)
1522 throw std::invalid_argument("Null image given");
1523
1524 if (iw != i2->width() || ih != i2->height() || ic != i2->channels())
1525 throw std::invalid_argument("Image dimensions do not match");
1526
1527 typename Image<T>::Ptr out(Image<T>::create());
1528 out->allocate(iw, ih, ic);
1529
1530 for (int64_t i = 0; i < i1->get_value_amount(); ++i)
1531 {
1532 if (i1->at(i) < i2->at(i))
1533 out->at(i) = i2->at(i) - i1->at(i);
1534 else
1535 out->at(i) = i1->at(i) - i2->at(i);
1536 }
1537
1538 return out;
1539}
1540
1541/* ---------------------------------------------------------------- */
1542
1543template <typename T>
1544void
1545gamma_correct (typename Image<T>::Ptr image, T const& power)
1546{
1548 std::for_each(image->begin(), image->end(), f);
1549}
1550
1551/* ---------------------------------------------------------------- */
1552
1553template <typename T>
1554void
1556{
1557 static_assert(std::is_floating_point<T>::value,
1558 "Only implemented for floating point images");
1559 int64_t const num_values = image->get_value_amount();
1560 for (int64_t i = 0; i < num_values; i++)
1561 {
1562 if (image->at(i) <= T(0.0031308))
1563 image->at(i) *= T(12.92);
1564 else
1565 {
1566 T corrected = std::pow(image->at(i), T(1.0)/T(2.4));
1567 image->at(i) = T(1.055) * corrected - T(0.055);
1568 }
1569 }
1570}
1571
1572/* ---------------------------------------------------------------- */
1573
1574template <typename T>
1575void
1577{
1578 static_assert(std::is_floating_point<T>::value,
1579 "Only implemented for floating point images");
1580 int64_t const num_values = image->get_value_amount();
1581 for (int64_t i = 0; i < num_values; i++)
1582 {
1583 if (image->at(i) <= T(0.04045))
1584 image->at(i) /= T(12.92);
1585 else
1586 {
1587 T base = (image->at(i) + T(0.055)) / T(1.055);
1588 image->at(i) = std::pow(base, T(2.4));
1589 }
1590 }
1591}
1592
1593/* ---------------------------------------------------------------- */
1594
1595template <typename T_IN, typename T_OUT>
1596typename Image<T_OUT>::Ptr
1598{
1599 if (image == nullptr)
1600 throw std::invalid_argument("Null image given");
1601
1602 int64_t const width = image->width();
1603 int64_t const height = image->height();
1604 int64_t const chans = image->channels();
1605 int64_t const row_stride = width * chans;
1606
1608 ret->allocate(width, height, chans);
1609
1610 /* Input image row and destination image rows. */
1611 std::vector<T_OUT> zeros(row_stride, T_OUT(0));
1612 T_IN const* inrow = image->get_data_pointer();
1613 T_OUT* dest = ret->get_data_pointer();
1614 T_OUT* prev = &zeros[0];
1615
1616 /*
1617 * I(x,y) = i(x,y) + I(x-1,y) + I(x,y-1) - I(x-1,y-1)
1618 */
1619 for (int64_t y = 0; y < height; ++y)
1620 {
1621 /* Calculate first pixel in row. */
1622 for (int64_t cc = 0; cc < chans; ++cc)
1623 dest[cc] = static_cast<T_OUT>(inrow[cc]) + prev[cc];
1624 /* Calculate all following pixels in row. */
1625 for (int64_t i = chans; i < row_stride; ++i)
1626 dest[i] = inrow[i] + prev[i] + dest[i - chans] - prev[i - chans];
1627
1628 prev = dest;
1629 dest += row_stride;
1630 inrow += row_stride;
1631 }
1632
1633 return ret;
1634}
1635
1636/* ---------------------------------------------------------------- */
1637
1638template <typename T>
1639T
1641 int64_t x1, int64_t y1, int64_t x2, int64_t y2, int64_t cc)
1642{
1643 int64_t const row_stride = sat->width() * sat->channels();
1644 int64_t const channels = sat->channels();
1645 T ret = sat->at(y2 * row_stride + x2 * channels + cc);
1646 if (x1 > 0)
1647 ret -= sat->at(y2 * row_stride + (x1-1) * channels + cc);
1648 if (y1 > 0)
1649 ret -= sat->at((y1-1) * row_stride + x2 * channels + cc);
1650 if (x1 > 0 && y1 > 0)
1651 ret += sat->at((y1-1) * row_stride + (x1-1) * channels + cc);
1652 return ret;
1653}
1654
1655/* ---------------------------------------------------------------- */
1656
1657template <typename T>
1658typename Image<T>::Ptr
1660 int64_t thumb_width, int64_t thumb_height)
1661{
1662 int64_t const width = image->width();
1663 int64_t const height = image->height();
1664 float const image_aspect = static_cast<float>(width) / height;
1665 float const thumb_aspect = static_cast<float>(thumb_width) / thumb_height;
1666
1667 int64_t rescale_width, rescale_height;
1668 int64_t crop_left, crop_top;
1669 if (image_aspect > thumb_aspect)
1670 {
1671 rescale_width = std::ceil(thumb_height * image_aspect);
1672 rescale_height = thumb_height;
1673 crop_left = (rescale_width - thumb_width) / 2;
1674 crop_top = 0;
1675 }
1676 else
1677 {
1678 rescale_width = thumb_width;
1679 rescale_height = std::ceil(thumb_width / image_aspect);
1680 crop_left = 0;
1681 crop_top = (rescale_height - thumb_height) / 2;
1682 }
1683
1684 typename mve::Image<T>::Ptr thumb = mve::image::rescale<T>(image,
1685 mve::image::RESCALE_LINEAR, rescale_width, rescale_height);
1686 thumb = mve::image::crop<T>(thumb, thumb_width, thumb_height,
1687 crop_left, crop_top, nullptr);
1688
1689 return thumb;
1690}
1691
1692/* ---------------------------------------------------------------- */
1693
1694template <typename T>
1695typename Image<T>::Ptr
1696image_undistort_msps (typename Image<T>::ConstPtr img, double k0, double k1)
1697{
1698 int64_t const width = img->width();
1699 int64_t const height = img->height();
1700 int64_t const chans = img->channels();
1701 int64_t const D = std::max(width, height);
1702
1703 double const width_half = static_cast<double>(width) / 2.0;
1704 double const height_half = static_cast<double>(height) / 2.0;
1705
1706 typename Image<T>::Ptr out = Image<T>::create(width, height, chans);
1707 out->fill(T(0));
1708 T* out_ptr = out->get_data_pointer();
1709
1710 for (int64_t y = 0; y < height; y++)
1711 for (int64_t x = 0; x < width; x++, out_ptr += chans)
1712 {
1713 double fx = static_cast<double>(x) - width_half;
1714 double fy = static_cast<double>(y) - height_half;
1715 double const s1 = D * D + k1 * (fx * fx + fy * fy);
1716 double const s2 = D * D + k0 * (fx * fx + fy * fy);
1717 double const factor = s1 / s2;
1718 fx = fx * factor + width_half;
1719 fy = fy * factor + height_half;
1720
1721 if (fx < -0.5 || fx > width - 0.5 || fy < -0.5 || fy > height - 0.5)
1722 continue;
1723 img->linear_at(fx, fy, out_ptr);
1724 }
1725
1726 return out;
1727}
1728
1729/* ---------------------------------------------------------------- */
1730
1731template <typename T>
1732typename Image<T>::Ptr
1734 double focal_length, double k2, double k4)
1735{
1736 if (img == nullptr)
1737 throw std::invalid_argument("Null image given");
1738
1739 if (k2 == 0.0 && k4 == 0.0)
1740 return img->duplicate();
1741
1742 int64_t const width = img->width();
1743 int64_t const height = img->height();
1744 int64_t const chans = img->channels();
1745 typename Image<T>::Ptr out = Image<T>::create(width, height, chans);
1746 out->fill(T(0));
1747 T* out_ptr = out->get_data_pointer();
1748
1749 double const fwidth2 = static_cast<double>(width) / 2.0;
1750 double const fheight2 = static_cast<double>(height) / 2.0;
1751 double const fnorm = static_cast<double>(std::max(width, height));
1752 for (int64_t y = 0; y < height; ++y)
1753 for (int64_t x = 0; x < width; ++x, out_ptr += chans)
1754 {
1755 double const fx = (static_cast<double>(x) + 0.5 - fwidth2) / fnorm;
1756 double const fy = (static_cast<double>(y) + 0.5 - fheight2) / fnorm;
1757 double const rd = (fx * fx + fy * fy) / MATH_POW2(focal_length);
1758 double const rd_factor = 1.0 + k2 * rd + k4 * rd * rd;
1759 double const dist_x = fx * rd_factor;
1760 double const dist_y = fy * rd_factor;
1761 float const ix = dist_x * fnorm + fwidth2 - 0.5;
1762 float const iy = dist_y * fnorm + fheight2 - 0.5;
1763 if (ix < -0.5 || ix > width - 0.5 || iy < -0.5 || iy > height - 0.5)
1764 continue;
1765 img->linear_at(ix, iy, out_ptr);
1766 }
1767
1768 return out;
1769}
1770
1771/* ---------------------------------------------------------------- */
1772
1773template <typename T>
1774typename Image<T>::Ptr
1776 double focal_length, double k1)
1777{
1778 if (img == nullptr)
1779 throw std::invalid_argument("Null image given");
1780
1781 if (k1 == 0.0)
1782 return img->duplicate();
1783
1784 int64_t const width = img->width();
1785 int64_t const height = img->height();
1786 int64_t const chans = img->channels();
1787
1788 /*
1789 * The image coordinates must be normalized before the distortion
1790 * model is applied. The image coordinates are first centered at
1791 * the origin and then scaled w.r.t. the focal length in pixel.
1792 */
1793 double const norm = focal_length * std::max(width, height);
1794 double const width_half = static_cast<double>(width) / 2.0;
1795 double const height_half = static_cast<double>(height) / 2.0;
1796
1797 typename Image<T>::Ptr out = Image<T>::create(width, height, chans);
1798 out->fill(T(0));
1799 T* out_ptr = out->begin();
1800
1801 for (int64_t y = 0; y < height; ++y)
1802 {
1803 for (int64_t x = 0; x < width; ++x, out_ptr += chans)
1804 {
1805 double fx = (static_cast<double>(x) - width_half) / norm;
1806 double fy = (static_cast<double>(y) - height_half) / norm;
1807 if (fy == 0.0)
1808 fy = 1e-10;
1809
1810 double const t2 = fy * fy;
1811 double const t3 = t2 * t2 * t2;
1812 double const t4 = fx * fx;
1813 double const t7 = k1 * (t2 + t4);
1814
1815 if (k1 > 0.0)
1816 {
1817 double const t8 = 1.0 / t7;
1818 double const t10 = t3 / (t7 * t7);
1819 double const t14 = std::sqrt(t10 * (0.25 + t8 / 27.0));
1820 double const t15 = t2 * t8 * fy * 0.5;
1821 double const t17 = std::pow(t14 + t15, 1.0/3.0);
1822 double const t18 = t17 - t2 * t8 / (t17 * 3.0);
1823 fx = t18 * fx / fy;
1824 fy = t18;
1825 }
1826 else
1827 {
1828 double const t9 = t3 / (t7 * t7 * 4.0);
1829 double const t11 = t3 / (t7 * t7 * t7 * 27.0);
1830 std::complex<double> const t12 = t9 + t11;
1831 std::complex<double> const t13 = std::sqrt(t12);
1832 double const t14 = t2 / t7;
1833 double const t15 = t14 * fy * 0.5;
1834 std::complex<double> const t16 = t13 + t15;
1835 std::complex<double> const t17 = std::pow(t16, 1.0/3.0);
1836 std::complex<double> const t18 = (t17 + t14 / (t17 * 3.0))
1837 * std::complex<double>(0.0, std::sqrt(3.0));
1838 std::complex<double> const t19 = -0.5 * (t17 + t18)
1839 + t14 / (t17 * 6.0);
1840 fx = t19.real() * fx / fy;
1841 fy = t19.real();
1842 }
1843
1844 fx = fx * norm + width_half;
1845 fy = fy * norm + height_half;
1846
1847 if (fx < -0.5 || fx > width - 0.5 || fy < -0.5 || fy > height - 0.5)
1848 continue;
1849 img->linear_at(fx, fy, out_ptr);
1850 }
1851 }
1852
1853 return out;
1854}
1855
1858
1859#endif /* MVE_IMAGE_TOOLS_HEADER */
Accumulator class that operates on arbitrary types.
Definition accum.h:35
void add(T const &value, float weight)
Adds the weighted given value to the internal value.
Definition accum.h:86
T normalized(float weight) const
Returns a normalized version of the internal value, i.e.
Definition accum.h:102
int64_t height(void) const
Returns the height of the image.
Definition image_base.h:207
int64_t channels(void) const
Returns the amount of channels in the image.
Definition image_base.h:213
int64_t width(void) const
Returns the width of the image.
Definition image_base.h:201
Multi-channel image class of arbitrary but homogenous data type.
Definition image.h:40
std::shared_ptr< Image< T > > Ptr
Definition image.h:42
std::shared_ptr< Image< T > const > ConstPtr
Definition image.h:43
T const & at(int64_t index) const
Linear indexing of image data.
Definition image.h:307
T linear_at(float x, float y, int64_t channel) const
Linear interpolation (more expensive) for a single color channel.
Definition image.h:409
Ptr duplicate(void) const
Duplicates the image.
Definition image.h:212
static Ptr create(void)
Smart pointer image constructor.
Definition image.h:191
void fill_color(T const *color)
Fills every pixel of the image with the given color.
Definition image.h:219
void delete_channel(int64_t channel)
Deletes a channel from the image.
Definition image.h:288
int64_t get_pixel_amount(void) const
Returns the amount of pixels in the image (w * h).
Definition image_base.h:500
int64_t get_value_amount(void) const
Returns the amount of values in the image (w * h * c).
Definition image_base.h:507
void fill(T const &value)
Fills the data with a constant value.
Definition image_base.h:423
T * begin(void)
Returns data pointer to beginning.
Definition image_base.h:472
T * end(void)
Returns data pointer to end.
Definition image_base.h:486
T const * get_data_pointer(void) const
Returns the data pointer.
Definition image_base.h:454
void allocate(int64_t width, int64_t height, int64_t chans)
Allocates new image space, clearing previous content.
Definition image_base.h:395
#define MATH_EPSILON_EQ(x, v, eps)
Definition defines.h:96
#define MATH_POW2(x)
Definition defines.h:68
#define MVE_IMAGE_NAMESPACE_END
Definition defines.h:17
#define MVE_NAMESPACE_BEGIN
Definition defines.h:13
#define MVE_IMAGE_NAMESPACE_BEGIN
Definition defines.h:16
#define MVE_NAMESPACE_END
Definition defines.h:14
T const & clamp(T const &v, T const &min=T(0), T const &max=T(1))
Returns value 'v' clamped to the interval specified by 'min' and 'max'.
Definition functions.h:204
T gaussian(T const &x, T const &sigma)
Gaussian function g(x) = exp( -1/2 * (x/sigma)^2 ).
Definition functions.h:36
T interpolate(T const &v1, float w1)
Generic interpolation (weighting) of a single value.
Definition functions.h:80
T gaussian_xx(T const &xx, T const &sigma)
Gaussian function that expects x to be squared.
Definition functions.h:49
Image< T >::Ptr rescale(typename Image< T >::ConstPtr image, RescaleInterpolation interp, int64_t width, int64_t height)
Returns a rescaled version of 'image' with dimensions 'width' times 'height' using 'interp' for value...
RotateType
Image rotation type.
@ ROTATE_SWAP
Exchanges x- and y-axis.
@ ROTATE_CW
Clock wise rotation.
@ ROTATE_180
180 degree rotation
@ ROTATE_CCW
Counter-clock wise rotation.
T desaturate_luminance(T const *v)
T desaturate_lightness(T const *v)
Image< T >::Ptr rescale_half_size(typename Image< T >::ConstPtr image)
Returns a rescaled version of image, scaled by factor 1/2, by grouping blocks of 2x2 pixel into one p...
FloatImage::Ptr byte_to_float_image(ByteImage::ConstPtr image)
Converts a given byte image to a float image.
void gamma_correct_inv_srgb(typename Image< T >::Ptr image)
Applies inverse gamma correction to float/double (in-place) images with nonlinear R'G'B' values in th...
ByteImage::Ptr int_to_byte_image(IntImage::ConstPtr image)
Convertes a given int image to a byte image.
Image< T >::Ptr difference(typename Image< T >::ConstPtr i1, typename Image< T >::ConstPtr i2)
Creates a difference image by computing the absolute difference per value.
void rescale_nearest(typename Image< T >::ConstPtr in, typename Image< T >::Ptr out)
Rescales image 'in' using nearest neighbor.
Image< T >::Ptr rescale_double_size(typename Image< T >::ConstPtr img)
Returns a rescaled version of the image, upscaled with linear interpolation by factor 2.
Image< T >::Ptr subtract(typename Image< T >::ConstPtr i1, typename Image< T >::ConstPtr i2)
Subtracts two images to create the signed difference between the values.
T desaturate_maximum(T const *v)
Image< T >::Ptr rescale_half_size_gaussian(typename Image< T >::ConstPtr image, float sigma=0.866025403784439f)
Returns a rescaled version of the image, scaled with a gaussian approximation by factor 1/2.
DoubleImage::Ptr byte_to_double_image(ByteImage::ConstPtr image)
Converts a given byte image to a double image.
ByteImage::Ptr float_to_byte_image(FloatImage::ConstPtr image, float vmin, float vmax)
Converts a given float image to a byte image.
Image< T >::Ptr blur_boxfilter(typename Image< T >::ConstPtr in, int ks)
Blurs the image using a box filter of integer size 'ks'.
Image< T >::Ptr desaturate(typename Image< T >::ConstPtr image, DesaturateType type)
Desaturates an RGB or RGBA image to G or GA respectively.
T desaturate_average(T const *v)
Image< T >::Ptr blur_gaussian(typename Image< T >::ConstPtr in, float sigma)
Blurs the image using a gaussian convolution kernel.
Image< DST >::Ptr type_to_type_image(typename Image< SRC >::ConstPtr image)
Generic conversion between image types without scaling or clamping.
void rescale_gaussian(typename Image< T >::ConstPtr in, typename Image< T >::Ptr out, float sigma_factor=1.0f)
Rescales image 'in' using a gaussian kernel mask.
Image< T >::Ptr image_undistort_k2k4(typename Image< T >::ConstPtr img, double focal_length, double k2, double k4)
Undistorts the input image given the focal length of the image and two undistortion parameters.
void rescale_linear(typename Image< T >::ConstPtr in, typename Image< T >::Ptr out)
Rescales image 'in' using linear interpolation.
mve::Image< T >::Ptr sobel_edge(typename mve::Image< T >::ConstPtr img)
Implementation of the Sobel operator.
RescaleInterpolation
Rescale interpolation type.
@ RESCALE_GAUSSIAN
Not suited for byte images.
FlipType
Image flipping type.
void float_image_normalize(FloatImage::Ptr image)
Normalizes a float image IN-PLACE such that all values are [0, 1].
Image< T >::Ptr image_undistort_msps(typename Image< T >::ConstPtr img, double k0, double k1)
Undistorts the input image given the two undistortion parameters.
ByteImage::Ptr double_to_byte_image(DoubleImage::ConstPtr image, double vmin, double vmax)
Converts a given double image to a byte image.
Image< T >::Ptr rescale_half_size_subsample(typename Image< T >::ConstPtr image)
Returns a rescaled version of the image by subsampling every second column and row.
void find_min_max_value(typename mve::Image< T >::ConstPtr image, T *vmin, T *vmax)
Finds the smallest and largest value in the given image.
void gamma_correct(ByteImage::Ptr image, float power)
Applies fast gamma correction to byte image using a lookup table.
DesaturateType
Desaturaturation type.
@ DESATURATE_LIGHTNESS
(max(R,G,B) + min(R,G,B)) * 1/2
@ DESATURATE_AVERAGE
(R + G + B) * 1/3
@ DESATURATE_LUMINANCE
0.30 * R + 0.59 * G + 0.11 * B
@ DESATURATE_LUMINOSITY
0.21 * R + 0.72 * G + 0.07 * B
@ DESATURATE_MAXIMUM
max(R,G,B)
Image< T >::Ptr crop(typename Image< T >::ConstPtr image, int64_t width, int64_t height, int64_t left, int64_t top, T const *fill_color)
Returns a sub-image by cropping against a rectangular region.
Image< T >::Ptr rescale_double_size_supersample(typename Image< T >::ConstPtr img)
Returns a rescaled version of the image, upscaled with linear interpolation.
Image< T >::Ptr create_thumbnail(typename Image< T >::ConstPtr image, int64_t thumb_width, int64_t thumb_height)
Creates a thumbnail of the given size by first rescaling the image and then cropping to fill the thum...
Image< T >::Ptr rotate(typename Image< T >::ConstPtr image, RotateType type)
Returns a rotated copy of the given image either rotated clock wise, counter-clock wise,...
void reduce_alpha(typename mve::Image< T >::Ptr img)
Reduce alpha: Reduces RGBA or GA images to RGB or G images.
void gamma_correct_srgb(typename Image< T >::Ptr image)
Applies gamma correction to float/double images (in-place) with linear RGB values in range [0,...
T desaturate_luminosity(T const *v)
void flip(typename Image< T >::Ptr image, FlipType type)
Flips the given image either horizontally, vertically or both IN-PLACE.
ByteImage::Ptr raw_to_byte_image(RawImage::ConstPtr image, uint16_t vmin, uint16_t vmax)
Converts a given raw image to a byte image.
Image< T_OUT >::Ptr integral_image(typename Image< T_IN >::ConstPtr image)
Calculates the integral image (or summed area table) for the input image.
Image< T >::Ptr expand_grayscale(typename Image< T >::ConstPtr image)
Expands a gray image (one or two channels) to an RGB or RGBA image.
FloatImage::Ptr raw_to_float_image(RawImage::ConstPtr image)
Converts a given raw image to a float image.
Image< T >::Ptr image_undistort_vsfm(typename Image< T >::ConstPtr img, double focal_length, double k1)
Undistorts the input image given the focal length of the image and a single distortion parameter.
T integral_image_area(typename Image< T >::ConstPtr sat, int64_t x1, int64_t y1, int64_t x2, int64_t y2, int64_t cc=0)
Sums over the rectangle defined by A=(x1,y1) and B=(x2,y2) on the given SAT for channel cc.
T gaussian_kernel(typename Image< T >::ConstPtr img, float x, float y, int64_t c, float sigma)
void swap(mve::Image< T > &a, mve::Image< T > &b)
Specialization of std::swap for efficient image swapping.
Definition image.h:478
for-each functor: raises each operand to the power of constant value.
Definition algo.h:259