Artifacts after DCT Transformation
I am getting artifacts after applying a DCT transformation. Here is a link with screenshots showing the expected output vs. the artifacts that appear.
Expected Fame:
:
Actual Frame:
:
I believe the problem lies in the matrix multiplication step inside either dct.hpp or quantize.hpp
What I have done to narrow down the problem: Compressor: Before writing into the bitstream (after applying DCT->Quantization) I have asserted that all values range between -128 and 127.
Decompressor: when reading from the bitstream (before applying the inverse of Quantization -> DCT) all values are static_casted to ensure negative values are correctly read. I also assert that all values read from the bistream lie between -128 and 127 before applying the inverse operations.
Additionally, blocks emitted from the decompressor match the blocks read by the compressor.
I have verified that commenting out the forward transformations (compressor.cpp) and the inverse transformations (decompressor.cpp) produces the correct output.
dct.hpp:
//dct.hpp
static const Matrix8d C {
{0.3535533905932738,0.3535533905932738,0.3535533905932738,0.3535533905932738,0.3535533905932738,0.3535533905932738,0.3535533905932738,0.3535533905932738,},
{0.4903926402016152,0.41573480615127256,0.2777851165098011,0.0975451610080641,-0.0975451610080641,-0.2777851165098011,-0.41573480615127256,-0.4903926402016152,},
{0.46193976625564337,0.19134171618254492,-0.19134171618254492,-0.46193976625564337,-0.46193976625564337,-0.19134171618254492,0.19134171618254492,0.46193976625564337,},
{0.4157348061512727,-0.09754516100806418,-0.4903926402016151,-0.2777851165098011,0.2777851165098011,0.4903926402016151,0.09754516100806418,-0.4157348061512727,},
{0.35355339059327373,-0.35355339059327373,-0.35355339059327373,0.35355339059327373,0.35355339059327373,-0.35355339059327373,-0.35355339059327373,0.35355339059327373,},
{0.2777851165098011,-0.4903926402016151,0.09754516100806418,0.4157348061512727,-0.4157348061512727,-0.09754516100806418,0.4903926402016151,-0.2777851165098011,},
{0.19134171618254492,-0.46193976625564337,0.46193976625564337,-0.19134171618254492,-0.19134171618254492,0.46193976625564337,-0.46193976625564337,0.19134171618254492,},
{0.0975451610080641,-0.2777851165098011,0.41573480615127256,-0.4903926402016152,0.4903926402016152,-0.41573480615127256,0.2777851165098011,-0.0975451610080641,},
};
class DCT {
private:
inline static Matrix8d tmp{8, 8};
public:
static Matrix8d forward(Matrix8d A) {
A = C * A * C.transpose();
return A;
}
static Matrix8d inverse(Matrix8d D) {
D = C.transpose() * D * C;
D = D.array().round();
return D;
}
}
quantize.hpp:
//quantize.hpp
static const Matrix8d QYb { // quantization matrix for Yb plane
{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},
{25, 35, 55, 64, 81, 104, 113, 92},
{49, 64, 78, 87, 103, 121, 120, 101},
{72, 92, 95, 98, 112, 100, 103, 99},
};
static const Matrix8d QCb { //quantization matrix for Cr and Cb planes
{17, 18, 24, 47, 99, 99, 99, 99},
{18, 21, 26, 66, 99, 99, 99, 99},
{24, 26, 56, 99, 99, 99, 99, 99},
{47, 66, 99, 99, 99, 99, 99, 99},
{99, 99, 99, 99, 99, 99, 99, 99},
{99, 99, 99, 99, 99, 99, 99, 99},
{99, 99, 99, 99, 99, 99, 99, 99},
{99, 99, 99, 99, 99, 99, 99, 99},
};
template
class Quantize {
public:
static Matrix8d forward(Matrix8d D) {
D = D.cwiseQuotient(QMatrix);
D = D.array().round(); // round before writing to bitstream
if(!((-128 <= D.array()).all() && (D.array() <= 127).all())){
std::cerr << "offender" << std::endl;
std::cerr << D;
assert((-128 <= D.array()).all() && (D.array() <= 127).all());
}
assert((-128 <= D.array()).all() && (D.array() <= 127).all());
return D;
}
static Matrix8d inverse(Matrix8d T) {
assert((-128 <= T.array()).all() && (T.array() <= 127).all());
T = T.cwiseProduct(QMatrix);
return T;
}
};
Appendix
This is the entire processing done to the video in both the decompressor and compressor.
// compressor.cpp
while (reader.read_next_frame()){
output_stream.push_byte(1); //Use a one byte flag to indicate whether there is a frame here
YUVFrame420& frame = reader.frame();
// 8x8 blocks for each frame
Matrix8d Yb(8, 8), Cbb(8,8), Crb(8,8);
// process Y frame
for (u32 y0 = 0; y0 < height; y0 += 8) {
for(u32 x0 = 0; x0 < width; x0 += 8) {
// fill up 8x8 block
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// original values range from 0 - 255
Yb(x % 8, y % 8) = static_cast(frame.Y(x, y));
Yb = Codec::DCT::forward(Yb);
Yb = Codec::QuantizeYb::forward(Yb);
// write 8x8 into the bitstream
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// emit transformed values into bitstream
output_stream.push_byte(static_cast(Yb(x % 8, y % 8)));
}
}
// process Cb frame
for (u32 y0 = 0; y0 < height/2; y0 += 8) {
for(u32 x0 = 0; x0 < width/2; x0 += 8) {
// fill up 8x8 block
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// original values range from 0 - 255
Cbb(x % 8, y % 8) = static_cast(frame.Cb(x, y));
Cbb = Codec::DCT::forward(Cbb);
Cbb = Codec::QuantizeYb::forward(Cbb);
// write 8x8 into the bitstream
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
output_stream.push_byte(static_cast(Cbb(x % 8, y % 8)));
}
}
// process Cr frame
for (u32 y0 = 0; y0 < height/2; y0 += 8) {
for(u32 x0 = 0; x0 < width/2; x0 += 8) {
// fill up 8x8 block
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// original values range from 0 - 255
Crb(x % 8, y % 8) = static_cast(frame.Cr(x, y));
Crb = Codec::DCT::forward(Crb);
Crb = Codec::QuantizeYb::forward(Crb);
// write 8x8 into the bitstream
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
output_stream.push_byte(static_cast(Crb(x % 8, y % 8)));
}
}
output_stream.push_byte(frame.Y(x,y));
}
//decompressor.cpp
while (input_stream.read_byte()){
YUVFrame420& frame = writer.frame();
// 8x8 blocks for each frame
Matrix8d Yb(8, 8), Cbb(8, 8), Crb(8, 8);
for (u32 y0 = 0; y0 < height; y0 += 8) {
for(u32 x0 = 0; x0 < width; x0 += 8) {
// fill up 8x8 block
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// transformed values range fom -128 - 127
Yb(x % 8, y % 8) = static_cast(input_stream.read_byte());
// write 8x8 into the bitstream
Yb = Codec::QuantizeYb::inverse(Yb);
Yb = Codec::DCT::inverse(Yb);
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
frame.Y(x, y) = static_cast(Yb(x % 8, y % 8));
}
}
for (u32 y0 = 0; y0 < height/2; y0 += 8) {
for(u32 x0 = 0; x0 < width/2; x0 += 8) {
// fill up 8x8 block
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// transformed values range fom -128 - 127
Cbb(x % 8, y % 8) = static_cast(input_stream.read_byte());
// write 8x8 into the bitstream
Cbb = Codec::QuantizeYb::inverse(Cbb);
Cbb = Codec::DCT::inverse(Cbb);
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
frame.Cb(x, y) = static_cast(Cbb(x % 8, y % 8));
}
}
for (u32 y0 = 0; y0 < height/2; y0 += 8) {
for(u32 x0 = 0; x0 < width/2; x0 += 8) {
// fill up 8x8 block
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
// transformed values range fom -128 - 127
Crb(x % 8, y % 8) = static_cast(input_stream.read_byte());
// write 8x8 into the bitstream
Crb = Codec::QuantizeYb::inverse(Crb);
Crb = Codec::DCT::inverse(Crb);
for(u32 y = y0; y < y0 + 8; y++)
for(u32 x = x0; x < x0 + 8; x++)
frame.Cr(x, y) = static_cast(Crb(x % 8, y % 8));
}
}
writer.write_frame();
}