[XviD-devel] New lumimasking code from ReferenceDivX
Dirk Knop
xvid-devel@xvid.org
Fri, 22 Nov 2002 12:49:50 +0100
This is a multi-part message in MIME format.
--------------060201070208060604030505
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
Hi,
ReferenceDivX (only reachable via IRC: irc.openprojects.net / channel
#xvid) sent me an adapt_quant.c which is very promising. It's algorithms
are based on HVS (e.g. edge-detection and stuff) and produce a _very_
nice output on my tests here (i.e. less blocking in dark areas).
The attached file needs the custom inter_matrix-code (last function in
adapt_quant.c) to be commented out as it doesn't work with dev-api-3
(but it works with stable-tree).
ReferenceDivX has some other ideas which can help lumi-masking a bit
more, but I forced him to set this code free ;)
Please test it and send reports!
Best regards,
Koepi
(in the name of refdivx)
--------------060201070208060604030505
Content-Type: text/plain;
name="adapt_quant.c"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="adapt_quant.c"
#include "../portab.h"
#include "adapt_quant.h"
#include "../dct/fdct.h"
#include "../dct/idct.h"
#include "../quant/quant_mpeg4.h"
#include "../quant/quant_h263.h"
#include "../encoder.h"
#include <math.h>
#include <stdlib.h> /* free, malloc */
#include <stdio.h>
#define M_PI 3.1415926536
#define RDIFF(a,b) ((int)(a+0.5)-(int)(b+0.5))
#define NJPEG 8 //'NJPEG' defines the JPEG's DCT block size (8x8)
#define sqr(X) ((X) * (X))
double **alloc_double(int width, int height) {
int i;
double **a = (double **) malloc(height * sizeof(double *));
for (i = 0; i < height; i++)
a[i] = (double *) malloc(width * sizeof(double));
return a;
}
/*
double **create_visibility_threshold(double viewing_distance, double image_width, double image_height, int pixel_width, int pixel_height, int blocksize) {
int i, j;
// image_width .. horizontal image size in centimeters
// pixel_width .. no. of pixels in image_width
// viewing_distance .. viewing distance in centimeters
double L,LT,Lk,Lf,Wx,Wy;
double So,at,ko,ak,fo,af;
double fio,foj,f,r,R;
double Tmin;
double teta;
double parabola_stepness;
double fmin;
double **thres = alloc_double(blocksize, blocksize);
// compute pixel size in radians
Wx = (image_width / pixel_width) / viewing_distance;
Wy = (image_height / pixel_height) / viewing_distance;
// convert pixel size in radians to size in degrees
Wx = Wx * 180.0 / M_PI;
Wy = Wy * 180.0 / M_PI;
//Wx = Wx * 180.0 / 3.14;
//Wy = Wy * 180.0 / 3.14;
// compute mean luminance in candels per square meter according to an exponential law
L = 240.0; // mean luminance in gray-levels, was 128
L = pow(10, (7.0 / 255.0 * L - 4.0));
// Compute luminance functions i.e. Tmin, fmin and K.
LT = 13.45, So = 94.7, at = 0.649;
Lk = 300, ko = 3.125, ak = 0.0706;
Lf = 300, fo = 6.78, af = 0.182;
// Minimal threshold function i.e. Tmin.
if (L > LT)
Tmin = L/So;
else
Tmin = L/So*pow(LT / L, 1 - at);
// Parabola stepness function i.e. K.
if (L > Lk)
parabola_stepness = ko;
else
parabola_stepness = ko * pow(L / Lk, ak);
// minimal frequency function i.e. fmin
if (L > Lf)
fmin = fo;
else
fmin = fo * pow(L / Lf, af);
for (i = 0; i < blocksize; i++)
for (j = 0; j < blocksize; j++) {
fio = i;
foj = j;
fio /= (double) (2.0 * blocksize);
foj /= (double) (2.0 * blocksize);
fio /= Wx;
foj /= Wy;
// compute spatial frequencies
f = sqrt(sqr(fio) + sqr(foj));
// compute angular parameters
teta = 2.0 * fio * foj;
if (i > 0 && j > 0) {
teta /= sqr(f);
}
if (teta > 1.0) teta = 1.0;
teta = asin(teta);
// determine obliqueness effect values
R = 0.7;
r = R + (1 - R) * sqr(cos(teta));
// Compute luminance thresholds i.e. thres.
// f[0][0] is 0 but to don't have log(0) we make f[0][0] = 1.
if (i == 0 && j == 0)
f = 1;
thres[i][j] =
log10(Tmin / r) +
parabola_stepness * pow(log10(f / fmin), 2);
thres[i][j] = pow(10, thres[i][j]);
// According to our experiments we can multiply the thresholds
// until a factor of 3.5 so that they are not yet visibly after luminance
// masking corrections.
thres[i][j] = 2.2 * thres[i][j];
}
// Approximate luminance threshold for DC coefficient
thres[0][0] = MIN(thres[0][1], thres[1][0]);
return thres;
}
*/
/*
void contrast_masking(double block[][8], double **thres, int blocksize, double mask[][8] )
{
int i, j;
double actthres,lumthres;
// double **mask = (double **) alloc_double(blocksize, blocksize);
// Contrast masking degree i.e. w.
double w = 0.7;
// Compute the mean luminance in DCT domain.
double mean_lum = 128 * blocksize;
// Compute luminance masking threshods
double lum_thres = fabs(block[0][0]);
if (lum_thres == 0 || lum_thres == blocksize)
return NULL;
for (i = 0; i < blocksize; i++) {
for (j = 0; j < blocksize; j++) {
// the factor 2.2 was added to increase the luminance masking
lumthres = 3.2 * thres[i][j] * pow(lum_thres/mean_lum, 0.65);
if (lum_thres < 600)
lumthres = lumthres * 2.2; // very dark areas can be marked a bit more
actthres = pow(abs(block[i][j]), w) * pow(lumthres, (1 - w));
if (i == 0 && j == 0)
actthres = lumthres;
// Compute contrast masking thresholds
mask[i][j] = lumthres;
if (lumthres < actthres)
mask[i][j] = actthres;
}
}
}
*/
int **alloc_int(int width, int height) {
int i;
int **a = (int **) malloc(height * sizeof(int *));
for (i = 0; i < height; i++)
a[i] = (int *) malloc(width * sizeof(int));
return a;
}
void free_int(int **a, int height) {
int i;
for (i = 0; i < height; i++)
free(a[i]);
free(a);
}
int
normalize_quantizer_field(float *in,
int *out,
int num,
int min_quant,
int max_quant)
{
int i;
int finished;
do {
finished = 1;
for (i = 1; i < num; i++) {
if (RDIFF(in[i], in[i - 1]) > 2) {
in[i] -= (float) 0.5;
finished = 0;
} else if (RDIFF(in[i], in[i - 1]) < -2) {
in[i - 1] -= (float) 0.5;
finished = 0;
}
if (in[i] > max_quant) {
in[i] = (float) max_quant;
finished = 0;
}
if (in[i] < min_quant) {
in[i] = (float) min_quant;
finished = 0;
}
if (in[i - 1] > max_quant) {
in[i - 1] = (float) max_quant;
finished = 0;
}
if (in[i - 1] < min_quant) {
in[i - 1] = (float) min_quant;
finished = 0;
}
}
} while (!finished);
out[0] = 0;
for (i = 1; i < num; i++)
out[i] = RDIFF(in[i], in[i - 1]);
return (int) (in[0] + 0.5);
}
void microfilter(int* val1, int* val2)
{
int temp1=0;
temp1 = abs(*val1 - *val2);
temp1 /= 2;
temp1 = ((temp1 * .5)+.5);
if (temp1 > 3)
temp1 = 3;
if (val1 >= val2)
{
*val1 -= temp1;
*val2 += temp1;
}
else
{
*val1 += temp1;
*val2 -= temp1;
}
}
int
adaptive_quantization(unsigned char *buf,
int stride,
int *intquant,
int framequant,
int min_quant,
int max_quant,
int mb_width,
int mb_height) // no qstride because normalization
{
int i, j, k, l;
int n1=0, n2=0, n3=0, n4=0; //new code for lumblock counters
int x=0, y=0;
int thresh=100; // bigger number means less edges for sobel
int sx, sy; //stuff for sobel
int xx, yy; //stuff for sobel
const int kx[3][3] =
{ { -1, -2, -1},
{ 0, 0, 0},
{ 1, 2, 1}};
const int ky[3][3] =
{ { -1, 0, 1},
{ -2, 0, 2},
{ -1, 0, 1}};
unsigned char *ptr;
float *quant;
short *lumblock1;
short *lumblock2;
short *lumblock3;
short *lumblock4;
short *lumblock1a;
short *lumblock2a;
short *lumblock3a;
short *lumblock4a;
short *qcoeff1;
short *qcoeff2;
short *qcoeff3;
short *qcoeff4;
//int **imageIn = alloc_int(16, 16);
//int **imageOut = alloc_int(16, 16);
int imageIn[16][16];
int imageOut[16][16];
int val1a=0, val2a=0, lastquant=0;
int ji=0, jj=0, jl=0, jm=0, jk=0, tempquant=2, test=0;
float fmicroblock=0, fedgeerror=0.0;
int edge=0, iedgeerror=0, imicroblock=0;
if (!(quant = (float *) malloc(mb_width * mb_height * sizeof(float))))
return(-1);
if(!(lumblock1 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock2 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock3 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock4 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock1a = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock2a = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock3a = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(lumblock4a = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(qcoeff1 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(qcoeff2 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(qcoeff3 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
if(!(qcoeff4 = (short *) malloc(8 * 8 * sizeof(short))))
return(-1);
//start of macroblock loop
for (k = 0; k < mb_height; k++) {
for (l = 0; l < mb_width; l++) // do this for all macroblocks individually
{
quant[k * mb_width + l] = (float)framequant; //sets each quant to the constant
ptr = &buf[16 * k * stride + 16 * l]; // address of MB
n1=n2=n3=n4=0; //resets counters
for (i = 0; i < 16; i++) {
for (j = 0; j < 16; j++) {
imageIn[i][j] = ptr[i * stride + j]; //orginal image
if (j < 8)
{
if (i < 8)
{
lumblock1[n1] = ptr[i * stride + j];
lumblock1a[n1] = ptr[i * stride + j];
n1++;
}
else
{
lumblock3[n3] = ptr[i * stride + j];
lumblock3a[n3] = ptr[i * stride + j];
n3++;
}
}
else
{
if (i < 8)
{
lumblock2[n2] = ptr[i * stride + j];
lumblock2a[n2] = ptr[i * stride + j];
n2++;
}
else
{
lumblock4[n4] = ptr[i * stride + j];
lumblock4a[n4] = ptr[i * stride + j];
n4++;
}
}
}
}
// sobel edge finder to array using 3x3 matrix
//////////////////////////////////////////////////////////
thresh=100; //edge threshold
for(x=1;x<(16-1);x++)
for(y=1;y<(16-1);y++) {
sx = 0; //resets back to 0
for(xx=-1;xx<=1;xx++)
for(yy=-1;yy<=1;yy++)
sx+= imageIn[x+xx][y+yy] * kx[xx+1][yy+1];
sy = 0; //resets back to 0
for(xx=-1;xx<=1;xx++)
for(yy=-1;yy<=1;yy++)
sy+= imageIn[x+xx][y+yy] * ky[xx+1][yy+1];
if(sqrt(sx*sx+sy*sy)>thresh)
imageOut[x][y] = 1;
else
imageOut[x][y] = 255;
}
//////////////////////////////////////////////////////////
//start new code
//we should start the quant at the lowest possible value
tempquant = min_quant;
test = 0;
do {
n1=n2=n3=n4=0; //resets counters
fmicroblock=0;
imicroblock=0;
fdct(&lumblock1[0]);
quant4_intra(&qcoeff1[0], &lumblock1[0], tempquant, 1);
dequant4_intra(&lumblock1a[0], &qcoeff1[0], tempquant, 1);
idct(&lumblock1a[0]);
fdct(&lumblock2[0]);
quant4_intra(&qcoeff2[0], &lumblock2[0], tempquant, 1);
dequant4_intra(&lumblock2a[0], &qcoeff2[0], tempquant, 1);
idct(&lumblock2a[0]);
fdct(&lumblock3[0]);
quant4_intra(&qcoeff3[0], &lumblock3[0], tempquant, 1);
dequant4_intra(&lumblock3a[0], &qcoeff3[0], tempquant, 1);
idct(&lumblock3a[0]);
fdct(&lumblock4[0]);
quant4_intra(&qcoeff4[0], &lumblock4[0], tempquant, 1);
dequant4_intra(&lumblock4a[0], &qcoeff4[0], tempquant, 1);
idct(&lumblock4a[0]);
idct(&lumblock1[0]);
idct(&lumblock2[0]);
idct(&lumblock3[0]);
idct(&lumblock4[0]);
for(ji = 0; ji< 8; ji++)
{
imicroblock += abs(lumblock1a[(((ji +1)*8)-1)] - lumblock2a[(8*ji)]);
imicroblock += abs(lumblock3a[(((ji +1)*8)-1)] - lumblock4a[(8*ji)]);
imicroblock += abs(lumblock2a[(63-ji)] - lumblock4a[(7 - ji)]);
imicroblock += abs(lumblock1a[(63-ji)] - lumblock3a[(7 - ji)]);
}
fmicroblock = (float)(imicroblock / 32);
edge=1;
fedgeerror=0;
iedgeerror=0;
for(x=0;x<16;x++)
for(y=0;y<16;y++) {
if (imageOut[x][y] == 1)
{
edge++;
//then find the quad
if (y < 8)
{
if (x < 8)
{
//Quant 1 j(0-7) j(0-7)
iedgeerror += abs(lumblock1a[n1] - imageIn[x][y]);
n1++;
}
else
{
//Quant 3 j(8-15) i(0-7)
iedgeerror += abs(lumblock3a[n3] - imageIn[x][y]);
n3++;
}
}
else
{
if (x < 8)
{
//Quant 2 j(0-7) i(7-15)
iedgeerror += abs(lumblock2a[n2] - imageIn[x][y]);
n2++;
}
else
{
//Quant 4 j(7-15) i(7-15)
iedgeerror += abs(lumblock4a[n4] - imageIn[x][y]);
n4++;
}
}
}
else
{
if (y < 8)
{
if (x < 8)
{
//Quant 1 j(0-7) j(0-7)
n1++;
}
else
{
//Quant 3 j(8-15) i(0-7)
n3++;
}
}
else
{
if (x < 8)
{
//Quant 2 j(0-7) i(7-15)
n2++;
}
else
{
//Quant 4 j(7-15) i(7-15)
n4++;
}
}
}
}
fedgeerror = (float)(iedgeerror / edge);
if (((edge > 128) && (fedgeerror > 5)) || (fmicroblock > 5) || (tempquant == max_quant) || (((l!=0) && ((lastquant+2) > tempquant ) || (lastquant-2) < tempquant )) )
//
{
if ((tempquant -1) < 2)
lastquant=2;
else
lastquant= (tempquant -1);
quant[k * mb_width + l] = (float)(lastquant);
test = 1;
}
else
{
test = 0;
tempquant++;
}
} while (test == 0);
n1=0;
n2=0;
n3=0;
n4=0; //resets counters
//
for (i = 0; i < 16; i++) {
for (j = 0; j < 16; j++) {
if (j == 7)
{
if (i < 8)
{
val1a = lumblock1[(((n1 +1)*8)-1)];
val2a = lumblock2[(8*n1)];
n1++;
}
else
{
val1a = lumblock3[(((n2 +1)*8)-1)];
val2a = lumblock4[(8*n2)];
n2++;
}
microfilter(&val1a, &val2a);
ptr[i * stride + j] = val1a;
ptr[i * stride + (j+1)] = val2a;
}
if (i == 7)
{
if (j < 8)
{
val1a = lumblock1[(56+n4)];
val2a = lumblock3[(n4)];
n4++;
}
else
{
val1a = lumblock2[(56+n3)];
val2a = lumblock4[(n3)];
n3++;
}
microfilter(&val1a, &val2a);
ptr[i * stride + j] = val1a;
ptr[(i+1) * stride + j] = val2a;
}
}
}
//
}
}
i = normalize_quantizer_field(quant, intquant,
mb_width * mb_height,
min_quant, max_quant);
free(quant);
free(lumblock1);
free(lumblock2);
free(lumblock3);
free(lumblock4);
free(lumblock1a);
free(lumblock2a);
free(lumblock3a);
free(lumblock4a);
free(qcoeff1);
free(qcoeff2);
free(qcoeff3);
free(qcoeff4);
return(i);
}
//adapts inter matrix code
int adaptive_matrix2a(unsigned char *buf,
int stride,
int min_quant,
int mb_width,
int mb_height) // no qstride because normalization
{
int i, j, k, l;
float cvalue=18;
float *quant;
unsigned char *ptr;
float *val;
uint8_t *tempquantout;
float global = 0.;
uint32_t mid_range = 0;
const float DarkAmpl = 14 / 2;
const float BrightAmpl = 10 / 2;
const float DarkThres = 70;
const float BrightThres = 200;
const float GlobalDarkThres = 70; //adjusts one place up
const float GlobalBrightThres = 150; //adjusts one place up
const float GlobalDarkThres2 = 35; //adjusts two place up
const float GlobalBrightThres2 = 190; //adjusts two place up
const float MidRangeThres = 20;
const float UpperLimit = 200;
const float LowerLimit = 25;
double hvsfwm[64]={ 1,1,1,1,1,1,1,1,
1,1,.9599,1,.9571,1,.9599,
.8746,.9283,.8898,.8898,.9283,
.8746,.7684,.8404,.8192,.7617,
.8192,.8404,.7684,.6571,.7371,
.7371,.6669,.6669,.7371,.7371,
.6571,.6306,.6471,.5912,.5419,
.5912,.6471,.6306,.5558,.5196,
.4564,.4564,.5196,.5558,.4495,
.3930,.3598,.3930,.4495,.3393,
.2948,.2948,.3393,.2480,.2278,
.2480,.1828,.1828,.1391};
int order[64]={ 1,2,6,7,15,16,28,29,
3,5,8,14,17,27,30,43,
4,9,13,18,26,31,42,44,
10,12,19,25,32,41,45,54,
11,20,24,33,40,46,53,55,
21,23,34,39,47,52,56,61,
22,35,38,48,51,57,60,62,
36,37,49,50,58,59,63,64};
if (min_quant == 2)
cvalue=20.5;
else if (min_quant == 3)
cvalue=19.5;
else if (min_quant == 4)
cvalue=19;
else if (min_quant == 5)
cvalue=18;
else if (min_quant >= 6)
cvalue=17;
if(!(val = (float *) malloc(mb_width * mb_height * sizeof(float))))
return(-1);
if(!(tempquantout = (uint8_t *) malloc(8 * 8 * sizeof(uint8_t))))
return(-1);
for (k = 0; k < mb_height; k++) {
for (l = 0; l < mb_width; l++) // do this for all macroblocks individually
{
// quant[k * mb_width + l] = (float) framequant;
// calculate luminance-masking
ptr = &buf[16 * k * stride + 16 * l]; // address of MB
val[k * mb_width + l] = 0.;
for (i = 0; i < 16; i++)
for (j = 0; j < 16; j++)
val[k * mb_width + l] += ptr[i * stride + j];
val[k * mb_width + l] /= 256.; //divides total by 16*16 to get average
global +=val[k * mb_width + l];
}
}
global /=mb_width * mb_height;
if ((global > GlobalBrightThres) || (global < GlobalDarkThres))
cvalue++; // it is pretty brite then have factor go up a bit
if ((global > GlobalBrightThres2) || (global < GlobalDarkThres2))
cvalue=cvalue+2; // it is pretty brite then have factor of 2
if (cvalue > 20.5)
cvalue=20.5;
k=0;
for(i=0; i<8; i++)
{
for(j=0; j<8; j++)
{
tempquantout[k] = (int)((cvalue / hvsfwm[(order[k])-1])+.5);
//tempquantout[k] = (int)((cvalue / hvsfwm[(order[k])-1])+.5);
k++;
}
}
i = set_inter_matrix1a(tempquantout);
free(val);
free(tempquantout);
return(i);
}
--------------060201070208060604030505--