双三次插值及优化Word版

1.数学模型

对于一个目的像素，其坐标通过反向变换得到的在原图中的浮点坐标为(i+u,j+v)，其中i、j均为非负整数，u、v为[0,1)区间的浮点数，双三次插值考虑一个浮点坐标(i+u,j+v)周围的16个邻点，目的像素值f(i+u,j+v)可由如下插值公式得到：

f(i+u,j+v) = [A] * [B] * [C]

[A]=[ S(u + 1)S(u + 0)S(u - 1)S(u - 2) ]

┏f(i-1, j-1)f(i-1, j+0)f(i-1, j+1)f(i-1, j+2) ┓

[B]=┃f(i+0, j-1)f(i+0, j+0)f(i+0, j+1)f(i+0, j+2) ┃

┃f(i+1, j-1)f(i+1, j+0)f(i+1, j+1)f(i+1, j+2) ┃

┗f(i+2, j-1)f(i+2, j+0)f(i+2, j+1)f(i+2, j+2) ┛

┏S(v + 1) ┓

[C]=┃S(v + 0) ┃

┃S(v - 1) ┃

┗S(v - 2) ┛

┏1-2*Abs(x)^2+Abs(x)^3, 0<=Abs(x)<1

S(x)=｛4-8*Abs(x)+5*Abs(x)^2-Abs(x)^3, 1<=Abs(x)<2

┗0, Abs(x)>=2

S(x)是对Sin(x*Pi)/x 的逼近（Pi是圆周率——π），为插值核。

2.计算流程

1. 获取16个点的坐标P1、P2……P16

2. 由插值核计算公式S(x) 分别计算出x、y方向的插值核向量Su、Sv

3. 进行矩阵运算，得到插值结果

iTemp1 = Su0 * P1 + Su1 * P5 + Su2 * P9 + Su3 * P13

iTemp2 = Su0 * P2 + Su1 * P6 + Su2 * P10 + Su3 * P14

iTemp3 = Su0 * P3 + Su1 * P7 + Su2 * P11 + Su3 * P15

iTemp4 = Su0 * P4 + Su1 * P8 + Su2 * P12 + Su3 * P16

iResult = Sv1 * iTemp1 + Sv2 * iTemp2 + Sv3 * iTemp3 + Sv4 * iTemp4

4. 在得到插值结果图后，我们发现图像中有“毛刺”，因此对插值结果做了个后处理，即：设该点在原图中的像素值为pSrc，若abs(iResult - pSrc) 大于某阈值，我们认为插值后的点可能污染原图，因此用原像素值pSrc代替。

3. 算法优化

由于双三次插值计算一个点的坐标需要其周围16个点，更有多达20次的乘法及15次的加法，计算量可以说是非常大，势必要进行优化。

我们选择了Intel的SSE2优化技术，它只支持在P4及以上的机器。测试当前CPU是否支持SSE2，可由CPUID指令得到，代码为：

BOOL g_bSSE2 = FALSE;

__asm

{

mov eax, 1;

cpuid;

test edx, 0x04000000;

jz NotSupport;

mov g_bSSE2, 1

NotSupport:

}

支持SSE2的CPU引入了8个128位的寄存器，这样一个寄存器中就可以存放4个点(RGB)，有利于并行计算。

详细代码见Transform.cpp中函数Optimize_Bicubic。

优化中遇到的问题：

1. 图像每个点由RGB通道组成，由于1个SSE2寄存器有16个字节，这样读入4个像素点后，要浪费4个字节，同时要花费时间将数据对齐，即由BRGB | RGBR | GBRG | BRGB 对齐成0RGB | 0RGB | 0RGB | 0RGB ;

2. 读16字节数据到寄存器时，由于图像地址不能保证是16字节对齐，因此需用更多时钟周期的MOVDQU指令(6个以上时钟周期)；如能使地址16字节对齐，则可用MOVDQA 指令(1个时钟周期) ;

3. 为了消除除法及浮点运算，对权值放大256倍，这样在计算插值核时，必须用2Bytes 来表示1个系数，而图像数据都是1Byte，这样在对齐做乘法时，要浪费一半的SSE2寄存器的空间，导致运算时间变长；而若降低插值核的精度，使其在1Byte表示范围内时，运算的精度又大为下降；

4. 对各指令的周期以及若干行指令是否能够并行流水缺乏经验和认识。

附：SSE2指令整理

算术(Arithmetic)指令：

ADDPD--Packed Double-Precision Floating-Point Add SSE2 2个double对应相加

ADDPD xmm0, xmm1/m128

ADDPS--Packed Single-Precision Floating-Point Add SSE 4个float对应相加

ADDPS xmm0, xmm1/m128

ADDSD--Scalar Double-Precision Floating-Point Add

1个double(低端)对应相加SSE2 ADDSD xmm0, xmm1/m64

ADDSS--Scalar Single-Precision Floating-Point Add SSE 1个float(低端)对应相加

ADDSS xmm0, xmm1/m32

PADDQ--Packed Quadword Add

PSADBW--Packed Sum of Absolute Differences

Opcode Instruction Description

0F F6 /r

PSADBW mm1,

mm2/m64

Absolute difference of packed unsigned byte integers

from mm2 /m64 and mm1; differences are then summed

to produce an unsigned word integer result.

66 0F

F6 /r

PSADBW

xmm1,

xmm2/m128

Absolute difference of packed unsigned byte integers

from xmm2 /m128 and xmm1; the 8 low differences and

8 high differences are then summed separately to

produce two word integer results.

Opcode Instruction Description

0F F8 /r

PSUBB mm,

mm/m64

Subtract packed byte integers in mm/m64 from

packed byte integers in mm.

66 0F F8

/r

PSUBB xmm1,

xmm2/m128

Subtract packed byte integers in xmm2/m128 from

packed byte integers in xmm1.

0F F9 /r

PSUBW mm,

mm/m64

Subtract packed word integers in mm/m64 from

packed word integers in mm.

66 0F F9

/r

PSUBW xmm1,

xmm2/m128

Subtract packed word integers in xmm2/m128 from

packed word integers in xmm1.