algorithm – Three-dimensional DCT spatial frequency components illustration

I am attempting to make three-dimensional discrete cosine transformation spatial frequency components illustration. The 3D cubes is used to present the level of each coefficient. The more bright the cube, the higher value of its coefficient, and vice versa. The 3D inverse discrete cosine transformation is used here and the formula is given as below.

The 3D inverse discrete cosine transformation $$x(n_1, n_{2}, n_{3})$$ of size $$N_{1} times N_{2} times N_{3}$$ is

$$x(n_{1}, n_{2}, n_{3}) = {frac {8}{N_{1} N_{2} N_{3}}} sum_{{n_1 = 0}}^{N_1 – 1} sum_{{n_2 = 0}}^{N_2 – 1} sum_{{n_3 = 0}}^{N_3 – 1} epsilon_{k_{1}} epsilon_{k_{2}} epsilon_{k_{3}} X(k_{1}, k_{2}, k_{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}})
$$

where

$$
n_{1} = 0, 1, dots, N_{1} – 1
$$

$$
n_{2} = 0, 1, dots, N_{2} – 1
$$

$$
n_{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
$$

The experimental implementation

  • plot3dcube function implementation: The function plotcube is from here.

    function () = plot3dcube(input_array, alpha, title_text, xlabel_text, ylabel_text, zlabel_text)
        for z = 1:size(input_array, 3)
            for y = 1:size(input_array, 2)
                for x = 1:size(input_array, 1)
                    cubeColorR = input_array(x, y, z);
                    cubeColorG = input_array(x, y, z);
                    cubeColorB = input_array(x, y, z);
                    plotcube((1 1 1), (x - 1 y - 1 z - 1), alpha, (cubeColorR cubeColorG cubeColorB));
                end
            end
        end
        title(title_text);
        xlabel(xlabel_text);
        ylabel(ylabel_text);
        zlabel(zlabel_text);
    end
    
  • IDCT3D.m: The implementation of inverse DCT calculation.

    function IDCT3DOutput=IDCT3D(X)
    
    N1=size(X,1);
    N2=size(X,2);
    N3=size(X,3);
    
    for n1=0:N1-1
        for n2=0:N2-1
            for n3=0:N3-1
                sm=0;
                parfor k1=0:N1-1
                    for k2=0:N2-1
                        for k3=0:N3-1
                            if(k1==0)
                                alpha1 = 1/sqrt(2);
                            else
                                alpha1 = 1;
                            end
                            if(k2==0)
                                alpha2 = 1/sqrt(2);
                            else
                                alpha2 = 1;
                            end
                            if(k3==0)
                                alpha3 = 1/sqrt(2);
                            else
                                alpha3 = 1;
                            end
                            sm=sm+ alpha1*alpha2*alpha3*X(k1+1,k2+1,k3+1)*...  
                                cos(pi/(2*N1)*(2*n1+1)*k1)*cos(pi/(2*N2)*(2*n2+1)*k2)*cos(pi/(2*N3)*(2*n3+1)*k3);                        
                        end
                    end
                end
                IDCT3DOutput(n1+1,n2+1,n3+1)=sm;
            end
        end
    end
    

The main testing code:

sizex = 8;
sizey = 8;
sizez = 8;
OutputFigureLocation = "./Figures/";

OutputCollection = cell(sizex, sizey, sizez);

%%% calculation
for z = 1:sizez
    for y = 1:sizey
        for x = 1:sizex
            Input = zeros(sizex, sizey, sizez);
            Input(x, y, z) = 1;
            Output = IDCT3D(Input);
            Output = (Output + 1) ./ 2;         %%  Adjust the range of numbers
            OutputCollection{x, y, z} = {Output};
            fprintf('%d_%d_%d size 3D IDCT: The %d_%d_%d / %d_%d_%d block calculation has done. n' , sizex, sizey, sizez, x, y, z, sizex, sizey, sizez);
        end
    end
end

%%% plot
alpha = 0.3;
title_text = "Three-dimensional DCT spatial frequency components";
xlabel_text = 'x';
ylabel_text = 'y';
zlabel_text = 'z';
for z = 1:sizez
    for y = 1:sizey
        for x = 1:sizex
            close all;
            OutputFilename = OutputFigureLocation + x + '_' + y + '_' + z + '.png';
            if isfile(OutputFilename)
                continue;
            end
            plot3dcube(OutputCollection{x, y, z}{1}, alpha, title_text, xlabel_text, ylabel_text, zlabel_text);
            saveas(gcf, OutputFilename);
        end
    end
end


function () = plot3dcube(input_array, alpha, title_text, xlabel_text, ylabel_text, zlabel_text)
    for z = 1:size(input_array, 3)
        for y = 1:size(input_array, 2)
            for x = 1:size(input_array, 1)
                cubeColorR = input_array(x, y, z);
                cubeColorG = input_array(x, y, z);
                cubeColorB = input_array(x, y, z);
                plotcube((1 1 1), (x - 1 y - 1 z - 1), alpha, (cubeColorR cubeColorG cubeColorB));
            end
        end
    end
    title(title_text);
    xlabel(xlabel_text);
    ylabel(ylabel_text);
    zlabel(zlabel_text);
end

The several output examples:

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

Reference

  • Wikipedia – Discrete cosine transform
  • Malavika Bhaskaranand and Jerry D. Gibson, “Distributions of 3D DCT coefficients for video,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, 2009.
  • J. Augustin Jacob and N. Senthil Kumar, “Determining Optimal Cube for 3D-DCT Based Video Compression for Different Motion Levels,” ICTACT Journal on Image and Video Processing, Vol. 03, November 2012.

If there is any possible improvement, please let me know.