Issues with "zeroing out" DCT coefficients

28 views (last 30 days)
Ryder
Ryder on 11 Dec 2024 at 3:13
Commented: Walter Roberson on 11 Dec 2024 at 4:31
Hello Matlab!
I'm a bit new to Matlab but I'm working on a really important project and so I'm really desperate, which is why I'm posting this question.
I tried to make this post as concise as possible, but its still quite extensive, which is why I've tried to organize it the best I can. My key questions are bolded if you don't want to read the rest.
Context:
I have recently been working on a project which involves the JPEG algorithm. My goal is to create a program that takes an image and an input of desired file size, and compresses the image to that size. The two variables I'm looking at are: 'k', which represents a multiplier for the default quantization table; and 't', which represents the tolerance for which coefficients are set to zero after quantization (typically this value is 0.5 as the coefficients are rounded to the nearest integer).
I created a function 'g' to estimate the file size. I used a database to calculated an approximation for "average bits per character", which I then multiply by the number of pixels ('n') and the probability that a pixel is outside the tolerance of T (ie 1 - P(-T<x<T) ). This probability was calculated by integrating a laplacian distribution and using the identity b = √2/σ.
The functions are related by a lagrangian, which I use to minimize quality loss at the desired file size.
I put my code at the bottom if anyone wanted to try running it (you would need to have an image path for image.jpg).
Issue:
My issue is currently with the final step of using an inverse DCT to view the compressed image. What I was trying to do is divide each coefficient by its corresponding value in the quantization table, and then set it equal to 0 if the value is less than 'T'. How can I write this in code? Assuming I multiplied the array 'Q' by 'k', naming the new array 'Qk' I was thinking something like this ('C' is an array containing the dct coefficients after using dct() on the image matrix):
C(abs(C)./Qk<T) = 0;
But this doesn't work. Even something like:
C(abs(C)<0.2) = 0;
Doesn't actually change any values in 'C' (I checked manually), even though I saw something similar in a video. However, running:
C(C<0.2) = 0;
does actually work as intended, but doesn't work for my scenario as many values in 'C' are negative. So I'd appreciate if someone could explain how to set the the values below a certain magnitude to 0.
However, I noticed another glaring issue. When I was going through the DCT coefficients of a test image, I noticed that they all had a smaller magnitude than I expected (the magnitude of most of them is in between 0 and 0.5). The value of 'k' is 0.0644 no matter what size I put in. The value of 't' is -0.3500 when the input size is 4e7, but it varies as I change the input; the value of 'l' (lagrange multiplier) is always 0. Although 't' being negative shouldn't be an issue (can be fixed with an abs()), I'm concerned about the other values of 't' at different sizes. For example, inputting 3e7 yields 't = 8.35'. My issue is that this value is not realistically attainable by dividing a DCT coefficient (who's magnitude is typically in between 0 and 0.5) by k (0.0644) times Q (a value in the quantiztation table which can't be lower than 10). I'm pretty bad with words here so in code that would be:
C(i)/(k*Q(i))
Does this mean it's impossible to achieve this size without setting all coefficients to 0? Relatively speaking, a size of 4e7 to 3e7 shouldn't mean that the image should go completely black right? I was under the impression that most of the high frequency data (thats practically invisible to the eye) would be erased, but the lower frequency data would still remain, but even the smallest quantization divisor (being 10) doesn't result in a quotient that is anywhere near clearing the threshold of 't' for a non-zero coefficient (in this case).
Any help is enoumously appreciated, even if its just partly answering one of the questions, particularly the first one.
THANK YOU SO MUCH MATLAB COMMUNITY!
Code:
syms k t l
baka = imread('image.jpg');
img = im2gray(imread('image.jpg'));
C = dct2(double(img));
I = C';
I = I(:)'; %turns the matrix into an array for easier computation
Q = [16 11 10 16 24 40 51 61 12 12 14 19 26 58 60 55 14 13 16 24 40 57 69 56 14 17 22 29 51 87 80 62 18 22 37 56 68 109 103 77 24 35 55 64 81 104 113 92 49 64 78 87 103 121 120 101 72 92 95 98 112 100 103 99];
n=length(I)/length(Q);
Q = repmat(Q,1,n); %concatenate Q to be the same length as I
IQ = [I./Q];
n = length(I);
b = std(IQ);
b = sqrt(2)*k/sqrt(b);
f = (sum(I.^2) -2/k*sum((I.^2)./Q) + sum((I./Q).^2)/k^2)/n %function for quality loss (based on mean squared error)
g = (3.295)*n*(2.71^(-1*t*b))-40000000; %function for size (calculated 3.295 for the "average bits per coefficient stored")
F = f - l*g;
f0 = diff(F, k);
f1 = diff(F, t);
f2 = diff(F, l);
equations = @(k,t,l) [f0;f1;f2];
[sol_k, sol_t, sol_l] = solve([f0;f1;f2], [k,t,l]);
sols = [sol_k, sol_t, sol_l]
K = eval(sols(1))
T = eval(sols(2))
L = eval(sols(3))
%C(abs(C)./Qk<T) = 0; this doesn't work
%C(C<0.2) = 0; this works but not what I intended
Ct = idct2(C);
%figure
cla reset
imshowpair(img, [Ct, []], 'montage') %display before and after next to eachother
  4 Comments
Walter Roberson
Walter Roberson on 11 Dec 2024 at 4:29
Mathworks does not document any definition of the result of eval() of a symbolic expression.
In practice, it is equivalent to
eval(char(sols(1)))
However, char() of a symbolic expression gives a result that is not exactly MATLAB expression and is not internal MuPAD code. char() of symbolic expressions generates something that is intended to be human readable instead of being mechanically processed. char() of a symbolic expression produces some expressions that have parameters in a different order than the corresponding MATLAB code, and char() of a symbolic expression produces expressions that do not have valid syntax for MATLAB or for MuPAD.
Therefore you should never eval() a symbolic expression. You should subs() as necessary, and then finally double() the result.
Walter Roberson
Walter Roberson on 11 Dec 2024 at 4:31
C(abs(C)./Qk<T) = 0;
Qk is not defined in your code.

Sign in to comment.

Answers (0)

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!