图形图像处理－之－任意角度的高质量的快速的图像旋转 下篇 补充话题

HouSisong@GMail.com 2007.06.29
tag:图像旋转,任意角度,图像缩放,速度优化,定点数优化,近邻取样插值,二次线性插值,
三次卷积插值,MipMap链,三次线性插值,MMX/SSE优化,CPU缓存优化,AlphaBlend,颜色混合,并行

摘要: 该文章是《任意角度的高质量的快速的图像旋转》的一些高级补充话题；
给出了一个完整的Alpha混合的插值旋转实现；并尝试将旋转函数并行化，从而在多核电脑上获得更快的速度；

任意角度的高质量的快速的图像旋转 全文 分为:
上篇 纯软件的任意角度的快速旋转
中篇 高质量的旋转
下篇 补充话题

正文：
为了便于讨论，这里只处理32bit的ARGB颜色；
代码使用C++;涉及到汇编优化的时候假定为x86平台;使用的编译器为vc2005;
为了代码的可读性,没有加入异常处理代码;
测试使用的CPU为AMD64x2 4200+;
(基础代码参考《图形图像处理－之－任意角度的高质量的快速的图像旋转》系列前面的文章)

A:完整的Alpha混合的双线性插值旋转实现
《高质量的旋转》中已经涉及到了边界的AlphaBlend的问题，这里顺水推舟的实现一个支持全图片Alpha通道Blend混合的双线性插值旋转函数；
首先给出带完整Alpha通道的源图片:




这张图片是带有8比特Alpha的32比特RGB真彩bmp图片;
带的Alpha通道在工具里可能显示不出来，单独提取出来的图示：




函数实现：

void BilInear_BlendBorder_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
{
unsigned long x0=(x_16>>16);
unsigned long y0=(y_16>>16);

TARGB32 pixel[4];
bool IsInPic;
pixel[0]=Pixels_Bound(pic,x0,y0,IsInPic);
if (!IsInPic) pixel[0].a=0;
pixel[2]=Pixels_Bound(pic,x0,y0+1,IsInPic);
if (!IsInPic) pixel[2].a=0;
pixel[1]=Pixels_Bound(pic,x0+1,y0,IsInPic);
if (!IsInPic) pixel[1].a=0;
pixel[3]=Pixels_Bound(pic,x0+1,y0+1,IsInPic);
if (!IsInPic) pixel[3].a=0;

TPicRegion npic;
npic.pdata =&pixel[0];
npic.byte_width=2*sizeof(TARGB32);
//npic.width =2;
//npic.height =2;
BilInear_Fast_MMX(npic,(unsigned short)x_16,(unsigned short)y_16,result);
}

void PicRotary_BilInear_BlendLine_MMX(TARGB32* pDstLine,long dst_border_x0,long dst_in_x0,long dst_in_x1,long dst_border_x1,
const TPicRegion& SrcPic,long srcx0_16,long srcy0_16,long Ax_16,long Ay_16)
{
long x;
for (x=dst_border_x0;x<dst_in_x0;++x)
{
TARGB32 src_color;
BilInear_BlendBorder_MMX(SrcPic,srcx0_16,srcy0_16,&src_color);
if (src_color.a>0)
pDstLine[x]=AlphaBlend_MMX(pDstLine[x],src_color);
srcx0_16+=Ax_16;
srcy0_16+=Ay_16;
}
for (x=dst_in_x0;x<dst_in_x1;++x)
{
TARGB32 src_color;
BilInear_Fast_MMX(SrcPic,srcx0_16,srcy0_16,&src_color);
if (src_color.a==255)
pDstLine[x]=src_color;
else if (src_color.a>0)
pDstLine[x]=AlphaBlend_MMX(pDstLine[x],src_color);
srcx0_16+=Ax_16;
srcy0_16+=Ay_16;
}
for (x=dst_in_x1;x<dst_border_x1;++x)
{
TARGB32 src_color;
BilInear_BlendBorder_MMX(SrcPic,srcx0_16,srcy0_16,&src_color);
if (src_color.a>0)
pDstLine[x]=AlphaBlend_MMX(pDstLine[x],src_color);
srcx0_16+=Ax_16;
srcy0_16+=Ay_16;
}
asm emms
}

void PicRotaryBlendBilInear_MMX(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y)
{
if ( (fabs(ZoomX*Src.width)<1.0e-4) || (fabs(ZoomY*Src.height)<1.0e-4) ) return; //太小的缩放比例认为已经不可见
double tmprZoomXY=1.0/(ZoomX*ZoomY);
double rZoomX=tmprZoomXY*ZoomY;
double rZoomY=tmprZoomXY*ZoomX;
double sinA,cosA;
SinCos(RotaryAngle,sinA,cosA);
long Ax_16=(long)(rZoomX*cosA*(1<<16));
long Ay_16=(long)(rZoomX*sinA*(1<<16));
long Bx_16=(long)(-rZoomY*sinA*(1<<16));
long By_16=(long)(rZoomY*cosA*(1<<16));
double rx0=Src.width*0.5; //(rx0,ry0)为旋转中心
double ry0=Src.height*0.5;
long Cx_16=(long)((-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0)*(1<<16));
long Cy_16=(long)((-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0)*(1<<16));

TRotaryClipData rcData;
rcData.Ax_16=Ax_16;
rcData.Bx_16=Bx_16;
rcData.Cx_16=Cx_16;
rcData.Ay_16=Ay_16;
rcData.By_16=By_16;
rcData.Cy_16=Cy_16;
rcData.dst_width=Dst.width;
rcData.dst_height=Dst.height;
rcData.src_width=Src.width;
rcData.src_height=Src.height;
if (!rcData.inti_clip(move_x,move_y,1)) return;

TARGB32* pDstLine=Dst.pdata;
((TUInt8*&)pDstLine)+=(Dst.byte_width*rcData.out_dst_down_y);
while (true) //to down
{
long y=rcData.out_dst_down_y;
if (y>=Dst.height) break;
if (y>=0)
{
PicRotary_BilInear_BlendLine_MMX(pDstLine,rcData.out_dst_x0_boder,rcData.out_dst_x0_in,
rcData.out_dst_x1_in,rcData.out_dst_x1_boder,Src,rcData.out_src_x0_16,rcData.out_src_y0_16,Ax_16,Ay_16);
}
if (!rcData.next_clip_line_down()) break;
((TUInt8*&)pDstLine)+=Dst.byte_width;
}

pDstLine=Dst.pdata;
((TUInt8*&)pDstLine)+=(Dst.byte_width*rcData.out_dst_up_y);
while (rcData.next_clip_line_up()) //to up
{
long y=rcData.out_dst_up_y;
if (y<0) break;
((TUInt8*&)pDstLine)-=Dst.byte_width;
if (y<Dst.height)
{
PicRotary_BilInear_BlendLine_MMX(pDstLine,rcData.out_dst_x0_boder,rcData.out_dst_x0_in,
rcData.out_dst_x1_in,rcData.out_dst_x1_boder,Src,rcData.out_src_x0_16,rcData.out_src_y0_16,Ax_16,Ay_16);
}
}
}

效果图:



B:在双核上并行三次卷积插值旋转的一个简单实现
(假设图片旋转绘制到目的图片的中间)
这里利用CWorkThreadPool来并行执行任务;
(参见我的文章《并行计算简介和多核CPU编程Demo》，里面有CWorkThreadPool类的完整源代码)
最容易想到的方案就是分成上下两部分分别调用PicRotaryThreeOrder_MMX，从而并行执行；

struct TRotaryThreeOrder_WorkData
{
const TPicRegion* Dst;
const TPicRegion* Src;
double RotaryAngle;
double ZoomX;
double ZoomY;
double move_x;
double move_y;
};

void RotaryThreeOrder_callback(void* wd)
{
TRotaryThreeOrder_WorkData* WorkData=(TRotaryThreeOrder_WorkData*)wd;
PicRotaryThreeOrder_MMX(*WorkData->Dst,*WorkData->Src,WorkData->RotaryAngle,WorkData->ZoomX,WorkData->ZoomY,WorkData->move_x,WorkData->move_y);
}

void PicRotaryThreeOrder_MMX_parallel2(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y)
{
TRotaryThreeOrder_WorkData work_list[2];
TRotaryThreeOrder_WorkData* pwork_list[2];
for (long i=0;i<2;++i)
{
work_list[i].Src=&Src;
work_list[i].RotaryAngle=RotaryAngle;
work_list[i].ZoomX=ZoomX;
work_list[i].ZoomY=ZoomY;
work_list[i].move_x=move_x;
work_list[i].move_y=move_y;
pwork_list[i]=&work_list[i];
}
TPicRegion dst_up=Dst;
dst_up.height=Dst.height/2;
work_list[0].Dst=&dst_up;
TPicRegion dst_down=Dst;
dst_down.pdata=&Pixels(Dst,0,dst_up.height);
dst_down.height=Dst.height-dst_up.height;
work_list[1].Dst=&dst_down;
work_list[1].move_y=move_y-Dst.height/2;
CWorkThreadPool::work_execute(RotaryThreeOrder_callback,(void**)&pwork_list,2);
}

//注：测试图片都是800*600的图片旋转到1004*1004的图片中心 测试成绩取各个旋转角度的平均速度值
//CPU: AMD64x2 4200+
////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicRotaryThreeOrder_MMX_parallel2 68.1 fps
////////////////////////////////////////////////////////////////////////////////

并行化的实现PicRotaryThreeOrder_MMX_parallel2比PicRotaryThreeOrder_MMX的34.4fps快了98.0%!
在双核CPU上执行速度几乎是单核上的2倍!

C:一个通用的针对任意多核并行的一个简单实现
有了上面的并行基础，我们来实现一个更加通用一些的版本；根据CPU核心数来动态分配任务；
实现方式为直接按照扫描行来分配(但这样处理可能不利于内存的高效访问)，就懒得去估算任务量了:)

void PicRotaryThreeOrder_MMX_part(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y,long part_i,long part_count)
{
if ( (fabs(ZoomX*Src.width)<1.0e-4) || (fabs(ZoomY*Src.height)<1.0e-4) ) return; //太小的缩放比例认为已经不可见
double tmprZoomXY=1.0/(ZoomX*ZoomY);
double rZoomX=tmprZoomXY*ZoomY;
double rZoomY=tmprZoomXY*ZoomX;
double sinA,cosA;
SinCos(RotaryAngle,sinA,cosA);
long Ax_16=(long)(rZoomX*cosA*(1<<16));
long Ay_16=(long)(rZoomX*sinA*(1<<16));
long Bx_16=(long)(-rZoomY*sinA*(1<<16));
long By_16=(long)(rZoomY*cosA*(1<<16));
double rx0=Src.width*0.5; //(rx0,ry0)为旋转中心
double ry0=Src.height*0.5;
long Cx_16=(long)((-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0)*(1<<16));
long Cy_16=(long)((-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0)*(1<<16));

TRotaryClipData rcData;
rcData.Ax_16=Ax_16;
rcData.Bx_16=Bx_16;
rcData.Cx_16=Cx_16;
rcData.Ay_16=Ay_16;
rcData.By_16=By_16;
rcData.Cy_16=Cy_16;
rcData.dst_width=Dst.width;
rcData.dst_height=Dst.height;
rcData.src_width=Src.width;
rcData.src_height=Src.height;
if (!rcData.inti_clip(move_x,move_y,2)) return;

TARGB32* pDstLine=Dst.pdata;
((TUInt8*&)pDstLine)+=(Dst.byte_width*rcData.out_dst_down_y);
long run_part_i=0;
while (true) //to down
{
long y=rcData.out_dst_down_y;
if (y>=Dst.height) break;
if (y>=0)
{
if (run_part_i%part_count==part_i)
PicRotary_ThreeOrder_CopyLine_MMX(pDstLine,rcData.out_dst_x0_boder,rcData.out_dst_x0_in,
rcData.out_dst_x1_in,rcData.out_dst_x1_boder,Src,rcData.out_src_x0_16,rcData.out_src_y0_16,Ax_16,Ay_16);
++run_part_i;
}
if (!rcData.next_clip_line_down()) break;
((TUInt8*&)pDstLine)+=Dst.byte_width;
}

pDstLine=Dst.pdata;
((TUInt8*&)pDstLine)+=(Dst.byte_width*rcData.out_dst_up_y);
while (rcData.next_clip_line_up()) //to up
{
long y=rcData.out_dst_up_y;
if (y<0) break;
((TUInt8*&)pDstLine)-=Dst.byte_width;
if (y<Dst.height)
{
if (run_part_i%part_count==part_i)
PicRotary_ThreeOrder_CopyLine_MMX(pDstLine,rcData.out_dst_x0_boder,rcData.out_dst_x0_in,
rcData.out_dst_x1_in,rcData.out_dst_x1_boder,Src,rcData.out_src_x0_16,rcData.out_src_y0_16,Ax_16,Ay_16);
++run_part_i;
}
}
}

struct TRotaryThreeOrder_part_WorkData
{
const TPicRegion* Dst;
const TPicRegion* Src;
double RotaryAngle;
double ZoomX;
double ZoomY;
double move_x;
double move_y;
long part_i;
long part_count;
};

void RotaryThreeOrder_part_callback(void* wd)
{
TRotaryThreeOrder_part_WorkData* WorkData=(TRotaryThreeOrder_part_WorkData*)wd;
PicRotaryThreeOrder_MMX_part(*WorkData->Dst,*WorkData->Src,WorkData->RotaryAngle,WorkData->ZoomX,WorkData->ZoomY,
WorkData->move_x,WorkData->move_y,WorkData->part_i,WorkData->part_count);
}

void PicRotaryThreeOrder_MMX_parallel(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y)
{
long work_count=CWorkThreadPool::best_work_count();
std::vector<TRotaryThreeOrder_part_WorkData> work_list(work_count);
std::vector<TRotaryThreeOrder_part_WorkData*> pwork_list(work_count);
long i;
for (i=0;i<work_count;++i)
{
work_list[i].Dst=&Dst;
work_list[i].Src=&Src;
work_list[i].RotaryAngle=RotaryAngle;
work_list[i].ZoomX=ZoomX;
work_list[i].ZoomY=ZoomY;
work_list[i].move_x=move_x;
work_list[i].move_y=move_y;
work_list[i].part_i=i;
work_list[i].part_count=work_count;
pwork_list[i]=&work_list[i];
}
CWorkThreadPool::work_execute(RotaryThreeOrder_part_callback,(void**)&pwork_list[0],work_count);
}

//注：测试图片都是800*600的图片旋转到1004*1004的图片中心 测试成绩取各个旋转角度的平均速度值
//CPU: AMD64x2 4200+
////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicRotaryThreeOrder_MMX_parallel 65.2 fps
////////////////////////////////////////////////////////////////////////////////



这个实现能应付大多数时候的并行需求了，包括以后的4核8核...;
要是谁有4核CPU来跑一下，速度应该能达到120多fps .

(注意:这里的并行任务分割方案仅仅是简单的举例(用了代码改动最小的方案)，你应该根据你的需求来更好的并行化你的任务； 如果分割后的单个任务太小，并行的优势可能就体现不出来，甚至于更慢；)