📄 dct.c
字号:
/*********
we must transmit the quant-tables if we customize them.
This will crap out performance on small files.
----
tail smashing.
preliminary experiments on checa show that tail smashing
only helps for ratio ~ 100 , and is irrelevant around ratio ~ 1
quality Q @ no tail smash
equiv
quality (Q+1) @ tail smash J *= 60
----
notes:
the "rateD" computation is *NEGATIVE* (indicating that a lower
rate gives *better* quality) and some sub-sections of the
image! This indicates that some locally-adaptive quantization
could be smashing.
*********/
#define SMASH_SOUGHT_MIN 2
#define SMASH_DISTORTION_DIV 15
#define SMASH_SOUGHT_MULT (0.4)
//#define SMASH_SOUGHT_MULT (99)
//#define SMASH_SOUGHT_MULT (tune_param/10.0)
/** surprisingly, small changes in these params
don't affect the R-D curves at all! **/
#define VERBOSE
#define DO_LOG
#define DEBUG
//#define FIDDLE_QUANTIZERS
//#define NO_CODING
#define VISUAL_QUANTIZERS // else uniform, ** hurts the raw numbers!
#define DPCM_ORDER_M1 // else order0, ** order(-1) is safer
//#define DPCM_ORDER_0 // else order1 : order1 is best on Lena
#define FLAT_ORDER_M1 // else jpeg : flat is best for o(-1)
//#define JPEG_ORDER_M1 // else custom
#define USE_EOBS // ** helps about 0.02 bpb
#define CODE_SIGNS // helps on files bigger than 256x256
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <crbinc/inc.h>
#include <crbinc/fileutil.h>
#include <crbinc/cindcatr.h>
#include <crbinc/timer.h>
#include <crbinc/bbitio.h>
#include <crbinc/arithc.h>
#include <crbinc/scontext.h>
#include <crbinc/context.h>
#include "dct.h"
static const int qtable_jpeg[DCTBLOCK] =
{ 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 };
static const int qtable_jpeg_flattened[DCTBLOCK] =
{ 18, 23, 22, 28, 34, 44, 50, 55,
24, 24, 26, 31, 36, 53, 54, 52,
26, 25, 28, 34, 44, 53, 58, 52,
26, 29, 33, 38, 50, 65, 63, 55,
30, 33, 43, 52, 58, 73, 71, 61,
34, 41, 52, 56, 63, 71, 74, 67,
49, 56, 62, 65, 71, 77, 77, 70,
59, 67, 68, 69, 74, 70, 71, 70};
static const int qtable_semiuniform[DCTBLOCK] =
{ 16, 32, 64, 64, 64, 64, 64, 64,
32, 64, 64, 64, 64, 64, 64, 64,
64, 64, 64, 64, 64, 64, 64, 64,
64, 64, 64, 64, 64, 64, 64, 64,
64, 64, 64, 64, 64, 64, 64, 64,
64, 64, 64, 64, 64, 64, 64, 64,
64, 64, 64, 64, 64, 64, 64, 64,
64, 64, 64, 64, 64, 64, 64, 64 };
static const int qtable_uniform[DCTBLOCK] =
{ 100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 100, 100, 100, 100, 100, 100 };
static const int qtable_none[DCTBLOCK] =
{ 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1 };
static const int qtable_jpeg_squashed[DCTBLOCK] = { /** 4*sqrt(jpeg_flat) **/
18, 19, 19, 21, 23, 27, 28, 30,
20, 20, 21, 22, 24, 29, 30, 29,
21, 20, 21, 23, 27, 29, 31, 29,
21, 22, 23, 25, 28, 32, 32, 30,
22, 23, 26, 29, 31, 34, 34, 31,
23, 26, 29, 30, 32, 34, 35, 33,
28, 30, 32, 32, 34, 35, 35, 34,
31, 33, 33, 33, 35, 34, 34, 34 };
#ifdef VISUAL_QUANTIZERS
//#define quant_table qtable_jpeg
#define quant_table qtable_jpeg_squashed
#else // UNIFORM QUANTIZERS
#define quant_table qtable_semiuniform
//#define quant_table qtable_uniform /** smashing is a clear win here
#endif
static const int zigzag[DCTBLOCK] = {
0, 1, 8, 16, 9, 2, 3, 10,
17, 24, 32, 25, 18, 11, 4, 5,
12, 19, 26, 33, 40, 48, 41, 34,
27, 20, 13, 6, 7, 14, 21, 28,
35, 42, 49, 56, 57, 50, 43, 36,
29, 22, 15, 23, 30, 37, 44, 51,
58, 59, 52, 45, 38, 31, 39, 46,
53, 60, 61, 54, 47, 55, 62, 63 };
/*** unzigzag : (never used)
{ 0, 1, 5, 6, 14, 15, 27, 28, 2, 4, 7, 13, 16, 26, 29, 42, 3, 8, 12, 17,
25, 30, 41, 43, 9, 11, 18, 24, 31, 40, 44, 53, 10, 19, 23, 32, 39, 45, 52, 54,
20, 22, 33, 38, 46, 51, 55, 60, 21, 34, 37, 47, 50, 56, 59, 61, 35, 36, 48, 49,
57, 58, 62, 63 };
*******/
static int zag_x[DCTBLOCK] = { // (zigzag & 7)
0, 1, 0,
0, 1, 2, 3, 2, 1, 0,
0, 1, 2, 3, 4, 5, 4, 3, 2, 1, 0,
0, 1, 2, 3, 4, 5, 6, 7, 6, 5, 4, 3, 2, 1, 0,
1, 2, 3, 4, 5, 6, 7, 7, 6, 5, 4, 3, 2,
3, 4, 5, 6, 7, 7, 6, 5, 4,
5, 6, 7, 7, 6, 7 };
static int zag_y[DCTBLOCK] = { // (zigzag>>3)
0,
0, 1, 2, 1, 0,
0, 1, 2, 3, 4, 3, 2, 1, 0,
0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1, 0,
0, 1, 2, 3, 4, 5, 6, 7, 7, 6, 5, 4, 3, 2, 1,
2, 3, 4, 5, 6, 7, 7, 6, 5, 4, 3,
4, 5, 6, 7, 7, 6, 5,
6, 7, 7 };
#define cleanup(str) if ( 0 ) ; else { errputs(str); exit(10); }
#ifdef DO_LOG
#define LOG(x) if (0) ; else x;
#else
#define LOG(x) if (1) ; else x;
#endif
#define TRANS_MAX (0xFFF)
#define TRANS_MAX_EV (0xFFE)
#define TRANS_TOPBIT (0x800)
#define TRANS_LOWBIT (0x002) // since 001 is the sign
#define VAL_MAX (TRANS_TOPBIT)
#define VAL_EOB (VAL_MAX)
#define TRANS_EOB (TRANS_MAX+1)
#define TRANS_DONT_CODE (TRANS_MAX+2)
#define data_trans(val) ( ((val) < 0 ) ? ( -1+min(TRANS_MAX_EV,(-val-val))) : (min(TRANS_MAX_EV,(val+val))) )
#define trans_data(val) (( (val) == TRANS_EOB || (val) == TRANS_DONT_CODE ) ? 0 : ( ((val) & 1 ) ? (-((val+1)>>1)) : (((val)>>1)) ))
int log_n0,log_n1,log_nm1;
int log_zero_len,log_zero_len_smashed,log_zero_len_gained;
int log_0p0,log_0p1,log_0pbig;
int log_nsign_noc,log_nsign_same,log_nsign_diff;
int log_nsign,log_nsign_dpcm;
int log_max_RoDs5,log_min_RoDs5;
#define RESET_LOG() if ( 0 ) ; else { log_n0=log_n1=log_nm1=log_zero_len=log_zero_len_smashed=log_zero_len_gained=log_0p0=log_0p1=log_0pbig=log_nsign_noc=log_nsign_same=log_nsign_diff=log_nsign=log_nsign_dpcm=log_max_RoDs5=log_min_RoDs5=0; };
#define CODE_CONTEXTS 9
static const int context_buckets[CODE_CONTEXTS] =
{ 0,1,2,3,4,8,16,32, TRANS_MAX + 10 }; // doesn't make too much difference
#define ORDER1_TOTMAX 6000
#define ORDER1_ALPHABET 30 // infinity is optimum, but slow for scontext
#ifdef USE_EOBS
#define ORDER1_ESCAPE (ORDER1_ALPHABET-2)
#define ORDER1_EOB (ORDER1_ALPHABET-1)
#else
#define ORDER1_ESCAPE (ORDER1_ALPHABET-1)
#endif // USE_EOBS
#define ORDER1_INC 15 // and 30 is only 0.001 bpb worse than infinity
#define ORDER0_TOTMAX 2500
#define ORDER0_ALPHABET (VAL_MAX + 2 - ORDER1_ALPHABET)
#ifdef CODE_SIGNS
#define SIGN_CONTEXTS 3
#define SIGNORDER0_TOTMAX 1000
#define SIGNORDER0_INC 5
#endif
#define FLAT_ORDER_M1_ESCAPE (200)
#define M1_THRESH_1 8
#define M1_THRESH_2 24
#define M1_THRESH_3 40
#define M1_THRESH_4 100
#define PSNR_MAX (48.165) //(10*log10(256^2))
int tune_param = 1;
ubyte **raw = NULL,**rawout = NULL;
uword **trans = NULL;
ubyte *comp = NULL,*signcomp = NULL,*rawcomp = NULL;
FILE *rawF = NULL;
struct FAI * FAI = NULL;
struct FAI * signFAI = NULL;
struct FAI * rawFAI = NULL;
struct BitIOInfo * BII = NULL;
struct BitIOInfo * signBII = NULL;
struct BitIOInfo * rawBII = NULL;
context *order0 = NULL;
#ifdef CODE_SIGNS
scontext **signorder0 = NULL;
#endif
scontext **order1 = NULL;
int num_planes;
void ExitFunc(void);
void encode(int sym,int context);
void encode_m1_jpeg(int sym);
void encode_m1_flat(int sym);
void encode_m1_custom(int sym);
void encodeSign(bool sym,int prev);
int code_image(ubyte **raw,ubyte **rawout,int width,int height,int num_planes,
int * qtable,int J_sought_s5,int max_distortion,
int * comp_dpcm_len_ptr,int *comp_ac_len_ptr,int *comp_sign_len_ptr,int *comp_raw_len_ptr);
void dct_image(ubyte **raw,uword **trans,int width,int height,int num_planes,
int * quant_table);
void idct_image(ubyte **rawout,uword **trans,int width,int height,int num_planes,
int * quant_table);
void TheImageAnalyzer(ubyte **original,ubyte **vs,
int num_planes,int width,int height,
float ratio,FILE *sio);
void init_allocs(int num_planes,int rawsize,int complen);
void init_coders(void);
int done_coders(void);
void free_coders(void);
int main(int argc,char *argv[])
{
char *t1,*t2;
int width,height;
int rawsize,complen,totsize;
int comp_tot_len,comp_dpcm_len,comp_ac_len,comp_sign_len,comp_raw_len,i;
clock_t startclock,stopclock;
int quality,quantizer,J_sought_s5,max_distortion;
float J_sought;
int qtable[DCTBLOCK];
errputs("dct v0.5, (c) 1998 by Charles Bloom");
if ( argc <= 3 ) {
errputs("dct");
errputs("dct <raw file> <width> <height> [planes] [quality] [tune param]");
exit(10);
}
if ( atexit(ExitFunc) )
cleanup("Couldn't install exit trap!");
if ( ( rawF = fopen(argv[1],"rb") ) == NULL )
cleanup("open raw failed");
width = atol(argv[2]);
height = atol(argv[3]);
if ( ((width/DCTLINE)*DCTLINE) != width )
cleanup("size must be divable by 8");
if ( ((height/DCTLINE)*DCTLINE) != height )
cleanup("size must be divable by 8");
if ( argc > 4 )
num_planes = atol(argv[4]);
else
num_planes = 1;
if ( argc > 5 )
quality = min(max(1,atol(argv[5])),100);
else
quality = 50;
if ( argc > 6 )
tune_param = atol(argv[6]);
else
tune_param = 1;
rawsize = width * height;
totsize = rawsize*num_planes;
complen = (totsize*9)>>3;
init_allocs(num_planes,rawsize,complen);
for(i=0;i<num_planes;i++) {
if ( !FReadOk(rawF,raw[i],rawsize) )
errputs("fread short! continuing..");
}
fclose(rawF); rawF = NULL;
if ( quality >= 50 ) quantizer = (100 - quality)*2;
else quantizer = 5000/quality;
for(i=0;i<DCTBLOCK;i++) {
qtable[i] = (quant_table[i]*quantizer + 40)/100;
if ( qtable[i] < 1 ) qtable[i] = 1;
zag_y[i] *= width;
}
/**** measure the R_D slope for a change of "quality" ****/
#ifdef USE_EOBS // otherwise this is irrelevant
do {
int qtn[DCTBLOCK];
int base_len,new_len,base_err,new_err,i,j,d;
double dRoD;
if ( quality > 99 )
{ J_sought_s5 = 9999; max_distortion=0; break; }
if ( quality >= 50 ) quantizer = (100 - (quality-1))*2;
else quantizer = 5000/(quality-1);
for(i=0;i<DCTBLOCK;i++) {
qtn[i] = (quant_table[i]*quantizer + 40)/100;
if ( qtn[i] < 1 ) qtn[i] = 1;
}
base_len = code_image(raw,rawout,width,height,num_planes,qtable,1000,0,
NULL,NULL,NULL,NULL);
base_err=0;
for(i=0;i<num_planes;i++)
for(j=0;j<rawsize;j++) {
d = raw[i][j] - rawout[i][j];
base_err += d*d;
}
new_len = code_image(raw,rawout,width,height,num_planes,qtn,1000,0,
NULL,NULL,NULL,NULL);
new_err=0;
for(i=0;i<num_planes;i++)
for(j=0;j<rawsize;j++) {
d = raw[i][j] - rawout[i][j];
new_err += d*d;
}
max_distortion = base_err - new_err;
if ( max_distortion == 0 ) max_distortion++;
if ( max_distortion < 0 ) max_distortion = - max_distortion;
dRoD = 8*(double)(new_len - base_len)/sqrt((double)max_distortion);
if ( dRoD < 0 ) dRoD = - dRoD ;
#ifdef VERBOSE
/*#*/{
float base_rmse,new_rmse;
base_rmse= sqrt((float)base_err/totsize);
new_rmse= sqrt((float)new_err/totsize);
printf("base: len= %d , err= %.2f (rmse=%.2f) , new: len= %d , err= %.2f (%.2f)\n",base_len,(float)base_err/totsize,base_rmse,new_len,(float)new_err/totsize,new_rmse);
printf("dRoD = %f , perfo(rmse): base =%.3f, new=%.3f \n",dRoD,((float)totsize/base_len)/base_rmse,((float)totsize/new_len)/new_rmse);
/*#*/}
#endif
J_sought_s5 = (dRoD*32.0*SMASH_SOUGHT_MULT);
if ( J_sought_s5 < SMASH_SOUGHT_MIN ) J_sought_s5 = SMASH_SOUGHT_MIN;
max_distortion /= SMASH_DISTORTION_DIV;
} while(0);
#endif // USE_EOBS
/*** fiddle quantizers ***/
#ifdef FIDDLE_QUANTIZERS
do {
int b,d,j;
int r_c,r_p,r_m,d_c,d_p,d_m,inc;
int tot_wiggle=0,num_wiggles=0;
bool active[DCTBLOCK];
int nactive=64;
/** make free changes which don't affect distortion,
but decrease the rate ;
very slow and helps just a tad. (Ratio/RMSE += 0.01) **/
for(b=0;b<DCTBLOCK;b++)
active[b] = true;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -