[实用向] JPEG图像压缩详解

ywy_c_asm · 2021-09-06 23:47:17 · 个人记录

前言

最简单的图像存储格式是BMP位图，它直接把每个位置的像素RGB值存起来，体积会很大。

那么如何把一个图像进行压缩呢？我们要知道，图像本质是连续的二维信息（当然，除非这图像全是乱码，不是给人看的……），像素值只是对它的离散采样而已。BMP位图只是把图像简简单单地当成一个二维数组，并没有很好的利用其连续性。

举个例子，假如现在有一幅2000×2000的被红色填满的图像，用BMP存储，图像大小至少是2000×2000×3Bytes=11MB，可是像素值全是(255,0,0)，这样显然不优美。当然，更多的例子则是，图像的相邻像素颜色值的差别一般都不大。

JPEG图像压缩的原理就是利用这种连续性，再结合人眼的一些特性，去除掉一些人眼不敏感的信息，把原图用尽量小的整数表示，最后用变长的二进制编码存储它们。这种压缩是有损的。

本文共分为两部分，第一部分是JPEG压缩算法流程的介绍，第二部分是JPEG文件存储格式的介绍。毕竟，我最近在手动实现JPEG编/解码器的时候找了大量资料，发现这两个重要的内容往往是割裂的，光知道这算法咋搞的也没用，它在实际上会以一种十分诡异的方式存储，错一点就挂掉了，导致我花了大量时间去探索。

JPEG图像压缩算法

1. 颜色的变换

大家都知道，计算机内一般以RGB值表示一个像素的颜色，我们把这个称作RGB颜色空间，在其中三种颜色分量被视作地位相等的（所以，位图必须以相同的形式把3种颜色分量都存起来）。

然而事实真的如此吗？并不。由于人体的视神经结构，人眼对于亮度更加敏感。而当颜色值相同时，绿色的亮度要比红色和蓝色都要明显的高。所以就有这样一种不同于RGB的基于亮度的YUV颜色空间（一般和YCbCr空间指同一个东西），Y分量表示亮度，U和V分量共同表示色度。RGB与YUV的转换关系如下：

（注意：由于U和V可能是负的，在实际操作中，需要把U和V分别加上128）

这两种颜色空间各自都有长处，RGB对机器十分友好，易于存储，而YUV更加人性化，考虑到了人的视觉特性。

JPEG在压缩时，先把颜色变换到YUV，由于人眼对Y的敏感程度远高于U和V，所以一般使用Y:U:V=2:1:1的采样方式。啥意思呢？每个格子都得提供Y，但U和V则是把图像划分为若干2×2的连续子矩阵，每个矩阵里仅取一个U和V（一般是平均值），这样就利用了图像的连续性以及人眼的敏感特性，使得U和V少存了一半。这也是其“有损”的第一个来源。

对于Y、U、V三个分量，我们需要分别把它们划分成若干个8×8的连续子矩阵。也就是说，我们需要把图像的宽度和高度补成16的倍数（毕竟U和V是两行/列一采样）。每个8×8的子矩阵是JPEG存储的基本单位，我们后续的工作就是如何把这些8×8矩阵压缩并存储，基本上就不考虑它是哪个分量提供的了。

2. DCT变换以及量化

众所周知，一维的信号可以通过傅里叶变换转成频域，这样就能把一个连续信号拆成若干正/余弦波的叠加组合。其实二维的图像也可以如此，只不过拆成的是形如\sin(ax+by)这样的二维平面波而已。低频波反映了图像的主体，高频波反映了图像的细节，所以人眼对低频波的敏感程度要高于高频波。

（感兴趣的同学可以去尝试用二维FFT对图像做变换，然后尝试把低频波或高频波去掉，再IDFT回去看看图像变成了啥，挺有意思的，这里就不展开讲了，毕竟JPEG也不需要FFT……）

而DCT（离散余弦变换）是离散傅里叶变换的一种特殊形式，它相当于把原信号做了偶延拓，仅用余弦表示原信号，并不会涉及复数（在DFT中，复数的实部是余弦的系数，虚部是正弦的系数）。

现在我们需要对这些8×8矩阵分别做DCT，做法是用矩阵乘法，T_{8×8}[i][j]=\left\{\begin{aligned}\frac 1{\sqrt 8}\space(i=0) \\ \frac{\cos(\frac{\pi(2j+1)i}{16})}{2}\space(i>0)\end{aligned}\right.，那么DCT(A)=TAT^T。

下面是一个实际的DCT后的矩阵的例子：

可以看到，左上角第一个元素的值是明显高于右下方元素的，基本上所有正常的图像都是如此，低频波远高于高频波。那个第一个元素其实就是两个维度振幅都为0的余弦波，其实也就是常数，我们把它称为直流系数，其余元素都表示有振幅的余弦波，被称为交流系数。可以发现，直流系数是对图像影响较大的因子，于是我们可以采取和之前的YUV类似的思想，把一些不太重要的高频信息过滤掉，于是就需要量化操作。

所谓量化操作其实就是一个8×8常量矩阵，让DCT后的矩阵里的所有数去与这个量化矩阵进行对应位相除，再四舍五入取整（这也是JPEG的第二个“有损”的来源），得到一个整数矩阵，我们实际压缩的就是这个整数矩阵。

一般的，对于Y分量和UV分量的矩阵，会采取两个不同的量化矩阵，以下是Windows画图采取的Y量化矩阵和UV量化矩阵（这个是我在画图保存的jpg文件的二进制信息里提取的）：

可以发现，量化表的右下部分的值会比左上值大，这其实就相当于对高频的滤波。除完后矩阵内的整数会缩小，并且有很多位置会变成0，这就给数据的压缩存储带来了便利。也就是说，量化表里的数越大，损失越大，但压缩效果越好，可以根据损失程度的需要对量化表的数值大小进行调整。

3. 基于Huffman编码的压缩

现在我们得到了一堆8×8整数矩阵，并且其中有相当多的0，现在要对它们进行压缩。不过在此之前，由于直流系数数值较大，我们可以继续利用图像的连续性，对相邻矩阵的直流系数进行差分，这样就大大缩小了直流系数的值。这个差分的顺序比较奇怪，我们在后面的实现部分再提。

首先，我们需要把二维的矩阵转为一维的整数序列，在这里我们使用Z字形扫描对矩阵进行遍历，这样能够充分的利用连续性（其实目的是为了让连续的0尽可能多），如下图所示：

此时，序列的第一个数是（经过差分后的）直流系数，我们先不管它。剩下的63个交流系数，对于每个非0值y，假设它前面有x个0，那么我们用二元组(x,y)表示，这就是典型的游程长度编码。

那么你可能就会问了，我原来存一个y只需要一个字节，现在存这个二元组就需要两个字节，这不是一种浪费吗？并非如此，因为JPEG的压缩方式是变长编码，而非一般意义上的定长编码（例如byte,short,int这样的长度固定的整型），所以，真正的压缩数据并非一串字节序列，而是一串二进制位序列！

对于y，我们求出它的二进制长度len。但是，y有可能是负的，在这里JPEG使用一种特殊规则存储变长二进制整数，如下表所示：

这样，原来的二元组就变成了三元组(x,len,y)，实际上我们把x和len视作一个元素，它还是一个二元组(x/len,y)。在实际存储的时候，y直接存长为len的二进制序列，而x/len则使用Huffman编码表示。一般的，我们需要4棵Huffman树，分别对应Y、UV的直流（直流其实就是序列的第一个元素，不带x的(len,y)）和交流部分。这样，给定一个二进制序列，我们就能还原出所有的三元组（Huffman码在前，y在后）了。

注意，如果这个序列以0结尾，那么只需要在最后一个非0数后插一个三元组(0/0,0)表示后面直到结尾都是0。注意：如果序列结尾不是0，那么不需要非得以这个标记结尾。（这个把我坑惨了qaq……）

所以，JPEG中，Huffman编码并不是用来压缩实际的数的，而是用来压缩数的长度信息的。毕竟，数可能不少，但是我们已经让它们缩的很小，长度值也就那几个，再者，这其中还蕴含了中间被压掉的0。

JPEG图像存储格式

下面我们来考虑更加实际的问题，这些东西具体是如何存在JPG图像文件里的？我们就来看看JPG文件的字节码里面都有些什么东西。

FBI Warning：由于种种历史原因，JPEG图像文件中的所有整数全部使用大端方式存储，也就是说整数的高位字节在地址较低的位置，这与平常x86平台上使用的小端方式相反。

JPEG图像一般使用JIFF格式进行存储，它会在文件中划分出若干个段，每段的开头都是“0xFF+段类型字节码”的双字节标识，并且，0xFF是一个很特殊的字节，它一般是用来标识段的，如果它出现在实际数据里，后面需要加一个0字节表示这个0xFF是实际数据的一部分。

以下是一些主要的段，不重要的段我直接忽略了。

1. 文件起始（SOI）

2. 应用程序保留标记0（APP0）

3. 量化表（DQT）

注意：一个DQT段可以放多个量化表，或者多个量化表也可以放在多个DQT段里（Windows画图貌似就这么干的）。

最后两个字段可以多次出现。

4. 帧图像头（SOFO）

5. Huffman表（DHT）

这个和量化表一样，一个DHT可以有多个Huffman表，或者也可以有多个DHT段。

在这里我不得不说一下JPEG中诡异的Huffman存储方式，通过上面的两个字段我们知道了每个字节的Huffman长度，然后我需要这样构造编码：第一个元素是全0，由上一个元素转移到下一个元素时，Huffman码要+1，如果长度不够就在最后补0。这样造出来的Huffman树是满足左子树的叶子深度不超过右子树的。

那么为什么Huffman编码的都是字节呢？是这样的，直流的len​当然是一个字节，交流的x/len​要把x​放到高8位，如果x>15​，那么我们需要在前面放(15/0,0)表示连续的16个0​，这样所有的x​都是不超过0xF的。

6. 扫描头（SOS）

7. 图像压缩数据

紧接在SOS段之后，这一大串二进制序列从第一个字节的最高位开始一直到最后一个字节。

在这个序列中，被压缩后的矩阵的排列方式是这样的，首先将图像划分为若干16×16​子矩阵，每个子矩阵里有4个Y分量的8×8子矩阵，U和V分量的8×8子矩阵各一个，如下图所示：

我们在最终编码时，把这些16×16子矩阵视作一个基本单位，基本单位的顺序就是一行一行的来，而基本单位内部是6个8×8子矩阵的编码，它们的顺序是Y_1,Y_2,Y_3,Y_4,U,V，注意：直流系数差分的顺序就是相同分量的子矩阵在最终的矩阵序列中的出现顺序！

8. 文件结尾（EOI）

以下是我花了4天时间写的一个JPEG编码/解码示例程序，写起来真的比较恶心，而且还存在着两个我无法解决或解释的问题：

①不知道为什么，我在解码最后需要对颜色值进行一个r+=295; g-=18; b+=343;的变换才能显示出正常图像，而且这3个神秘的常数对于所有图像都是成立的，编码时也需要在最开始倒着来一遍。

②大概是我对于颜色变换的中间某一步的理解出现了些许偏差，我编码后的jpg图像虽然跟原图差不多，但在正常的图像查看软件上会存在网格纹的干扰，但是我用自己的解码器显示就十分正常，而且解码器显示画图保存的jpg也会在颜色细节上有一些不好的地方，应该是编码器和解码器在相同的地方都出了点偏差。

code：

COLORREF cols[N][N]; //用来保存位图的颜色值，属于Encode的输入，Decode的输出
const double PI=acos(-1);
struct Matrix{
    double v[8][8];
    int n;
    Matrix(int _n) {  memset(v,0,sizeof(v)); n=_n;  }
    friend Matrix operator *(const Matrix &a,const Matrix &b){
        Matrix c(a.n);
        for(int i=0;i<a.n;i++)
            for(int k=0;k<a.n;k++)
                for(int j=0;j<a.n;j++) c.v[i][j]+=a.v[i][k]*b.v[k][j];
        return c;
    }
    Matrix Inv(){
        double me[8][8],res[8][8];
        memset(res,0,sizeof(res)); memcpy(me,v,sizeof(v));
        for(int i=0;i<n;i++) res[i][i]=1;
        for(int i=0;i<n;i++){
            int lst=-1;
            for(int j=i;j<n;j++) if(fabs(me[j][i])>1e-7) lst=j;
            for(int j=0;j<n;j++) swap(me[i][j],me[lst][j]),swap(res[i][j],res[lst][j]);
            double dk=me[i][i];
            for(int j=0;j<n;j++) me[i][j]/=dk,res[i][j]/=dk;
            for(int j=0;j<n;j++){
                if(j==i) continue;
                double x=me[j][i];
                for(int k=0;k<n;k++) me[j][k]-=x*me[i][k],res[j][k]-=x*res[i][k];
            }
        }
        Matrix c(n); memcpy(c.v,res,sizeof(res)); return c;
    }
    static Matrix DCT(){
        Matrix dct(8);
        for(int i=0;i<8;i++)
            for(int j=0;j<8;j++){
                if(i==0) dct.v[i][j]=1/sqrt(8);
                else dct.v[i][j]=cos(PI*(2*j+1)*i/16)/2;
            }
        return dct;
    }
    static Matrix YUV(){
        Matrix yuv(3);
        yuv.v[0][0]=0.299; yuv.v[0][1]=0.578; yuv.v[0][2]=0.114;
        yuv.v[1][0]=-0.1687; yuv.v[1][1]=-0.3313; yuv.v[1][2]=0.5;
        yuv.v[2][0]=0.5; yuv.v[2][1]=-0.4187; yuv.v[2][2]=-0.0813;
        return yuv;
    }
    Matrix T(){
        Matrix c(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++) c.v[i][j]=v[j][i];
        return c;
    }
};
int realheight,realwidth; //图像的宽度和高度，属于Encode的输入和Decode的输出
unsigned char buf[10000000];
double val[3][N][N];
struct BitIterator{
    unsigned char **p;
    int *pos;
    BitIterator(unsigned char *ptr){
        p=new unsigned char*; *p=ptr;
        pos=new int; *pos=7;
    }
    int Next(){
        int x=((*(*p))>>(*pos))&1;
        if(*pos==0){
            *pos=7;
            if(*(*p)==0xff) (*p)++;
            (*p)++;
        }else (*pos)--;
        return x;
    }
    int GetVal(int len){
        if(!len) return 0;
        int x=0; for(int i=0;i<len;i++) x=(x<<1)|Next();
        if(x>=(1<<(len-1))) return x;
        return x-(1<<len)+1;
    }
    void WriteByte(unsigned char byte){
        if(*pos!=7) (*p)++,*pos=7;
        *(*p)=byte; (*p)++;
    }
    void Out(int bit){
        if(*pos==7) *(*p)=0;
        *(*p)|=(bit<<(*pos));
        if(*pos==0){
            *pos=7;
            if(*(*p)==0xff) (*p)++,*(*p)=0;
            (*p)++;
        }else (*pos)--;
    }
    void OutVal(int x,int len){
        if(x<0) x=(1<<len)+x-1;
        for(int i=len-1;i>=0;i--) Out((x>>i)&1);
    }
};
int calc(int x){
    if(x==0) return 1;
    int ans=0; while(x) x>>=1,ans++; return ans;
}
struct Huffman{
    vector<int> ch[2],val,vec;
    int cnt[17];
    map<int,int> mp_len,mp_code;
    Huffman(){
        ch[0].clear(); ch[1].clear(); val.clear();
        ch[0].push_back(0); ch[1].push_back(0); val.push_back(0);
    }
    typedef struct _n{
        int cnt; vector<int> vec;
        friend bool operator <(_n a,_n b) { return a.cnt>b.cnt; }
    }node;
    void dfs(int pt,int deep,int x){
        if(!ch[0][pt]&&!ch[1][pt]) mp_len[val[pt]]=deep,mp_code[val[pt]]=x; else{
            if(ch[0][pt]) dfs(ch[0][pt],deep+1,x<<1);
            if(ch[1][pt]) dfs(ch[1][pt],deep+1,(x<<1)|1);
        }
    }
    void Build(map<int,int> &mp){
         for(auto i:mp) vec.push_back(i.first);
         memset(cnt,0,sizeof(cnt));
         if(mp.size()==12){
             cnt[1]=0; cnt[2]=1; cnt[3]=5; cnt[4]=1; cnt[5]=1; cnt[6]=1; cnt[7]=1;
             cnt[8]=1; cnt[9]=1;
         }else{
             cnt[2]=2; cnt[3]=1; cnt[4]=3; cnt[5]=3; cnt[6]=2; cnt[7]=4; cnt[8]=3;
             cnt[9]=5; cnt[10]=5; cnt[11]=4; cnt[12]=4; cnt[15]=1; cnt[16]=0x7d;
         }
         sort(vec.begin(),vec.end(),[&mp](int a,int b)->bool{
             return mp[a]>mp[b];
             });
         int lst=-1,lstlen=-1;
         int ptr=1,tot=0;
         for(int i:vec){
             while(tot==cnt[ptr]) tot=0,ptr++;
             if(lst==-1) lst=0,lstlen=ptr; else{
                 lst++; lst<<=(ptr-max(calc(lst),lstlen)); lstlen=ptr;
             }
             wstring str;
             for(int j=lstlen-1;j>=0;j--) str.push_back(L'0'+((lst>>j)&1));
             Insert(str,i);
             tot++;
         }
         dfs(0,0,0);
    }
    void Insert(wstring str,int x){
        int cur=0;
        for(int i=0;i<str.length();i++){
            if(!ch[str[i]-'0'][cur]) ch[str[i]-'0'][cur]=val.size(),ch[0].push_back(0),ch[1].push_back(0),val.push_back(0);
            cur=ch[str[i]-'0'][cur];
        }
        val[cur]=x;
    }
    int Read(BitIterator it){
        int cur=0; 
        while(ch[0][cur]||ch[1][cur]){
            cur=ch[it.Next()][cur];
        }
        return val[cur];
    }
};
unsigned short get16(unsigned char *p){
    unsigned short x=*p; return (x<<8)|(*(p+1));
}
void write16(unsigned char *p,unsigned short x){
    *p=((x>>8)&0xff); *(p+1)=x&0xff;
}
void Move(int &x,int &y){
    if((x+y)&1){ // down left
        if(x-1>=0&&y+1<8) x--,y++; else{
            if(x==0) y++; else x++;
        }
    }else{ // up right
        if(x+1<8&&y-1>=0) x++,y--; else{
            if(x==7) y++; else x++;
        }
    }
}
Matrix GetData(BitIterator it,Huffman &huf_dc,Huffman &huf_ac){
    Matrix res(8);
    int x=0,y=0;
    res.v[0][0]=it.GetVal(huf_dc.Read(it));
    Move(x,y);
    int ptr=1; while(ptr<64){
        int head=huf_ac.Read(it);
        if(head==0){
            while(ptr<64) res.v[y][x]=0,Move(x,y),ptr++;
            break;
        }
        int len0=(head>>4)&0xf,lenx=head&0xf;
        for(int i=0;i<len0;i++) res.v[y][x]=0,Move(x,y),ptr++;
        res.v[y][x]=it.GetVal(lenx); Move(x,y); ptr++;
    }
    return res;
}
void Decode(unsigned char *ptr){ //change realheight,realwidth,cols
    Matrix dct=Matrix::DCT().Inv(),dctT=dct.T();
    Matrix yuv=Matrix::YUV().Inv();
    Huffman ht[2][2]; //[dc/ac][id]
    int tab[4][8][8]; memset(tab,0,sizeof(tab));
    int used_tab[4],id[4],x_sample[4],y_sample[4];
    ptr+=2; while(1){
        if(*(ptr+1)==0xDB){ // div table
            unsigned char *dest=ptr+2+get16(ptr+2);
            ptr+=4;
            while(ptr<dest){
                unsigned char x=*ptr; ptr++;
                int len=1+((x>>4)&0xf),id=x&0xf; 
                for(int i=0;i<8;i++)
                    for(int j=0;j<8;j++){
                        if(len==1) tab[id][i][j]=*ptr; else tab[id][i][j]=get16(ptr);
                        ptr+=len;
                    }
            } continue;
        }
        if(*(ptr+1)==0xc0){ // info header
            realheight=get16(ptr+5);
            realwidth=get16(ptr+7);
            ptr+=10;
            for(int i=0;i<3;i++,ptr+=3){
                id[i]=*ptr; x_sample[i]=((*(ptr+1))>>4)&0xf;
                y_sample[i]=(*(ptr+1))&0xf; used_tab[i]=*(ptr+2);
            } 
            continue;
        }
        if(*(ptr+1)==0xc4){ // Huffman Table
            unsigned char *dest=ptr+2+get16(ptr+2);
            ptr+=4;
            while(ptr<dest){
                int cnt[17]; memset(cnt,0,sizeof(cnt));
                int dc_ac=((*ptr)>>4)&0xf,id=(*ptr)&0xf; ptr++;
                for(int i=1;i<=16;i++,ptr++) cnt[i]=*ptr;
                int lst=-1,lstlen=-1;
                for(int i=1;i<=16;i++){
                    while(cnt[i]){
                        if(lst==-1) lstlen=i,lst=0; else{
                            lst++; if(i>lstlen) lst<<=(i-max(lstlen,calc(lst))),lstlen=i;
                        }
                        wstring str; int x=lst;
                        if(!x) str=L"0";
                        else while(x) str.push_back(x%2+L'0'),x>>=1;
                        while(str.length()<i) str.push_back(L'0');
                        reverse(str.begin(),str.end());
                        ht[dc_ac][id].Insert(str,*ptr);
                        ptr++; cnt[i]--;
                    }
                }
            } 
            continue;
        }
        if(*(ptr+1)==0xda){
            int dc_id[4],ac_id[4];
            ptr+=5;
            for(int i=0;i<3;i++,ptr+=2){
                dc_id[i]=((*(ptr+1))>>4)&0xf;
                ac_id[i]=(*(ptr+1))&0xf;
            }
            ptr+=3;
            BitIterator it(ptr);
            int height=(realheight/16+(realheight%16!=0))*16;
            int width=(realwidth/16+(realwidth%16!=0))*16;
            double dCr=0,dCv=0,dY=0;
            for(int y=0;y<height;y+=16){
                for(int x=0;x<width;x+=16){
                    Matrix Y1=GetData(it,ht[0][dc_id[0]],ht[1][ac_id[0]]);
                    dY+=Y1.v[0][0]; Y1.v[0][0]=dY;
                    Matrix Y2=GetData(it,ht[0][dc_id[0]],ht[1][ac_id[0]]);
                    dY+=Y2.v[0][0]; Y2.v[0][0]=dY;
                    Matrix Y3=GetData(it,ht[0][dc_id[0]],ht[1][ac_id[0]]);
                    dY+=Y3.v[0][0]; Y3.v[0][0]=dY;
                    Matrix Y4=GetData(it,ht[0][dc_id[0]],ht[1][ac_id[0]]);
                    dY+=Y4.v[0][0]; Y4.v[0][0]=dY;
                    Matrix Cr=GetData(it,ht[0][dc_id[1]],ht[1][ac_id[1]]);
                    Matrix Cv=GetData(it,ht[0][dc_id[2]],ht[1][ac_id[2]]);
                    dCr+=Cr.v[0][0]; dCv+=Cv.v[0][0];
                    Cr.v[0][0]=dCr; Cv.v[0][0]=dCv;
                    for(int i=0;i<8;i++){
                        for(int j=0;j<8;j++) 
                            Cr.v[i][j]*=tab[used_tab[1]][i][j],Cv.v[i][j]*=tab[used_tab[2]][i][j];
                    }
                    Cr=dct*Cr*dctT; Cv=dct*Cv*dctT;
                    for(int i=0;i<8;i++){
                        for(int j=0;j<8;j++){
                            val[0][y+i][x+j]=Y1.v[i][j];
                            val[0][y+8+i][x+j]=Y3.v[i][j];
                            val[0][y+i][x+8+j]=Y2.v[i][j];
                            val[0][y+8+i][x+8+j]=Y4.v[i][j];
                            for(int a=0;a<2;a++){
                                for(int b=0;b<2;b++){
                                    val[1][y+i*2+a][x+j*2+b]=Cr.v[i][j];
                                    val[2][y+i*2+a][x+j*2+b]=Cv.v[i][j];
                                }
                            }
                        }
                    }
                }
            } 
            for(int y=0;y<height;y+=8){
                for(int x=0;x<width;x+=8){
                    Matrix Y(8);
                    for(int i=0;i<8;i++)
                        for(int j=0;j<8;j++) Y.v[i][j]=val[0][y+i][x+j]*tab[used_tab[0]][i][j];
                    Y=dct*Y*dctT;
                    for(int i=0;i<8;i++)
                        for(int j=0;j<8;j++) val[0][y+i][x+j]=Y.v[i][j];
                }
            }
            for(int i=0;i<realheight;i++){
                for(int j=0;j<realwidth;j++){
                    val[1][i][j]-=128; val[2][i][j]-=128;
                    int r=yuv.v[0][0]*val[0][i][j]+yuv.v[0][1]*val[1][i][j]+yuv.v[0][2]*val[2][i][j];
                    int g=yuv.v[1][0]*val[0][i][j]+yuv.v[1][1]*val[1][i][j]+yuv.v[1][2]*val[2][i][j];
                    int b=yuv.v[2][0]*val[0][i][j]+yuv.v[2][1]*val[1][i][j]+yuv.v[2][2]*val[2][i][j];
                    r+=295; g-=18; b+=343;
                    r=min(255,max(0,r)); g=min(255,max(0,g)); b=min(255,max(0,b));
                    cols[i][j]=RGB(r,g,b);
                }
            }
            return;
        }
        ptr+=(2+get16(ptr+2));
    }
}
struct IntMatrix{
    int v[8][8];
};
IntMatrix Ys[N/8][N/8],Us[N/16][N/16],Vs[N/16][N/16];
struct Node{
    unsigned char head;
    int x,len;
    bool isac;
    bool htid;
};
int D(int &source,int dx){
    int rest=source; source-=dx; return rest;
}
int GetLen(int x){
    if(x==0) return 0; if(x<0) x=-x;
    int ans=0; while(x) ans++,x>>=1; return ans;
}
void PushMatrix(IntMatrix &mat,vector<Node> &flow,int htid){
    Node dc; dc.head=dc.len=GetLen(mat.v[0][0]); dc.isac=0; dc.x=mat.v[0][0];
    dc.htid=htid; flow.push_back(dc);
    vector<int> vec;
    int x=0,y=0; Move(x,y);
    for(int i=1;i<64;i++) vec.push_back(mat.v[y][x]),Move(x,y);
    for(int i=0;i<vec.size();i++){
        int j=i; while(j+1<vec.size()&&!vec[j]) j++;
        if(j+1==vec.size()&&!vec[j]){
            Node ac; ac.head=ac.len=ac.x=0; ac.isac=1; ac.htid=htid; flow.push_back(ac); return;
        }
        j=min(j,i+15);
        Node ac; ac.head=((j-i)<<4)|GetLen(vec[j]); ac.len=GetLen(vec[j]);
        ac.isac=1; ac.htid=htid; ac.x=vec[j]; flow.push_back(ac);
        i=j;
    }
}
int Encode(unsigned char *ptr){ // use realheight,realwidth,cols[][], write jpeg flow in ptr
    unsigned char *begin=ptr;
    int table[2][8][8]={
        {
        {2,1,1,2,1,1,2,2},
        {2,2,2,2,2,2,3,5},
        {3,3,3,3,3,6,4,4},
        {3,5,7,6,7,7,7,6},
        {7,7,8,9,0xb,9,8,8},
        {0xa,8,7,7,0xa,0xd,0xa,0xa},
        {0xb,0xc,0xc,0xc,0xc,7,9,0xe},
        {0xf,0xd,0xc,0xe,0xb,0xc,0xc,0xc}
        },
        {
        {2,2,2,3,3,3,6,3},
        {3,6,0xc,8,7,8,0xc,0xc},
        {0xc,0xc,0xc,0xc,0xc,0xc,0xc,0xc},
        {0xc,0xc,0xc,0xc,0xc,0xc,0xc,0xc},
        {0xc,0xc,0xc,0xc,0xc,0xc,0xc,0xc},
        {0xc,0xc,0xc,0xc,0xc,0xc,0xc,0xc},
        {0xc,0xc,0xc,0xc,0xc,0xc,0xc,0xc},
        {0xc,0xc,0xc,0xc,0xc,0xc,0xc,0xc}
        }   
    };
    *ptr=0xff; ptr++; *ptr=0xd8; ptr++;
    *ptr=0xff; ptr++; *ptr=0xe0; ptr++; //app0 block
    write16(ptr,16); ptr+=2;
    *ptr=0x4a; ptr++; *ptr=0x46; ptr++; *ptr=0x49; ptr++; *ptr=0x46; ptr++; *ptr=0; ptr++; //JFIF\0
    *ptr=1; ptr++; *ptr=2; ptr++; //ver
    *ptr=1; ptr++; //1 pixel/inch
    write16(ptr,120); ptr+=2; write16(ptr,120); ptr+=2; //120dpi
    *ptr=0; ptr++; *ptr=0; ptr++; //no small image
    for(int op=0;op<2;op++){
        *ptr=0xff; ptr++; *ptr=0xdb; ptr++;
        write16(ptr,67); ptr+=2;
        *ptr=op; ptr++;
        for(int i=0;i<8;i++)
            for(int j=0;j<8;j++) *ptr=table[op][i][j],ptr++;
    }
    *ptr=0xff; ptr++; *ptr=0xc0; ptr++; //SOFO block
    write16(ptr,17); ptr+=2;
    *ptr=8; ptr++;
    write16(ptr,realheight); ptr+=2; write16(ptr,realwidth); ptr+=2;
    *ptr=3; ptr++; //YUV
    *ptr=1; ptr++; *ptr=0x22; ptr++; *ptr=0; ptr++; // Y comp
    *ptr=2; ptr++; *ptr=0x11; ptr++; *ptr=1; ptr++; // U comp
    *ptr=3; ptr++; *ptr=0x11; ptr++; *ptr=1; ptr++; // V comp
    vector<Node> vec;
    int height=(realheight/16+(realheight%16!=0))*16,width=(realwidth/16+(realwidth%16!=0))*16;
    Matrix yuv=Matrix::YUV(),dct=Matrix::DCT(),dctT=dct.T();
    for(int y=0;y<height;y+=8){
        for(int x=0;x<width;x+=8){
            Matrix me(8);
            for(int i=0;i<8;i++){
                for(int j=0;j<8;j++){
                    int r=cols[y+i][x+j]&255,g=(cols[y+i][x+j]>>8)&255,b=(cols[y+i][x+j]>>16)&255;
                    r-=295; g+=18; b-=343; //?
                    me.v[i][j]=r*yuv.v[0][0]+g*yuv.v[0][1]+b*yuv.v[0][2];
                }
            }
            me=dct*me*dctT;
            for(int i=0;i<8;i++)
                for(int j=0;j<8;j++) Ys[y/8][x/8].v[i][j]=round(me.v[i][j]/table[0][i][j]);
        }
    }
    int dY=0,dU=0,dV=0;
    for(int y=0;y<height;y+=16){
        for(int x=0;x<width;x+=16){
            Matrix rU(8),rV(8);
            for(int i=0;i<16;i++){
                for(int j=0;j<16;j++){
                    int r=cols[y+i][x+j]&255,g=(cols[y+i][x+j]>>8)&255,b=(cols[y+i][x+j]>>16)&255;
                    r-=295; g+=18; b-=343; //?
                    rU.v[i/2][j/2]+=(128+r*yuv.v[1][0]+g*yuv.v[1][1]+b*yuv.v[1][2]);
                    rV.v[i/2][j/2]+=(128+r*yuv.v[2][0]+g*yuv.v[2][1]+b*yuv.v[2][2]);
                }
            }
            for(int i=0;i<8;i++)
                for(int j=0;j<8;j++) rU.v[i][j]/=4,rV.v[i][j]/=4;
            rU=dct*rU*dctT; rV=dct*rV*dctT;
            for(int i=0;i<8;i++)
                for(int j=0;j<8;j++){
                    Us[y/16][x/16].v[i][j]=round(rU.v[i][j]/table[1][i][j]);
                    Vs[y/16][x/16].v[i][j]=round(rV.v[i][j]/table[1][i][j]);
                }
            dY=D(Ys[y/8][x/8].v[0][0],dY);
            dY=D(Ys[y/8][x/8+1].v[0][0],dY);
            dY=D(Ys[y/8+1][x/8].v[0][0],dY);
            dY=D(Ys[y/8+1][x/8+1].v[0][0],dY);
            dU=D(Us[y/16][x/16].v[0][0],dU);
            dV=D(Vs[y/16][x/16].v[0][0],dV);
            PushMatrix(Ys[y/8][x/8],vec,0); PushMatrix(Ys[y/8][x/8+1],vec,0);
            PushMatrix(Ys[y/8+1][x/8],vec,0); PushMatrix(Ys[y/8+1][x/8+1],vec,0);
            PushMatrix(Us[y/16][x/16],vec,1); PushMatrix(Vs[y/16][x/16],vec,1);
        }
    }
    map<int,int> cnt[2][2];
    for(int i=0;i<=11;i++) cnt[0][0][i]=cnt[0][1][i]=0;
    for(int i=0;i<=0xf;i++)
        for(int j=1;j<=0xa;j++) cnt[1][0][(i<<4)|j]=cnt[1][1][(i<<4)|j]=0;
    cnt[1][0][0]=cnt[1][1][0]=0;
    cnt[1][0][0xf0]=cnt[1][1][0xf0]=0;
    Huffman ht[2][2];
    for(int i=0;i<vec.size();i++){
        cnt[vec[i].isac][vec[i].htid][vec[i].head]++;
    }
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++) ht[i][j].Build(cnt[i][j]);
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            *ptr=0xff; ptr++; *ptr=0xc4; ptr++;
            unsigned char *last=ptr;
            ptr+=2;
            *ptr=(i<<4)|j; ptr++;
            for(int k=1;k<=16;k++) *ptr=ht[i][j].cnt[k],ptr++;
            for(int k=0;k<ht[i][j].vec.size();k++) *ptr=ht[i][j].vec[k],ptr++;
            write16(last,ptr-last);
        }
    }
    *ptr=0xff; ptr++; *ptr=0xda; ptr++; // SOS block
    write16(ptr,12); ptr+=2;
    *ptr=3; ptr++;
    *ptr=1; ptr++; *ptr=0; ptr++;
    *ptr=2; ptr++; *ptr=0x11; ptr++;
    *ptr=3; ptr++; *ptr=0x11; ptr++;
    *ptr=0; ptr++; *ptr=0x3f; ptr++; *ptr=0; ptr++;
    BitIterator it(ptr);
    for(int i=0;i<vec.size();i++){
        int len=ht[vec[i].isac][vec[i].htid].mp_len[vec[i].head];
        int code=ht[vec[i].isac][vec[i].htid].mp_code[vec[i].head];
        for(int j=len-1;j>=0;j--) it.Out((code>>j)&1);
        it.OutVal(vec[i].x,vec[i].len);
    }
    it.WriteByte(0xff); it.WriteByte(0xd9);
    return (*it.p)-begin;
}