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}) = \sum_{{k_1 = 0}}^{N_1 - 1} \sum_{{k_2 = 0}}^{N_2 - 1} \sum_{{k_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:

Another approach - with cellfun function

Besides using nested for for IDCT calculation, there is another way to simplify with cellfun function:

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

%%%    construct input blocks
InputCollection = cell(sizex, sizey, sizez);
for z = 1:sizez
for y = 1:sizey
for x = 1:sizex
Input = zeros(sizex, sizey, sizez);
Input(x, y, z) = 1;
InputCollection{x, y, z} = Input;
end
end
end

%%% calculation
OutputCollection = cell(sizex, sizey, sizez);
OutputCollection = cellfun(@IDCT3D, InputCollection, 'UniformOutput', false);
OutputCollection = cellfun(@myoffset, OutputCollection, 'UniformOutput', false);
save("OutputCollection_" + sizex + "_" + sizey + "_" + sizez + ".mat", "OutputCollection");

%%% plot
%load("OutputCollection_" + sizex + "_" + sizey + "_" + sizez + ".mat");
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}, alpha, title_text, xlabel_text, ylabel_text, zlabel_text);
saveas(gcf, OutputFilename);
end
end
end

function [output_array] = myoffset(input_array)
output_array = (input_array + 1) ./ 2;         %%  Adjust the range of numbers
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 syntax OutputCollection = cellfun(@IDCT3D, InputCollection, 'UniformOutput', false) used above is for calculating each cell in InputCollection with IDCT3D function and return the result to OutputCollection. However, as far as I know, there is one disadvantage which is that the intermediate status of calculation cannot be passed out. In other words, not like the for loop version which the saving operation can be inserted in to save the intermediate status, the object OutputCollection isn't updated until cellfun finish all the calculation. Based on this reason, I prefer to use the for version than the cellfun version. I am looking forward to some feedback, such as any other pros and cons, for the for loop version and the cellfun version.

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.

• In the equation at the top of the post, you probably meant to sum over k, not x. You also have two copies of the plot3dcube function. I maybe have time tonight to write a review, I see lots of opportunities for improvement, writing a review will take time! :) Commented Sep 4, 2021 at 15:50
• @Cris Luengo Thank you for letting me know the problem. The formula has been updated and it might be correct. If there is any problem, please let me know. Commented Sep 5, 2021 at 0:09

Your DCT computation can be made significantly more efficient, though for these small array sizes it doesn't matter a whole lot. I'll tackle this in the 2nd part of this answer. First I'll discuss some MATLAB syntax and style improvements.

OutputCollection{x, y, z} = {Output};


This puts Output into a 1x1 cell array, and assigns this cell array into the larger cell array. You needed to retrieve your array later with OutputCollection{x, y, z}{1} because of this. The {1} is surprising, and useless. Instead, assign with either of these two syntaxes:

OutputCollection{x, y, z} = Output;
OutputCollection(x, y, z) = {Output};


Either one of them assigns the array directly into one of the cells of the cell array.

function [] = plot3dcube(...)


It is not necessary to give []= here:

function plot3dcube(...)

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]);


Here you assigned the same value to three different variables. Simpler would be:

cubeColor = [1 1 1] * input_array(x, y, z);
plotcube([1 1 1], [x - 1 y - 1 z - 1], alpha, cubeColor);


The MATLAB Editor highlights this line:

IDCT3DOutput(n1+1,n2+1,n3+1)=sm;


and notes that the variable changes size every loop iteration. This is highly inefficient. You should always preallocate. Simply adding a line

IDCT3DOutput = zeros(N1,N2,N3);


at the top of the function, before the loops, significantly speeds up execution.

Your loops are sorted in reverse order. When looping over an array, indexed with X(k1,k2,k3), the loop over k1 should be the inner loop, and the one over k3 should be the outer loop. This causes the code to access the values in the array X in memory order, which is much faster because the cache will be used better. Likewise for the other three loops that iterate over the output array IDCT3DOutput. A 8x8x8 array is likely not large enough for this to matter a whole lot, but doing the same for a larger array that doesn't fit in the cache will make a big difference.

parfor should always be the outer loop. There is an overhead in distributing the work over the workers (and hopefully you have these configured as threaded, rather than separate processes as is the default), and doing so once or N1 * N2 * N3 times makes a big difference.

N1=size(X,1);
N2=size(X,2);
N3=size(X,3);


can be more easily written as

[N1,N2,N3] = size(X);


Finally, try to be consistent with spacing around operators. Some assignments you don't use any spacing around the operator, and some you have spaces on both sides. I'd suggest to pick one style, and use it consistently. It's worst when there's spaces on one side of an operator only:

sm=sm+ alpha1 ...


is not as easy to read as

sm = sm + alpha1 ...


In this second part I want to discuss the implementation of the algorithm itself, which can be much more efficient. Let's look at the three inner loops only:

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;


You are doing a terrible amount of useless work. For example, you test if(k1==0) N2*N3 times too many. For the k2 and k3 loops, this is a constant, so move the computation out of those loops. The same is true for the cos computations, which are quite expensive. This could look like this (additionally I've reversed the loop order as discussed earlier):

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


In my test this is about a 50% reduction in time. But we now sill compute alpha1 N3*N2 times too many, to avoid doing that we should compute all possible alpha1 values before the loops:

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


This yields another 3.5x speed improvement. Code now also has become simpler: we can do natural loops k1=1:N1, and index naturally X(k1,k2,k3). And it's now also much simpler to see how we can vectorize the computation. We use MATLAB's implicit singleton expansion here. This requires giving the arrays alpha1, alpha2 and alpha3 the right orientation:

alpha1 = cos(pi/(2*N1)*(2*n1+1)*(0:N1-1).'); % along 1st dimension
alpha1(1) = alpha1(1) / sqrt(2);
alpha2 = cos(pi/(2*N2)*(2*n2+1)*(0:N2-1)); % along 2nd dimension
alpha2(1) = alpha2(1) / sqrt(2);
alpha3 = cos(pi/(2*N3)*(2*n3+1)*(0:N3-1));
alpha3(1) = alpha3(1) / sqrt(2);
alpha3 = reshape(alpha3,1,1,N3); % along 3rd dimension
sm = ((X .* alpha1) .* alpha2) .* alpha3;
IDCT3DOutput(n1+1,n2+1,n3+1) = sum(sm(:));


This yields yet again a 3.5x speed improvement, for a total of about 14x over the original code. And yet again, the code was simplified.

When we now look at the outer loops, we again see that we can move some of these computations out of loops, which again significantly reduces the number of computations. But let's take a step back first. The DCT equation, $$\begin{split} x(n_{1}, n_{2}, n_{3}) = &\sum_{{k_1 = 0}}^{N_1 - 1} \sum_{{k_2 = 0}}^{N_2 - 1} \sum_{{k_3 = 0}}^{N_3 - 1} \epsilon_{k_{1}} \epsilon_{k_{2}} \epsilon_{k_{3}} X(k_{1}, k_{2}, k_{3}) \\ & \cos({\frac {\pi}{2N_{1}} (2n_{1} + 1)k_{1}}) \\ & \cos({\frac {\pi}{2N_{2}} (2n_{2} + 1)k_{2}}) \\ & \cos({\frac {\pi}{2N_{3}} (2n_{3} + 1)k_{3}}) \end{split}$$ can be rewritten to separate out summations, $$\begin{split} x(n_{1}, n_{2}, n_{3}) = & \sum_{{k_1 = 0}}^{N_1 - 1} \epsilon_{k_{1}} \cos({\frac {\pi}{2N_{1}} (2n_{1} + 1)k_{1}}) \Big[ \\ & \sum_{{k_2 = 0}}^{N_2 - 1} \epsilon_{k_{2}} \cos({\frac {\pi}{2N_{2}} (2n_{2} + 1)k_{2}}) \Big[ \\ & \sum_{{k_3 = 0}}^{N_3 - 1} \epsilon_{k_{3}} \cos({\frac {\pi}{2N_{3}} (2n_{3} + 1)k_{3}}) \Big[ X(k_{1}, k_{2}, k_{3}) \Big] \Big] \Big] \end{split}$$ I split the lines like this and added the big brackets to illustrate that we are applying the 1D DCT three times, on each of the axes of the image. We are basically doing:

out = DCT1D(in, 1);
out = DCT1D(out, 2);
out = DCT1D(out, 3);


(with the 2nd argument the axis over which it is applied.

Let's implement this:

function out = IDCT3D(X)
out = IDCT1D(X, 1);
out = IDCT1D(out, 2);
out = IDCT1D(out, 3);
end

function out = IDCT1D(X, dim)
N = size(X, dim);
sz = ones(1,3);
sz(dim) = N;
index = {':',':',':'};
out = zeros(size(X));
for n = 0:N-1
alpha = cos(pi/(2*N)*(2*n+1)*(0:N-1));
alpha(1) = alpha(1) / sqrt(2);
alpha = reshape(alpha, sz);
index{dim} = n + 1;
out(index{:}) = sum(X .* alpha, dim);
end
end


This gives me about 80x speedup over the previous fastest implementation, about 1200x speedup over the original one. (I hope I didn't mess up anything and the result is still correct...)