7
\$\begingroup\$

This is a follow-up question for Parallel 3D Discrete Cosine Transformation Implementation in Matlab and Operator overloading in Image class implementation in C++. I am trying to implement 3D Discrete Cosine Transformation calculation in C++. The formula of 3D Discrete Cosine Transformation is as follows.

The 3D discrete cosine transformation $$X(k_{1}, k_{2}, k_{3})$$ of size $$N_{1} \times N_{2} \times N_{3}$$ is

\begin{equation} \begin{split} {X(k_{1}, k_{2}, k_{3})} = {\frac {8}{N_{1} N_{2} N_{3}}} \epsilon_{k_{1}} \epsilon_{k_{2}} \epsilon_{k_{3}} \sum_{{n_1 = 0}}^{N_1 - 1} \sum_{{n_2 = 0}}^{N_2 - 1} \sum_{{n_3 = 0}}^{N_3 - 1} x(n_{1}, n_{2}, n_{3}) \\ \times \cos({\frac {\pi}{2N_{1}} (2n_{1} + 1)k_{1}}) \\ \times \cos({\frac {\pi}{2N_{2}} (2n_{2} + 1)k_{2}}) \\ \times \cos({\frac {\pi}{2N_{3}} (2n_{3} + 1)k_{3}}) \end{split} \label{eq:3DDCTMainFormula} \end{equation}

where

\begin{equation} \begin{split} k_{1} = 0, 1, \dots, N_{1} - 1 \\ k_{2} = 0, 1, \dots, N_{2} - 1 \\ k_{3} = 0, 1, \dots, N_{3} - 1 \\ \epsilon_{k_{i}} = \begin{cases} \frac{1}{\sqrt{2}} & \text{for $k_{i} = 0$} \\ 1 & \text{otherwise} \end{cases} i = 1, 2, 3 \end{split} \label{eq:3DDCTMainFormulaDetail} \end{equation}

The experimental implementation

  • dct3_detail template function implementation:

    template<arithmetic ElementT = double, arithmetic OutputT = ElementT>
    constexpr static Image<ElementT> dct3_detail(std::vector<Image<ElementT>> input, int plane_index)
    {
        std::size_t N3 = input.size();
        OutputT alpha1 = (plane_index == 0) ? (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))) : (static_cast<OutputT>(1.0));
        auto output = Image<OutputT>(input[plane_index].getWidth(), input[plane_index].getHeight());
        for (std::size_t y = 0; y < output.getHeight(); y++)
        {
            OutputT alpha2 = (y == 0) ? (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))) : (static_cast<OutputT>(1.0));
            for (std::size_t x = 0; x < output.getWidth(); x++)
            {
                OutputT sum{};
                OutputT alpha3 = (x == 0) ? (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))) : (static_cast<OutputT>(1.0));
                for (std::size_t inner_z = 0; inner_z < N3; inner_z++)
                {
                    auto plane = input[inner_z];
                    auto N1 = static_cast<OutputT>(plane.getWidth());
                    auto N2 = static_cast<OutputT>(plane.getHeight());
                    for (std::size_t inner_y = 0; inner_y < plane.getHeight(); inner_y++)
                    {
                        for (std::size_t inner_x = 0; inner_x < plane.getWidth(); inner_x++)
                        {
                            auto l1 = (std::numbers::pi / (2 * N1) * (2 * static_cast<OutputT>(inner_x) + 1) * x);
                            auto l2 = (std::numbers::pi / (2 * N2) * (2 * static_cast<OutputT>(inner_y) + 1) * y);
                            auto l3 = (std::numbers::pi / (2 * static_cast<OutputT>(N3)) * (2 * static_cast<OutputT>(inner_z) + 1) * static_cast<OutputT>(plane_index));
                            sum += static_cast<OutputT>(plane.at(inner_x, inner_y)) *
                                std::cos(l1) * std::cos(l2) * std::cos(l3);
                        }
                    }
                }
                auto N1 = static_cast<OutputT>(input[0].getWidth());
                auto N2 = static_cast<OutputT>(input[0].getHeight());
                output.at(x, y) = 8 * alpha1 * alpha2 * alpha3 * sum / (N1 * N2 * N3);
            }
        }
        return output;
    }
    
  • dct3 template function implementation:

    template<arithmetic ElementT = double, arithmetic OutputT = ElementT>
    constexpr static std::vector<Image<ElementT>> dct3(std::vector<Image<ElementT>> input)
    {
        std::vector<Image<ElementT>> output;
        for (std::size_t i = 0; i < input.size(); i++)
        {
            output.push_back(dct3_detail(input, i));
        }
        return output;
    }
    
  • arithmetic concept:

    template <typename T>
    concept arithmetic = std::is_arithmetic_v<T>;
    

Full Testing Code

The full tests for dct3 template function:

#include <algorithm>
#include <array>
#include <cassert>
#include <chrono>
#include <cmath>
#include <concepts>
#include <exception>
#include <execution>
#include <fstream>
#include <functional>
#include <iostream>
#include <iterator>
#include <numbers>
#include <numeric>
#include <ranges>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>

using BYTE = unsigned char;

struct RGB
{
    BYTE channels[3];
};

using GrayScale = BYTE;

namespace TinyDIP
{
    #define is_size_same(x, y) {assert(x.getWidth() == y.getWidth()); assert(x.getHeight() == y.getHeight());}

    //  Reference: https://stackoverflow.com/a/58067611/6667035
    template <typename T>
    concept arithmetic = std::is_arithmetic_v<T>;

    template <typename ElementT>
    class Image
    {
    public:
        Image() = default;

        Image(const std::size_t width, const std::size_t height):
            width(width),
            height(height),
            image_data(width * height) { }

        Image(const std::size_t width, const std::size_t height, const ElementT initVal):
            width(width),
            height(height),
            image_data(width * height, initVal) {}

        Image(std::vector<ElementT> input, std::size_t newWidth, std::size_t newHeight):
            width(newWidth),
            height(newHeight)
        {
            if (input.size() != newWidth * newHeight)
            {
                throw std::runtime_error("Image data input and the given size are mismatched!");
            }
            image_data = std::move(input);
        }

        constexpr ElementT& at(const unsigned int x, const unsigned int y)
        { 
            checkBoundary(x, y);
            return image_data[y * width + x];
        }

        constexpr ElementT const& at(const unsigned int x, const unsigned int y) const
        {
            checkBoundary(x, y);
            return image_data[y * width + x];
        }

        constexpr std::size_t getWidth() const
        {
            return width;
        }

        constexpr std::size_t getHeight() const noexcept
        {
            return height;
        }

        constexpr auto getSize() noexcept
        {
            return std::make_tuple(width, height);
        }

        std::vector<ElementT> const& getImageData() const { return image_data; }      //  expose the internal data

        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        {
            for (std::size_t y = 0; y < height; ++y)
            {
                for (std::size_t x = 0; x < width; ++x)
                {
                    //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                    os << +at(x, y) << separator;
                }
                os << "\n";
            }
            os << "\n";
            return;
        }

        //  Enable this function if ElementT = RGB
        void print(std::string separator = "\t", std::ostream& os = std::cout) const
        requires(std::same_as<ElementT, RGB>)
        {
            for (std::size_t y = 0; y < height; ++y)
            {
                for (std::size_t x = 0; x < width; ++x)
                {
                    os << "( ";
                    for (std::size_t channel_index = 0; channel_index < 3; ++channel_index)
                    {
                        //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                        os << +at(x, y).channels[channel_index] << separator;
                    }
                    os << ")" << separator;
                }
                os << "\n";
            }
            os << "\n";
            return;
        }

        friend std::ostream& operator<<(std::ostream& os, const Image<ElementT>& rhs)
        {
            const std::string separator = "\t";
            for (std::size_t y = 0; y < rhs.height; ++y)
            {
                for (std::size_t x = 0; x < rhs.width; ++x)
                {
                    //  Ref: https://isocpp.org/wiki/faq/input-output#print-char-or-ptr-as-number
                    os << +rhs.at(x, y) << separator;
                }
                os << "\n";
            }
            os << "\n";
            return os;
        }

        Image<ElementT>& operator+=(const Image<ElementT>& rhs)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::plus<>{});
            return *this;
        }

        Image<ElementT>& operator-=(const Image<ElementT>& rhs)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::minus<>{});
            return *this;
        }

        Image<ElementT>& operator*=(const Image<ElementT>& rhs)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::multiplies<>{});
            return *this;
        }

        Image<ElementT>& operator/=(const Image<ElementT>& rhs)
        {
            assert(rhs.width == this->width);
            assert(rhs.height == this->height);
            std::transform(std::ranges::cbegin(image_data), std::ranges::cend(image_data), std::ranges::cbegin(rhs.image_data),
                   std::ranges::begin(image_data), std::divides<>{});
            return *this;
        }

        friend bool operator==(Image<ElementT> const&, Image<ElementT> const&) = default;

        friend bool operator!=(Image<ElementT> const&, Image<ElementT> const&) = default;

        friend Image<ElementT> operator+(Image<ElementT> input1, const Image<ElementT>& input2)
        {
            return input1 += input2;
        }

        friend Image<ElementT> operator-(Image<ElementT> input1, const Image<ElementT>& input2)
        {
            return input1 -= input2;
        }

        Image<ElementT>& operator=(Image<ElementT> const& input) = default;  //  Copy Assign

        Image<ElementT>& operator=(Image<ElementT>&& other) = default;       //  Move Assign

        Image(const Image<ElementT> &input) = default;                       //  Copy Constructor

        Image(Image<ElementT> &&input) = default;                            //  Move Constructor
        
    private:
        std::size_t width;
        std::size_t height;
        std::vector<ElementT> image_data;

        void checkBoundary(const size_t x, const size_t y) const
        {
            assert(x < width);
            assert(y < height);
        }

    };

    template<arithmetic ElementT = double, arithmetic OutputT = ElementT>
    constexpr static Image<ElementT> dct3_detail(std::vector<Image<ElementT>> input, int plane_index)
    {
        std::size_t N3 = input.size();
        OutputT alpha1 = (plane_index == 0) ? (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))) : (static_cast<OutputT>(1.0));
        auto output = Image<OutputT>(input[plane_index].getWidth(), input[plane_index].getHeight());
        for (std::size_t y = 0; y < output.getHeight(); y++)
        {
            OutputT alpha2 = (y == 0) ? (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))) : (static_cast<OutputT>(1.0));
            for (std::size_t x = 0; x < output.getWidth(); x++)
            {
                OutputT sum{};
                OutputT alpha3 = (x == 0) ? (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))) : (static_cast<OutputT>(1.0));
                for (std::size_t inner_z = 0; inner_z < N3; inner_z++)
                {
                    auto plane = input[inner_z];
                    auto N1 = static_cast<OutputT>(plane.getWidth());
                    auto N2 = static_cast<OutputT>(plane.getHeight());
                    for (std::size_t inner_y = 0; inner_y < plane.getHeight(); inner_y++)
                    {
                        for (std::size_t inner_x = 0; inner_x < plane.getWidth(); inner_x++)
                        {
                            auto l1 = (std::numbers::pi / (2 * N1) * (2 * static_cast<OutputT>(inner_x) + 1) * x);
                            auto l2 = (std::numbers::pi / (2 * N2) * (2 * static_cast<OutputT>(inner_y) + 1) * y);
                            auto l3 = (std::numbers::pi / (2 * static_cast<OutputT>(N3)) * (2 * static_cast<OutputT>(inner_z) + 1) * static_cast<OutputT>(plane_index));
                            sum += static_cast<OutputT>(plane.at(inner_x, inner_y)) *
                                std::cos(l1) * std::cos(l2) * std::cos(l3);
                        }
                    }
                }
                auto N1 = static_cast<OutputT>(input[0].getWidth());
                auto N2 = static_cast<OutputT>(input[0].getHeight());
                output.at(x, y) = 8 * alpha1 * alpha2 * alpha3 * sum / (N1 * N2 * N3);
            }
        }
        return output;
    }

    template<arithmetic ElementT = double, arithmetic OutputT = ElementT>
    constexpr static std::vector<Image<ElementT>> dct3(std::vector<Image<ElementT>> input)
    {
        std::vector<Image<ElementT>> output;
        for (std::size_t i = 0; i < input.size(); i++)
        {
            output.push_back(dct3_detail(input, i));
        }
        return output;
    }
}

template<typename ElementT>
void print3(std::vector<TinyDIP::Image<ElementT>> input)
{
    for (std::size_t i = 0; i < input.size(); i++)
    {
        input[i].print();
        std::cout << "*******************\n";
    }
}

void dct3Test()
{
    std::size_t N1 = 10, N2 = 10, N3 = 10;
    std::vector<TinyDIP::Image<double>> test_input;
    for (std::size_t z = 0; z < N3; z++)
    {
        test_input.push_back(TinyDIP::Image<double>(N1, N2));
    }
    for (std::size_t z = 1; z <= N3; z++)
    {
        for (std::size_t y = 1; y <= N2; y++)
        {
            for (std::size_t x = 1; x <= N1; x++)
            {
                test_input[z - 1].at(y - 1, x - 1) = x * 100 + y * 10 + z;
            }
        }
    }
    print3(test_input);

    auto test_output = TinyDIP::dct3(test_input);
    print3(test_output);
}

int main()
{
    auto start = std::chrono::system_clock::now();
    dct3Test();
    auto end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end - start;
    std::time_t end_time = std::chrono::system_clock::to_time_t(end);
    std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
    return 0;
}

The output of the testing code above:

111 121 131 141 151 161 171 181 191 201 
211 221 231 241 251 261 271 281 291 301 
311 321 331 341 351 361 371 381 391 401 
411 421 431 441 451 461 471 481 491 501 
511 521 531 541 551 561 571 581 591 601 
611 621 631 641 651 661 671 681 691 701 
711 721 731 741 751 761 771 781 791 801 
811 821 831 841 851 861 871 881 891 901 
911 921 931 941 951 961 971 981 991 1001    
1011    1021    1031    1041    1051    1061    1071    1081    1091    1101    

*******************
112 122 132 142 152 162 172 182 192 202 
212 222 232 242 252 262 272 282 292 302 
312 322 332 342 352 362 372 382 392 402 
412 422 432 442 452 462 472 482 492 502 
512 522 532 542 552 562 572 582 592 602 
612 622 632 642 652 662 672 682 692 702 
712 722 732 742 752 762 772 782 792 802 
812 822 832 842 852 862 872 882 892 902 
912 922 932 942 952 962 972 982 992 1002    
1012    1022    1032    1042    1052    1062    1072    1082    1092    1102    

*******************
113 123 133 143 153 163 173 183 193 203 
213 223 233 243 253 263 273 283 293 303 
313 323 333 343 353 363 373 383 393 403 
413 423 433 443 453 463 473 483 493 503 
513 523 533 543 553 563 573 583 593 603 
613 623 633 643 653 663 673 683 693 703 
713 723 733 743 753 763 773 783 793 803 
813 823 833 843 853 863 873 883 893 903 
913 923 933 943 953 963 973 983 993 1003    
1013    1023    1033    1043    1053    1063    1073    1083    1093    1103    

*******************
114 124 134 144 154 164 174 184 194 204 
214 224 234 244 254 264 274 284 294 304 
314 324 334 344 354 364 374 384 394 404 
414 424 434 444 454 464 474 484 494 504 
514 524 534 544 554 564 574 584 594 604 
614 624 634 644 654 664 674 684 694 704 
714 724 734 744 754 764 774 784 794 804 
814 824 834 844 854 864 874 884 894 904 
914 924 934 944 954 964 974 984 994 1004    
1014    1024    1034    1044    1054    1064    1074    1084    1094    1104    

*******************
115 125 135 145 155 165 175 185 195 205 
215 225 235 245 255 265 275 285 295 305 
315 325 335 345 355 365 375 385 395 405 
415 425 435 445 455 465 475 485 495 505 
515 525 535 545 555 565 575 585 595 605 
615 625 635 645 655 665 675 685 695 705 
715 725 735 745 755 765 775 785 795 805 
815 825 835 845 855 865 875 885 895 905 
915 925 935 945 955 965 975 985 995 1005    
1015    1025    1035    1045    1055    1065    1075    1085    1095    1105    

*******************
116 126 136 146 156 166 176 186 196 206 
216 226 236 246 256 266 276 286 296 306 
316 326 336 346 356 366 376 386 396 406 
416 426 436 446 456 466 476 486 496 506 
516 526 536 546 556 566 576 586 596 606 
616 626 636 646 656 666 676 686 696 706 
716 726 736 746 756 766 776 786 796 806 
816 826 836 846 856 866 876 886 896 906 
916 926 936 946 956 966 976 986 996 1006    
1016    1026    1036    1046    1056    1066    1076    1086    1096    1106    

*******************
117 127 137 147 157 167 177 187 197 207 
217 227 237 247 257 267 277 287 297 307 
317 327 337 347 357 367 377 387 397 407 
417 427 437 447 457 467 477 487 497 507 
517 527 537 547 557 567 577 587 597 607 
617 627 637 647 657 667 677 687 697 707 
717 727 737 747 757 767 777 787 797 807 
817 827 837 847 857 867 877 887 897 907 
917 927 937 947 957 967 977 987 997 1007    
1017    1027    1037    1047    1057    1067    1077    1087    1097    1107    

*******************
118 128 138 148 158 168 178 188 198 208 
218 228 238 248 258 268 278 288 298 308 
318 328 338 348 358 368 378 388 398 408 
418 428 438 448 458 468 478 488 498 508 
518 528 538 548 558 568 578 588 598 608 
618 628 638 648 658 668 678 688 698 708 
718 728 738 748 758 768 778 788 798 808 
818 828 838 848 858 868 878 888 898 908 
918 928 938 948 958 968 978 988 998 1008    
1018    1028    1038    1048    1058    1068    1078    1088    1098    1108    

*******************
119 129 139 149 159 169 179 189 199 209 
219 229 239 249 259 269 279 289 299 309 
319 329 339 349 359 369 379 389 399 409 
419 429 439 449 459 469 479 489 499 509 
519 529 539 549 559 569 579 589 599 609 
619 629 639 649 659 669 679 689 699 709 
719 729 739 749 759 769 779 789 799 809 
819 829 839 849 859 869 879 889 899 909 
919 929 939 949 959 969 979 989 999 1009    
1019    1029    1039    1049    1059    1069    1079    1089    1099    1109    

*******************
120 130 140 150 160 170 180 190 200 210 
220 230 240 250 260 270 280 290 300 310 
320 330 340 350 360 370 380 390 400 410 
420 430 440 450 460 470 480 490 500 510 
520 530 540 550 560 570 580 590 600 610 
620 630 640 650 660 670 680 690 700 710 
720 730 740 750 760 770 780 790 800 810 
820 830 840 850 860 870 880 890 900 910 
920 930 940 950 960 970 980 990 1000    1010    
1020    1030    1040    1050    1060    1070    1080    1090    1100    1110    

*******************
1726.75 -80.7207    2.5284e-13  -8.64604    -7.77618e-14    -2.82843    1.59162e-13 -1.14371    -1.43473e-13    -0.320717   
-807.207    2.57244e-14 -1.24763e-13    -1.00325e-13    4.82332e-14 -7.91025e-14    -9.83958e-14    1.62707e-13 5.91661e-14 6.73658e-14 
3.81988e-13 -1.54346e-14    -1.28622e-14    -1.92933e-15    1.41484e-14 -3.85866e-15    -1.28622e-15    3.21555e-15 -5.46643e-15    8.84276e-15 
-86.4604    -1.22191e-14    -4.50177e-15    -8.36043e-15    5.14488e-15 1.92933e-15 7.07421e-15 7.71732e-15 2.57244e-15 -2.41166e-15    
-5.54792e-14    6.4311e-15  1.92933e-15 -1.28622e-15    -2.57244e-15    3.85866e-15 -4.50177e-15    7.71732e-15 -3.21555e-15    -2.21873e-14    
-28.2843    -5.14488e-15    1.28622e-15 -4.50177e-15    3.21555e-15 -5.78799e-15    -7.39576e-15    9.9682e-15  -1.28622e-15    9.64665e-15 
1.64164e-13 -7.71732e-15    4.50177e-15 2.57244e-14 -6.4311e-16 3.21555e-16 -1.28622e-15    3.85866e-15 -4.50177e-15    -3.85866e-15    
-11.4371    2.18657e-14 -1.57562e-14    -3.21555e-15    1.60777e-15 1.70424e-14 -3.21555e-16    -4.66255e-15    -1.60777e-15    -9.56626e-15    
-3.77668e-13    8.68198e-15 1.28622e-15 -5.78799e-15    -1.92933e-15    -1.447e-15  -5.14488e-15    2.73322e-15 2.81361e-15 1.00888e-14 
-3.20717    5.30566e-15 8.03887e-16 5.30566e-15 -1.68816e-14    -3.5371e-15 -1.63993e-14    7.39576e-15 1.2822e-14  1.18171e-14 

*******************
-8.07207    2.3152e-14  -7.71732e-15    -5.14488e-15    1.92933e-15 -6.4311e-16 -3.85866e-15    1.02898e-14 4.50177e-15 4.82332e-16 
-2.4824e-13 5.45697e-15 -5.45697e-15    1.18234e-14 -9.09495e-16    8.18545e-15 0   -4.54747e-15    4.09273e-15 -3.86535e-15    
-4.50177e-14    3.63798e-15 -1.81899e-15    6.36646e-15 1.81899e-15 6.36646e-15 3.63798e-15 -4.54747e-15    -4.54747e-15    -2.27374e-16    
6.4311e-15  -9.09495e-16    0   1.81899e-15 -9.09495e-16    0   -8.18545e-15    1.81899e-15 -4.54747e-15    4.09273e-15 
-1.54346e-14    -8.18545e-15    -1.81899e-15    -1.81899e-15    6.36646e-15 4.54747e-15 7.27596e-15 -1.81899e-15    -2.27374e-15    2.27374e-15 
-6.4311e-16 4.54747e-15 -3.63798e-15    -2.72848e-15    8.18545e-15 4.54747e-15 -5.91172e-15    -2.27374e-15    -6.13909e-15    -5.00222e-15    
1.02898e-14 6.36646e-15 1.81899e-15 -9.09495e-16    4.54747e-15 5.00222e-15 4.09273e-15 9.09495e-16 -2.27374e-16    8.6402e-15  
-1.31838e-14    -9.09495e-16    5.45697e-15 -1.36424e-15    -9.09495e-16    -2.72848e-15    -5.91172e-15    5.68434e-15 6.13909e-15 -8.18545e-15    
1.86502e-14 -5.00222e-15    -9.09495e-16    -4.54747e-16    -9.54969e-15    -4.77485e-15    -4.54747e-16    -9.09495e-16    5.00222e-15 3.41061e-15 
1.63993e-14 5.68434e-15 -6.82121e-15    -9.09495e-16    -9.77707e-15    -7.04858e-15    2.95586e-15 2.95586e-15 1.02318e-15 1.39266e-15 

*******************
3.31056e-13 -5.14488e-15    -1.41484e-14    -9.64665e-15    1.28622e-14 -6.4311e-15 6.4311e-15  3.21555e-15 -4.18021e-15    1.12544e-15 
-1.06756e-13    -1.09139e-14    -1.81899e-15    -9.09495e-16    -7.27596e-15    -9.09495e-16    -1.81899e-15    -4.54747e-15    4.54747e-15 -7.50333e-15    
-4.1159e-14 5.45697e-15 -3.63798e-15    -1.81899e-15    -5.45697e-15    -4.54747e-15    2.72848e-15 -2.27374e-15    9.09495e-16 -2.72848e-15    
-5.14488e-15    -1.45519e-14    -7.27596e-15    -9.09495e-16    -2.72848e-15    9.09495e-15 0   -2.72848e-15    -9.09495e-16    -1.11413e-14    
-3.21555e-15    3.63798e-15 1.00044e-14 -3.63798e-15    -4.54747e-15    -8.18545e-15    -4.54747e-15    -3.18323e-15    3.18323e-15 -5.22959e-15    
-1.41484e-14    6.36646e-15 0   9.09495e-16 -1.81899e-15    -8.18545e-15    -2.27374e-15    -7.7307e-15 1.59162e-15 -8.6402e-15 
1.92933e-14 1.81899e-15 -2.72848e-15    -1.81899e-15    2.27374e-15 -7.27596e-15    -2.72848e-15    2.27374e-15 1.81899e-15 3.18323e-15 
1.60777e-15 -9.09495e-16    4.54747e-15 2.27374e-15 0   5.91172e-15 -2.72848e-15    2.04636e-15 3.18323e-15 -5.91172e-15    
-1.54346e-14    3.18323e-15 -1.36424e-15    4.54747e-16 1.22782e-14 -4.3201e-15 8.41283e-15 0   -1.13687e-15    6.82121e-16 
-6.27032e-15    -4.54747e-16    6.82121e-16 -4.54747e-16    2.27374e-16 -2.38742e-15    -1.13687e-15    2.04636e-15 -6.82121e-16    -6.16751e-15    

*******************
-0.864604   -9.64665e-15    6.4311e-16  7.07421e-15 0   -7.07421e-15    1.35053e-14 -3.85866e-15    -6.4311e-16 -2.41166e-15    
-2.71392e-13    -8.18545e-15    -3.63798e-15    3.63798e-15 -8.18545e-15    -1.09139e-14    2.72848e-15 -3.63798e-15    -9.09495e-16    7.7307e-15  
-5.14488e-15    0   -9.09495e-15    -4.54747e-15    -2.72848e-15    -5.45697e-15    -2.72848e-15    -3.63798e-15    -6.36646e-15    -2.27374e-16    
-4.6947e-14 0   5.45697e-15 -4.54747e-15    -9.09495e-16    2.72848e-15 -4.54747e-15    -1.36424e-15    6.36646e-15 -1.18234e-14    
1.09329e-14 -8.18545e-15    5.45697e-15 -2.72848e-15    9.09495e-16 8.18545e-15 -1.09139e-14    3.18323e-15 2.27374e-15 5.45697e-15 
1.54346e-14 -1.81899e-15    -9.09495e-16    2.72848e-15 4.54747e-15 -2.72848e-15    -6.36646e-15    3.63798e-15 9.09495e-16 3.86535e-15 
5.78799e-15 -9.09495e-15    -7.27596e-15    -9.09495e-16    -3.63798e-15    4.54747e-15 1.22782e-14 -1.36424e-15    -2.27374e-16    5.68434e-16 
-1.28622e-14    5.45697e-15 5.91172e-15 6.82121e-15 -6.36646e-15    2.27374e-15 -2.27374e-15    1.36424e-15 3.86535e-15 2.6148e-15  
-5.46643e-15    2.72848e-15 -4.54747e-15    4.54747e-15 -2.04636e-15    -1.81899e-15    3.86535e-15 2.72848e-15 1.13687e-15 -1.47793e-15    
-1.78463e-14    8.6402e-15  -2.04636e-15    -2.27374e-15    1.28466e-14 2.38742e-15 -3.41061e-16    -1.02318e-15    -2.04636e-15    -9.37916e-16    

*******************
-4.09273e-14    1.1576e-14  1.09329e-14 3.85866e-15 -3.85866e-15    3.85866e-15 6.4311e-15  1.1576e-14  -5.14488e-15    -9.32509e-15    
-5.14488e-15    -6.36646e-15    -1.81899e-15    0   -4.54747e-15    -4.54747e-15    0   -5.45697e-15    -5.00222e-15    2.27374e-15 
-2.37951e-14    0   9.09495e-16 -3.63798e-15    -1.00044e-14    -9.09495e-16    3.63798e-15 5.91172e-15 5.00222e-15 -1.25056e-14    
6.4311e-16  9.09495e-16 7.27596e-15 -9.09495e-16    4.54747e-15 8.18545e-15 -1.63709e-14    5.00222e-15 -4.54747e-16    1.81899e-15 
-1.02898e-14    -4.54747e-15    -5.45697e-15    3.63798e-15 2.72848e-15 -1.81899e-15    1.09139e-14 -8.18545e-15    4.3201e-15  5.22959e-15 
-2.25088e-14    -1.81899e-15    -8.18545e-15    0   1.81899e-15 1.81899e-15 1.36424e-15 1.81899e-15 9.09495e-16 1.13687e-15 
1.80071e-14 -3.63798e-15    -5.00222e-15    9.09495e-16 5.45697e-15 6.82121e-15 -4.54747e-15    -9.54969e-15    3.86535e-15 -2.27374e-15    
-1.38269e-14    9.09495e-16 -3.63798e-15    -2.72848e-15    -1.81899e-15    2.72848e-15 -3.86535e-15    -1.81899e-15    -3.29692e-15    -7.56017e-15    
-2.95831e-14    5.45697e-15 -2.27374e-15    2.50111e-15 -1.36424e-15    1.11413e-14 9.54969e-15 5.34328e-15 -1.02318e-15    -1.04023e-14    
-1.09329e-14    -1.59162e-15    -9.09495e-16    -1.10276e-14    4.20641e-15 6.36646e-15 3.86535e-15 1.59162e-15 -7.04858e-15    1.29319e-14 

*******************
-0.282843   -9.00354e-15    -7.71732e-15    -3.21555e-15    9.64665e-15 -5.78799e-15    2.25088e-15 3.5371e-15  -1.92933e-15    4.50177e-15 
-1.08686e-13    2.72848e-15 -5.45697e-15    8.18545e-15 8.18545e-15 9.09495e-16 -4.54747e-16    3.18323e-15 -5.22959e-15    4.54747e-16 
4.75901e-14 -4.54747e-15    3.63798e-15 -9.09495e-16    0   -6.36646e-15    8.6402e-15  -1.81899e-15    -2.95586e-15    -8.6402e-15 
-3.08693e-14    -1.81899e-15    2.72848e-15 -9.09495e-16    1.00044e-14 0   6.36646e-15 -9.09495e-16    -1.36424e-15    1.29603e-14 
3.21555e-15 -1.81899e-15    2.72848e-15 1.81899e-15 3.63798e-15 0   2.27374e-15 -3.63798e-15    0   -6.13909e-15    
-3.47279e-14    -1.81899e-15    -5.45697e-15    0   9.09495e-16 5.00222e-15 4.54747e-16 1.13687e-15 -5.22959e-15    2.16005e-15 
-2.57244e-15    -2.27374e-15    5.00222e-15 5.45697e-15 2.27374e-15 7.27596e-15 9.09495e-16 -8.6402e-15 -9.09495e-16    -1.20508e-14    
-1.41484e-14    3.18323e-15 -5.91172e-15    4.54747e-16 8.18545e-15 -3.18323e-15    -2.27374e-16    -4.3201e-15 3.29692e-15 -1.81899e-15    
4.82332e-16 2.72848e-15 -1.13687e-14    5.22959e-15 6.13909e-15 4.09273e-15 5.00222e-15 -1.7053e-15 -2.38742e-15    1.47793e-15 
6.4311e-15  7.61702e-15 2.16005e-15 -4.54747e-16    2.27374e-16 -6.36646e-15    2.95586e-15 0   3.9222e-15  9.6918e-15  

*******************
1.64164e-13 -1.54346e-14    3.21555e-15 1.41484e-14 -3.21555e-15    5.46643e-15 -5.14488e-15    -1.28622e-15    -7.07421e-15    0   
-1.02898e-14    1.00044e-14 5.45697e-15 0   -7.27596e-15    -7.7307e-15 4.09273e-15 -4.54747e-15    2.95586e-15 1.59162e-14 
3.85866e-15 1.09139e-14 6.36646e-15 9.09495e-16 4.09273e-15 -3.63798e-15    9.09495e-16 4.54747e-16 9.09495e-16 1.36424e-15 
2.12226e-14 -3.63798e-15    1.81899e-15 -3.63798e-15    -1.00044e-14    -2.72848e-15    5.91172e-15 -3.63798e-15    5.68434e-15 1.87583e-14 
-4.88764e-14    3.63798e-15 -7.7307e-15 -5.45697e-15    -1.09139e-14    -9.54969e-15    -2.72848e-15    -7.7307e-15 9.32232e-15 -4.09273e-15    
-7.71732e-15    8.6402e-15  -9.54969e-15    -3.63798e-15    -8.6402e-15 -9.09495e-16    -9.09495e-16    -2.27374e-16    3.18323e-15 -4.77485e-15    
-2.57244e-15    5.91172e-15 -7.7307e-15 2.27374e-15 -5.91172e-15    5.91172e-15 -2.04636e-15    1.81899e-15 4.20641e-15 -5.85487e-15    
-2.89399e-15    -1.36424e-15    5.91172e-15 -7.7307e-15 7.50333e-15 -1.13687e-15    4.54747e-16 -1.81899e-15    6.36646e-15 -1.20508e-14    
-5.94877e-15    1.13687e-15 -6.82121e-16    9.54969e-15 2.72848e-15 5.91172e-15 -3.86535e-15    7.95808e-16 -2.04636e-15    -4.20641e-15    
-2.06599e-14    -6.48015e-15    4.3201e-15  3.97904e-15 2.72848e-15 -3.63798e-15    7.7307e-15  -5.68434e-17    2.70006e-15 6.96332e-15 

*******************
-0.114371   2.44382e-14 -3.21555e-16    -8.36043e-15    7.39576e-15 8.68198e-15 -4.50177e-15    -1.60777e-16    -1.92933e-15    -1.21387e-14    
2.08046e-13 -1.00044e-14    1.81899e-15 2.27374e-15 1.81899e-15 -4.54747e-15    -4.54747e-15    5.68434e-15 7.04858e-15 1.00044e-14 
-4.9841e-14 -2.72848e-15    9.09495e-15 9.09495e-16 -6.36646e-15    4.54747e-16 5.45697e-15 1.81899e-15 4.54747e-16 8.6402e-15  
3.85866e-14 5.45697e-15 -4.54747e-16    -4.09273e-15    -3.63798e-15    2.27374e-15 -5.00222e-15    -9.09495e-16    1.13687e-15 -1.37561e-14    
6.75265e-15 -2.72848e-15    -3.63798e-15    3.63798e-15 -9.09495e-16    -9.09495e-16    -2.04636e-15    1.13687e-15 2.6148e-15  1.53477e-15 
1.1576e-14  -4.09273e-15    -7.7307e-15 5.45697e-15 3.63798e-15 9.09495e-16 -2.95586e-15    2.50111e-15 -2.16005e-15    -1.00044e-14    
-1.06113e-14    5.91172e-15 6.82121e-15 -9.54969e-15    3.86535e-15 -1.13687e-15    9.09495e-15 -3.86535e-15    8.6402e-15  5.22959e-15 
-8.52121e-15    3.63798e-15 2.04636e-15 4.77485e-15 -3.18323e-15    1.59162e-15 -2.72848e-15    7.04858e-15 -5.22959e-15    1.33014e-14 
8.03887e-15 -5.22959e-15    5.00222e-15 -1.13687e-15    -2.72848e-15    1.19371e-14 -4.66116e-15    -7.95808e-16    7.21911e-15 -6.39488e-15    
-8.03887e-16    -8.07177e-15    2.27374e-15 5.68434e-16 -3.63798e-15    6.13909e-15 -9.09495e-16    8.69704e-15 -2.84217e-15    8.14282e-15 

*******************
-6.29825e-14    8.03887e-15 -8.36043e-15    -3.5371e-15 -6.4311e-15 4.34099e-15 -1.02898e-14    -6.4311e-16 -3.45672e-15    2.23079e-14 
6.75265e-14 7.7307e-15  -9.09495e-16    4.54747e-16 5.91172e-15 3.86535e-15 -8.6402e-15 6.82121e-15 -9.09495e-16    6.13909e-15 
-4.1159e-14 3.18323e-15 2.27374e-15 -2.27374e-15    9.54969e-15 -1.59162e-15    -4.3201e-15 2.50111e-15 1.36424e-15 5.22959e-15 
2.25088e-15 -1.00044e-14    -4.54747e-15    2.72848e-15 5.22959e-15 -7.27596e-15    9.32232e-15 -1.13687e-15    4.54747e-16 -1.47793e-15    
-1.67209e-14    -2.72848e-15    1.81899e-15 2.50111e-15 6.82121e-15 5.68434e-15 -1.59162e-15    1.93268e-15 -2.16005e-15    -1.67688e-14    
-4.66255e-15    6.36646e-15 4.54747e-15 1.59162e-15 7.95808e-15 -4.54747e-15    4.54747e-16 -2.50111e-15    -3.41061e-16    -9.43601e-15    
-7.23499e-15    7.04858e-15 9.32232e-15 -4.54747e-16    6.36646e-15 -1.04592e-14    7.50333e-15 8.07177e-15 -6.59384e-15    1.25056e-15 
2.47597e-14 1.38698e-14 -5.45697e-15    2.50111e-15 -6.82121e-16    1.28466e-14 -1.47793e-15    -3.75167e-15    3.12639e-15 -1.84741e-15    
1.96952e-14 -5.45697e-15    -3.75167e-15    -1.09139e-14    -3.63798e-15    -3.18323e-15    -1.00613e-14    -4.77485e-15    2.04636e-15 1.6172e-14  
9.36529e-15 9.09495e-16 3.52429e-15 4.66116e-15 -1.3415e-14 -1.02318e-14    -3.78009e-15    1.90425e-15 7.24754e-15 -1.09424e-14    

*******************
-0.0320717  1.43092e-14 4.34099e-15 6.27032e-15 -9.80743e-15    8.36043e-15 -5.94877e-15    -4.90371e-15    1.25004e-14 1.95345e-14 
5.75583e-14 -7.04858e-15    4.54747e-16 4.54747e-15 3.86535e-15 -1.13687e-14    6.59384e-15 2.95586e-15 -9.09495e-16    4.83169e-16 
1.94541e-14 -1.81899e-15    -5.22959e-15    0   3.41061e-15 -3.97904e-15    -3.18323e-15    -1.59162e-15    4.54747e-16 -6.62226e-15    
-1.12544e-15    -2.27374e-15    -2.95586e-15    5.00222e-15 3.97904e-15 -1.02318e-15    1.7053e-15  -4.20641e-15    -1.36424e-15    3.60956e-15 
1.92933e-15 -1.13687e-15    9.09495e-16 -1.93268e-15    9.89075e-15 -1.31877e-14    1.13687e-15 6.82121e-15 -7.50333e-15    -1.02602e-14    
-9.00354e-15    4.20641e-15 8.98126e-15 -4.09273e-15    -1.13687e-15    2.04636e-15 1.81899e-15 4.54747e-15 -1.53477e-15    8.78231e-15 
-1.80875e-14    -4.20641e-15    5.91172e-15 6.25278e-15 8.18545e-15 6.82121e-16 3.97904e-15 -2.84217e-16    -8.78231e-15    4.23483e-15 
-1.17368e-14    4.66116e-15 5.45697e-15 -6.0254e-15 -1.36424e-15    2.50111e-15 6.82121e-16 5.74119e-15 1.93268e-15 4.50484e-15 
4.22041e-15 1.93268e-15 -6.82121e-16    1.93268e-15 1.13687e-16 4.54747e-16 1.98952e-16 -5.3717e-15 1.15676e-14 -2.07478e-15    
9.305e-15   -1.19371e-15    4.66116e-15 4.74643e-15 3.32534e-15 3.35376e-15 6.79279e-15 3.90799e-15 -1.8332e-15 -3.41061e-16    

*******************
Computation finished at Wed Dec 29 00:34:06 2021
elapsed time: 0.0305631

A Godbolt link is here.

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

\$\endgroup\$

1 Answer 1

5
\$\begingroup\$

Here is what I noticed on a first viewing (of only the code of dct3 and dct3_detail, not the test code really). I haven't thought really about correctness or performance, just the things I noticed reading the code:

  1. dct3 and dct3_detail should not be static. Without very specific reasons, such as requiring different function pointers between translation units, there is no benefit to declaring a function template with internal linkage (i.e. static). It may even make the linker's work harder, since it forces it to consider instantiations with the same template arguments in different translation units as different functions.

  2. dct3 and dct3_detail should take the input parameter with type const std::vector<Image<ElementT>>& since they do not modify the input. With the current non-reference type, each call to these functions will make a complete copy of the input.

  3. In dct3: It is known up-front how large the output vector will be. Sufficient memory for that can be reserved up-front to avoid the reallocations which it would otherwise perform: output.reserve(input.size());

  4. OutputT is not used at all. dct3_detail is called without specifying it and because it doesn't appear in its parameters, the call dct3_detail(input, i) will always use the default OutputT = ElementT. If this is unintended, you should explicitly specify the template argument in the call: dct3_detail<ElementT, OutputT>(input, i) Also note that you could simplify this to dct3_detail<OutputT>(input, i) if you swapped the template parameters of dct3_detail. I suppose in dct3 it should then also be std::vector<Image<OutputT>> output; and the two return types should also be fixed accordingly.

  5. The type of plane_index should be std::size_t. At least I see no reason for it to be anything else.

  6. In std::size_t N3 = input.size(); you could also write auto N3 = input.size();. You seem to use the auto _ = declaration form everywhere else, so that would be more consistent. The same for OutputT alpha1 = .

  7. It seems very unlikely to me that any arithmetic type makes sense for OutputT. It seems to me, that OutputT should be constrained to floating type, i.e. std::floating_point or at least a non-floating-point type shouldn't be used in the intermediate steps of the calculation.

  8. (static_cast<OutputT>(1.0) / static_cast<OutputT>(std::sqrt(2))): C++ has a constant defined for sqrt(2) with highest possible precision for every floating point type and 1/sqrt(2) is just sqrt(2)/2, so assuming you followed point 7, you can use instead std::numbers::sqrt2_v<OutputT> / 2.

  9. static_cast<OutputT>(1.0) can be written simpler as OutputT{1.0} with the added benefit that this form is required to produce a diagnostic if the value in the braces cannot be represented in the OutputT type.

  10. It seems that you don't really access input in dct3_detail except as input[plane_index]. The only exception is input[0].getWidth(), but it appears that the result is required to be independent of the index anyway, right? If so, don't pass input and plane_index to dct3_detail. Just pass the element you currently want to use as a reference:

    constexpr Image<OutputT> dct3_detail(const Image<ElementT>& input_element)
    {
        // input[plane_index] -> input_element
        // input[0] -> input_element
    }
    
    //...
    
    for (auto& input_element : input)
    {
        output.push_back(dct3_detail(input_element));
    }
    
  11. std::numbers::pi should probably be std::numbers::pi_v<OutputT> to match the precision of the requested output type.

  12. In the test code the Image class is not constexpr friendly. It can never be used in a constexpr function. For templates this is not strictly checked, but because every specialization of dct3_detail and dct3 use Image, no specialization satisfies the constexpr conditions. Marking a function template with this property as constexpr anyway technically makes the program ill-formed, no diagnostic required. I suggest you either define Image in such a way that it can be used in constant expressions or you remove constexpr from the functions.

  13. std::cos is also not constexpr and for similar reasons as above, as long as you use it, you should not mark the functions constexpr.

  14. You are getting N1 and N2 twice, although my interpretation is that their values are fixed. You also seem to get them without type conversion multiple times in the loops. I would just get them once at the beginning of dct3_detail, once in unconverted form an once in converted form. It might be a performance benefit to get them even only once in the whole calculation.

\$\endgroup\$
4
  • 1
    \$\begingroup\$ About restricting OutputT to only floating point types: it was quite common to implement DCT algorithms using integers, as this was faster on a lot of hardware, sometimes even those with FPUs, before vector instructions became mainstream. But this requires carefully avoiding intermediate results smaller than 1 for example, or storing those intermediates and doing calculations with them as float or double, which defeats the purpose. \$\endgroup\$
    – G. Sliepen
    Dec 29, 2021 at 15:15
  • \$\begingroup\$ @G.Sliepen It might make sense to use three parameters, one for input type, one for ouput type used in the Image<...> return types and one intermediate type in which calculations are performed, which would replace all the other uses of OutputT in the current code. Then only the intermediate type really needs a constraint on floating point types. As you are saying, doing in the calculations in an integer type will require some major modifications and additional considerations. \$\endgroup\$ Dec 29, 2021 at 15:30
  • \$\begingroup\$ @G.Sliepen Thank you for the comments. As you mentioned, "it was quite common to implement DCT algorithms using integers". Yeah, but this still depends on the application / processing workflow, right? \$\endgroup\$
    – JimmyHu
    Dec 29, 2021 at 15:45
  • 1
    \$\begingroup\$ @JimmyHu Sure. For audio and video encoding/decoding, integer DCT is a useful technique (especially since the original values are all in the range 0..255). If you need more precision and are not doing real-time data processing, float and double would be used. \$\endgroup\$
    – G. Sliepen
    Dec 29, 2021 at 18:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.