SSE用いたアルファブレンドの高速化3 アルファチャンネルをもつ画像同士のアルファブレンド
Posted on Mon 28 November 2011 in ペイントツール
↓の続き
http://d.hatena.ne.jp/hiroki0/20090914/1252931917
http://d.hatena.ne.jp/hiroki0/20100320/1269074009
使用するアルファブレンド計算式
チャンネルサイズ8bitとして
ここで、aはアルファ値、dは画素値である。
この計算式をコードにすると以下のようになる。
255による除算は、高速化のために256の除算に変更しシフト演算で処理している。
また、計算量を減らすためにアルファ値により分岐を行なっている。
今回は、この処理をSSEにより実装して計測を行う。
if ((over_pix->a == 255) || (under_pix->a == 0)) {
dst_pix->value = over_pix->value;
} else if (over_pix->a == 0) {
dst_pix->value = under_pix->value;
} else if (under_pix->a == 255) {
uint8_t ra = ~over_pix->a;
dst_pix->b = (over_pix->b*over_pix->a + under_pix->b*ra)>>8;
dst_pix->g = (over_pix->g*over_pix->a + under_pix->g*ra)>>8;
dst_pix->r = (over_pix->r*over_pix->a + under_pix->r*ra)>>8;
dst_pix->a = 255;
} else {
uint8_t ra = ~over_pix->a;
uint8_t alpha = 255 - (ra*(255 - under_pix->a)>>8);
double inv_alpha = 1.0/alpha;
dst_pix->b = (over_pix->b*over_pix->a + ((under_pix->b*under_pix->a*ra)>>8))*inv_alpha;
dst_pix->g = (over_pix->g*over_pix->a + ((under_pix->g*under_pix->a*ra)>>8))*inv_alpha;
dst_pix->r = (over_pix->r*over_pix->a + ((under_pix->r*under_pix->a*ra)>>8))*inv_alpha;
dst_pix->a = alpha;
}
計測方法
32bit,64bit向けに2通りに実行ファイルを作成、それぞれについて
- SSEを使わない処理
- _mm_loadu_psでメモリ読み込み、_mm_storeu_psでメモリ書き込み
- _mm_load_psでメモリ読み込み、_mm_store_psでメモリ書き込み
- _mm_load_psでメモリ読み込み、_mm_stream_si128でメモリ書き込み
の4通りの処理の処理速度を計測する。
- チャンネルサイズ8bit 5700x5700BGR画像を2枚読み込み、アルファ値を指定してアルファブレンドを行う。
- 画像の保存、読み込みはopenCV 2.3を用いた。
- アルファブレンド部分の処理のみ処理時間を計測した。
- アライメントは常に合っている状態で処理を行う。
アルファブレンドの処理は、アルファ値により分岐を入れているので
ただコピーを行う部分、アルファ値1つによるアルファブレンドを行う部分、
アルファ値2つでアルファブレンドを行う部分の3箇所の計測を行いたい。
なので、与えるアルファチャンネルは、
- 上下ともにアルファ値255
- 下になる画像全てアルファ値255、上になる画像はアルファ値を均等に分布
- 上下ともにアルファ値を均等に分布
の三通りとした。
計測に使用した環境
- CPU: Core i7 2600 3.40GHz
- OS: windows 7 64bit
- コンパイラ: Visual Studio 2010
- 最適化オプション: 実行速度 (/O2)
計測結果
上下ともにアルファ値255
-
32bit向け実行ファイル
非SSE _mm_loadu_ps + _mm_storeu_ps _mm_load_ps + _mm_store_ps _mm_load_ps + _mm_stream_si128
1回目 111.002 msec 95.006 msec 90.9281 msec 90.853 msec 2回目 87.0272 msec 70.3019 msec 64.3998 msec 77.862 msec 3回目 86.3182 msec 66.67 msec 72.2693 msec 62.0172 msec
-
64bit向け実行ファイル
非SSE _mm_loadu_ps + _mm_storeu_ps _mm_load_ps + _mm_store_ps _mm_load_ps + _mm_stream_si128
1回目 85.2114 msec 82.4653 msec 85.074 msec 83.51 msec 2回目 56.3726 msec 63.0728 msec 71.2231 msec 56.5069 msec 3回目 55.2358 msec 62.6271 msec 62.9312 msec 58.8877 msec
このアルファ値では、分岐処理により上の画像の画素を出力メモリにコピーする処理だけになる。
SSEのロード・ストア命令による違いはバラツキがあり、殆ど無いと思われる。
32bitでは、SSEを使ったほうが1.1倍程度速くなっている。
32bitと64bitの結果を比較して、SSEを使わない処理では64bitの方が早くなっている。
下になる画像全てアルファ値255、上になる画像はアルファ値を均等に分布
-
32bit向け実行ファイル
非SSE _mm_loadu_ps + _mm_storeu_ps _mm_load_ps + _mm_store_ps _mm_load_ps + _mm_stream_si128
1回目 193.727 msec 114.325 msec 94.4707 msec 101.514 msec 2回目 187.955 msec 93.6102 msec 93.4529 msec 91.6304 msec 3回目 190.177 msec 104.997 msec 94.9949 msec 91.032 msec
-
64bit向け実行ファイル
非SSE _mm_loadu_ps + _mm_storeu_ps _mm_load_ps + _mm_store_ps _mm_load_ps + _mm_stream_si128
1回目 167.373 msec 99.5402 msec 89.3116 msec 100.947 msec 2回目 198.225 msec 91.6833 msec 99.7168 msec 89.8282 msec 3回目 186.376 msec 87.6821 msec 98.8608 msec 97.6576 msec
このアルファ値では、分岐処理により主に一つのアルファ値によるアルファブレンドを行う処理になる。
SSEのロード・ストア命令による違いはバラツキがあり、殆ど無いと思われる。
32bitと64bitによる違いは殆ど見られない。
また、SSEにより1.6倍〜2.0倍程度速くなっていることがわかる。
上下ともにアルファ値を均等に分布
-
32bit向け実行ファイル
非SSE _mm_loadu_ps + _mm_storeu_ps _mm_load_ps + _mm_store_ps _mm_load_ps + _mm_stream_si128
1回目 615.202 msec 182.197 msec 208.116 msec 176.766 msec 2回目 576.132 msec 182.21 msec 196.25 msec 172.828 msec 3回目 586.831 msec 176.329 msec 179.473 msec 171.598 msec
-
64bit向け実行ファイル
非SSE _mm_loadu_ps + _mm_storeu_ps _mm_load_ps + _mm_store_ps _mm_load_ps + _mm_stream_si128
1回目 327.867 msec 208.067 msec 204.649 msec 170.796 msec 2回目 328.621 msec 173.355 msec 176.605 msec 160.269 msec 3回目 336.624 msec 167.226 msec 168.419 msec 170.916 msec
このアルファ値では、分岐処理により主に2つのアルファ値によるアルファブレンドを行う処理になる。
SSEのロード・ストア命令による違いはバラツキがあり、殆ど無いと思われる。
非SSEの処理の32bitの方は64bitに対して1.8倍程度遅くなっており、かなり時間がかかっていることがわかる。
SSEによる速度の向上は32bitで約3倍64bitで約1.8倍である。
まとめ
測定条件のいずれの場合もSSEのロード・ストア命令による違いは見られなかった。
今回は、アライメントは常に合っている状態で計測を行ったのでアライメントをずらいた状態では違いが見られるかもしれない。
コード
SSEでの実装方針は、4ピクセル分のデータをロード、SSEのレジスタに16bitで8つのデータを持ち2ピクセルづつ処理し、
4ピクセル分ストアするようにする。
2つのアルファ値をもつアルファブレンドの処理では途中、除算処理が入る。
SSEの命令には整数除算が無いので、[16bit整数]→[32bit整数]→[32bit浮動小数]と変換して
_mm_rcp_ps()で逆数を求めて除算を行う。
_mm_rcp_ps()を用いることで_mm_div_ps()よりも高速に計算することができる。
_mm_rcp_ps()は11bit程度しか精度がないが、最後は8bit整数にするので問題ないはず。
#include <iostream>
#include <opencv/cv.h>
#include <opencv/highgui.h>
#include "IETimer.h"
union UCvPixel32
{
uint32_t value;
struct
{
uint8_t b,g,r,a;
};
};
typedef UCvPixel32 UCvPixel;
inline uint8_t* GetNextLineAddress(const IplImage* image, uint8_t* addr)
{
return addr + image->widthStep;
}
inline uint8_t* GetPixelAddress(const IplImage* image, int x, int y)
{
assert(image);
assert(image->depth == IPL_DEPTH_8U);
assert((0 <= x && x <= image->width-1) && (0 <= y && y <= image->height-1));
return (uint8_t*)(image->imageData + y*image->widthStep + x*image->nChannels);
}
inline UCvPixel32* GetNextLineAddress32(const IplImage* image, UCvPixel32* addr)
{
return (UCvPixel32*)((uint8_t*)addr + image->widthStep);
}
inline UCvPixel32* GetPixelAddress32(const IplImage* image, int x, int y)
{
assert(image);
assert(image->depth == IPL_DEPTH_8U);
assert((0 <= x && x <= image->width-1) && (0 <= y && y <= image->height-1));
return (UCvPixel32*)(image->imageData + y*image->widthStep + x*image->nChannels);
}
int main(void)
{
int x,y;
int r,g,b;
int r1, g1, b1;
int r2, g2, b2;
UCvPixel *over_pix;
UCvPixel *under_pix;
UCvPixel *dst_pix;
UCvPixel *over_line;
UCvPixel *under_line;
UCvPixel *dst_line;
IplImage* test1 = cvLoadImage("test1.bmp");
IplImage* test2 = cvLoadImage("test2.bmp");
IplImage* b_ch = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 1);
IplImage* g_ch = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 1);
IplImage* r_ch = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 1);
IplImage* a_ch1 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 1);
IplImage* a_ch2 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 1);
IplImage* img1 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 4);
IplImage* img2 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 4);
IplImage* img3 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 4);
IplImage* img4 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 4);
IplImage* img5 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 4);
IplImage* img6 = cvCreateImage(cvGetSize(test1), IPL_DEPTH_8U, 4);
//上になる画像のアルファ値をすべて255にする
cvSet(a_ch1, cvScalar(255));
//下になる画像のアルファ値をすべて255にする
cvSet(a_ch2, cvScalar(255));
//上になる画像のアルファ値を均等に分布させる
//uint8_t* line1 = GetPixelAddress(a_ch1, 0, 0);
//for (y = 0; y < a_ch1->height; y++) {
// uint8_t* pix = line1;
// for (x = 0; x < a_ch1->width; x++) {
// (*pix++) = ((double)x / a_ch1->width) * 255;
// }
// line1 = GetNextLineAddress(a_ch1, line1);
//}
//下になる画像のアルファ値を均等に分布させる
//uint8_t* line2 = GetPixelAddress(a_ch2, 0, 0);
//for (y = 0; y < a_ch2->height; y++) {
// uint8_t* pix = line2;
// for (x = 0; x < a_ch2->width; x++) {
// (*pix++) = ((double)y / a_ch2->height) * 255;
// }
// line2 = GetNextLineAddress(a_ch2, line2);
//}
cvSplit(test1, b_ch, g_ch, r_ch, NULL);
cvMerge(b_ch, g_ch, r_ch, a_ch1, img1);
cvSplit(test2, b_ch, g_ch, r_ch, NULL);
cvMerge(b_ch, g_ch, r_ch, a_ch2, img2);
//非SSE
IETimer timer;
over_line = GetPixelAddress32(img1, 0, 0);
under_line = GetPixelAddress32(img2, 0, 0);
dst_line = GetPixelAddress32(img3, 0, 0);
for(y=0; y<img1->height; y++){
over_pix = over_line;
under_pix = under_line;
dst_pix = dst_line;
for(x=0; x<img1->width; x++){
//alpha blend
if ((over_pix->a == 255) || (under_pix->a == 0)) {
dst_pix->value = over_pix->value;
} else if (over_pix->a == 0) {
dst_pix->value = under_pix->value;
} else if (under_pix->a == 255) {
uint8_t ra = ~over_pix->a;
dst_pix->b = (over_pix->b*over_pix->a + under_pix->b*ra)>>8;
dst_pix->g = (over_pix->g*over_pix->a + under_pix->g*ra)>>8;
dst_pix->r = (over_pix->r*over_pix->a + under_pix->r*ra)>>8;
dst_pix->a = 255;
} else {
uint8_t ra = ~over_pix->a;
uint8_t alpha = 255 - (ra*(255 - under_pix->a)>>8);
double inv_alpha = 1.0/alpha;
dst_pix->b = (over_pix->b*over_pix->a + ((under_pix->b*under_pix->a*ra)>>8))*inv_alpha;
dst_pix->g = (over_pix->g*over_pix->a + ((under_pix->g*under_pix->a*ra)>>8))*inv_alpha;
dst_pix->r = (over_pix->r*over_pix->a + ((under_pix->r*under_pix->a*ra)>>8))*inv_alpha;
dst_pix->a = alpha;
}
over_pix++;
under_pix++;
dst_pix++;
}
over_line = GetNextLineAddress32(img1, over_line);
under_line = GetNextLineAddress32(img2, under_line);
dst_line = GetNextLineAddress32(img3, dst_line);
}
std::cout << timer.elapsed() << " msec" << std::endl;
cvSaveImage("test3.bmp", img3);
//SSE: _mm_loadu_ps + _mm_storeu_ps
timer.restart();
{
over_line = GetPixelAddress32(img1, 0, 0);
under_line = GetPixelAddress32(img2, 0, 0);
dst_line = GetPixelAddress32(img4, 0, 0);
__m128i xmzero = _mm_setzero_si128();
__m128i xmff16 = _mm_set1_epi16(0xff);
for(y=0; y<img1->height; y++){
over_pix = over_line;
under_pix = under_line;
dst_pix = dst_line;
for(x=0; x<img1->width; x+=4){
//4ピクセルロード
__m128i under = _mm_castps_si128( _mm_loadu_ps((float*)under_pix) );
__m128i over = _mm_castps_si128( _mm_loadu_ps((float*)over_pix) );
__m128i xmo16, xmu16;
__m128i xmlo, xmhi;
__m128i xmover_alpha, xmover_ralpha, xmunder_alpha, xmunder_ralpha, xmalpha;
__m128 xmf_lo, xmf_hi, xmf_alpha_lo, xmf_alpha_hi;
//load lo double ward
xmo16 = _mm_unpacklo_epi8(over, xmzero);
xmu16 = _mm_unpacklo_epi8(under, xmzero);
//アルファ値を各個読み出し比較し分岐へ
int over_a1 = _mm_extract_epi16(xmo16, 3);
int over_a2 = _mm_extract_epi16(xmo16, 7);
int under_a1 = _mm_extract_epi16(xmu16, 3);
int under_a2 = _mm_extract_epi16(xmu16, 7);
if (((over_a1 == 255 && over_a2 == 255)) || (under_a1 == 0 && under_a2 == 0)) {
xmlo = xmo16;
} else if (under_a1 == 0 && under_a2 == 0) {
xmlo = xmu16;
} else if (under_a1 == 255 && under_a2 == 255) {
//16bitのアルファ値8個をシャッフルをつかってならべ、
//[a1,a1,a1,a1,a2,a2,a2]というようにする。
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha); //255 - over_alpha
//2ピクセル同時にアルファブレンドの計算
xmu16 = _mm_mullo_epi16(xmu16, xmover_ralpha); // under_pix * over_ralpha
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // (under_pix * over_ralpha + over_pix * over_alpha)
xmlo = _mm_srli_epi16(xmo16, 8); // (under_pix * over_ralpha + over_pix * over_alpha) >> 8
//アルファ値255を出力用のレジスタに各個書き込み
xmlo = _mm_insert_epi16(xmlo, 255, 3); //write alpha
xmlo = _mm_insert_epi16(xmlo, 255, 7); //write alpha
} else {
//16bitのアルファ値8個をシャッフルをつかってならべ
//[a1,a1,a1,a1,a2,a2,a2]というようにする。
//上になる画像と下になる画像の2つぶんを行なっておく
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha); // 255 - over_alpha
xmunder_alpha = _mm_shufflelo_epi16(xmu16, _MM_SHUFFLE(3,3,3,3));
xmunder_alpha = _mm_shufflehi_epi16(xmunder_alpha, _MM_SHUFFLE(3,3,3,3));
xmunder_ralpha = _mm_sub_epi16(xmff16, xmunder_alpha); // 255 - under_alpha
//新しくアルファ値になる値の計算
//ここのアルファ値も[a1,a1,a1,a1,a2,a2,a2]となるようにしている
xmalpha = _mm_sub_epi16(xmff16, _mm_srli_epi16( _mm_mullo_epi16(xmover_ralpha, xmunder_ralpha), 8)); //255 - ((over_ralpha * under_ralpha) >> 8)
//16bit幅での2ピクセル同時処理
//ここで 16bit乗算 → 8bit右シフト → 16bit乗算という順で計算しているので誤差が出る。
xmu16 = _mm_mullo_epi16( _mm_srli_epi16( _mm_mullo_epi16(xmunder_alpha, xmover_ralpha), 8), xmu16); //under_pix * ((under_alpha * over_ralpha)>>8)
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // over_pix * over_alpha + (under_pix * under_alpha * over_ralpha)>>8
//32bit浮動小数に変換してピクセル1の除算の計算
xmf_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmo16, xmzero) );
xmf_alpha_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmalpha, xmzero) );
xmf_lo = _mm_mul_ps(xmf_lo, _mm_rcp_ps(xmf_alpha_lo));
////32bit浮動小数に変換してピクセル2の除算の計算
xmf_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmo16, xmzero) );
xmf_alpha_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmalpha, xmzero) );
xmf_hi = _mm_mul_ps(xmf_hi, _mm_rcp_ps(xmf_alpha_hi));
//16bit整数に戻してレジスタに書き込み
xmlo = _mm_packs_epi32(_mm_cvtps_epi32(xmf_lo), _mm_cvtps_epi32(xmf_hi));
//新しくアルファ値になる値を取り出して各個書き込み
int a1 = _mm_extract_epi16(xmalpha, 3);
xmlo = _mm_insert_epi16(xmlo, a1, 3); //write alpha
int a2 = _mm_extract_epi16(xmalpha, 7);
xmlo = _mm_insert_epi16(xmlo, a2, 7); //write alpha
}
//load hi double ward
xmo16 = _mm_unpackhi_epi8(over, xmzero);
xmu16 = _mm_unpackhi_epi8(under, xmzero);
over_a1 = _mm_extract_epi16(xmo16, 3);
over_a2 = _mm_extract_epi16(xmo16, 7);
under_a1 = _mm_extract_epi16(xmu16, 3);
under_a2 = _mm_extract_epi16(xmu16, 7);
if (((over_a1 == 255 && over_a2 == 255)) || (under_a1 == 0 && under_a2 == 0)) {
xmhi = xmo16;
} else if (under_a1 == 0 && under_a2 == 0) {
xmhi = xmu16;
} else if (under_a1 == 255 && under_a2 == 255) {
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmu16 = _mm_mullo_epi16(xmu16, xmover_ralpha); // under_pix * over_ralpha
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // (under_pix * over_ralpha + over_pix * over_alpha)
xmhi = _mm_srli_epi16(xmo16, 8); // (under_pix * over_ralpha + over_pix * over_alpha) >> 8
xmhi = _mm_insert_epi16(xmhi, 255, 3); //write alpha
xmhi = _mm_insert_epi16(xmhi, 255, 7); //write alpha
} else {
//
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmunder_alpha = _mm_shufflelo_epi16(xmu16, _MM_SHUFFLE(3,3,3,3));
xmunder_alpha = _mm_shufflehi_epi16(xmunder_alpha, _MM_SHUFFLE(3,3,3,3));
xmunder_ralpha = _mm_sub_epi16(xmff16, xmunder_alpha);
xmalpha = _mm_sub_epi16(xmff16, _mm_srli_epi16( _mm_mullo_epi16(xmover_ralpha, xmunder_ralpha), 8)); //255 - ((over_ralpha * under_ralpha) >> 8)
xmu16 = _mm_mullo_epi16( _mm_srli_epi16( _mm_mullo_epi16(xmunder_alpha, xmover_ralpha), 8), xmu16); //under_pix * ((under_alpha * over_ralpha)>>8)
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // over_pix * over_alpha + (under_pix * ((under_alpha * over_ralpha)>>8)
//calc pixel1
xmf_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmo16, xmzero) );
xmf_alpha_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmalpha, xmzero) );
xmf_lo = _mm_mul_ps(xmf_lo, _mm_rcp_ps(xmf_alpha_lo));
//calc pixel 2
xmf_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmo16, xmzero) );
xmf_alpha_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmalpha, xmzero) );
xmf_hi = _mm_mul_ps(xmf_hi, _mm_rcp_ps(xmf_alpha_hi));
//pack to 16bit data
xmhi = _mm_packs_epi32(_mm_cvtps_epi32(xmf_lo), _mm_cvtps_epi32(xmf_hi));
int a3 = _mm_extract_epi16(xmalpha, 3);
xmhi = _mm_insert_epi16(xmhi, a3, 3); //write alpha
int a4 = _mm_extract_epi16(xmalpha, 7);
xmhi = _mm_insert_epi16(xmhi, a4, 7); //write alpha
}
//16bit幅、2ピクセルのデータが入った2つのレジスタを合わせて8bitへ
//4ピクセルのデータをメモリへストア
_mm_storeu_ps((float*)dst_pix,
_mm_castsi128_ps( _mm_packus_epi16(xmlo, xmhi)) );
//
over_pix += 4;
under_pix += 4;
dst_pix += 4;
}
over_line = GetNextLineAddress32(img1, over_line);
under_line = GetNextLineAddress32(img2, under_line);
dst_line = GetNextLineAddress32(img4, dst_line);
}
}
std::cout << timer.elapsed() << " msec" << std::endl;
cvSaveImage("test4.bmp", img4);
//SSE: _mm_load_ps + _mm_store_ps
timer.restart();
{
over_line = GetPixelAddress32(img1, 0, 0);
under_line = GetPixelAddress32(img2, 0, 0);
dst_line = GetPixelAddress32(img5, 0, 0);
__m128i xmzero = _mm_setzero_si128();
__m128i xmff16 = _mm_set1_epi16(0xff);
for(y=0; y<img1->height; y++){
over_pix = over_line;
under_pix = under_line;
dst_pix = dst_line;
for(x=0; x<img1->width; x+=4){
__m128i under = _mm_castps_si128( _mm_load_ps((float*)under_pix) );
__m128i over = _mm_castps_si128( _mm_load_ps((float*)over_pix) );
__m128i xmo16, xmu16;
__m128i xmlo, xmhi;
__m128i xmover_alpha, xmover_ralpha, xmunder_alpha, xmunder_ralpha, xmalpha;
__m128 xmf_lo, xmf_hi, xmf_alpha_lo, xmf_alpha_hi;
//load lo double ward
xmo16 = _mm_unpacklo_epi8(over, xmzero);
xmu16 = _mm_unpacklo_epi8(under, xmzero);
int over_a1 = _mm_extract_epi16(xmo16, 3);
int over_a2 = _mm_extract_epi16(xmo16, 7);
int under_a1 = _mm_extract_epi16(xmu16, 3);
int under_a2 = _mm_extract_epi16(xmu16, 7);
if (((over_a1 == 255 && over_a2 == 255)) || (under_a1 == 0 && under_a2 == 0)) {
xmlo = xmo16;
} else if (under_a1 == 0 && under_a2 == 0) {
xmlo = xmu16;
} else if (under_a1 == 255 && under_a2 == 255) {
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmu16 = _mm_mullo_epi16(xmu16, xmover_ralpha); // under_pix * over_ralpha
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // (under_pix * over_ralpha + over_pix * over_alpha)
xmlo = _mm_srli_epi16(xmo16, 8); // (under_pix * over_ralpha + over_pix * over_alpha) >> 8
xmlo = _mm_insert_epi16(xmlo, 255, 3); //write alpha
xmlo = _mm_insert_epi16(xmlo, 255, 7); //write alpha
} else {
//
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmunder_alpha = _mm_shufflelo_epi16(xmu16, _MM_SHUFFLE(3,3,3,3));
xmunder_alpha = _mm_shufflehi_epi16(xmunder_alpha, _MM_SHUFFLE(3,3,3,3));
xmunder_ralpha = _mm_sub_epi16(xmff16, xmunder_alpha);
xmalpha = _mm_sub_epi16(xmff16, _mm_srli_epi16( _mm_mullo_epi16(xmover_ralpha, xmunder_ralpha), 8)); //255 - ((over_ralpha * under_ralpha) >> 8)
xmu16 = _mm_mullo_epi16( _mm_srli_epi16( _mm_mullo_epi16(xmunder_alpha, xmover_ralpha), 8), xmu16); //under_pix * ((under_alpha * over_ralpha)>>8)
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // over_pix * over_alpha + (under_pix * under_alpha * over_ralpha)>>8
//calc pixel 1
xmf_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmo16, xmzero) );
xmf_alpha_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmalpha, xmzero) );
xmf_lo = _mm_mul_ps(xmf_lo, _mm_rcp_ps(xmf_alpha_lo));
//calc pixel 2
xmf_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmo16, xmzero) );
xmf_alpha_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmalpha, xmzero) );
xmf_hi = _mm_mul_ps(xmf_hi, _mm_rcp_ps(xmf_alpha_hi));
//pack to 16bit data
xmlo = _mm_packs_epi32(_mm_cvtps_epi32(xmf_lo), _mm_cvtps_epi32(xmf_hi));
int a1 = _mm_extract_epi16(xmalpha, 3);
xmlo = _mm_insert_epi16(xmlo, a1, 3); //write alpha
int a2 = _mm_extract_epi16(xmalpha, 7);
xmlo = _mm_insert_epi16(xmlo, a2, 7); //write alpha
}
//load hi double ward
xmo16 = _mm_unpackhi_epi8(over, xmzero);
xmu16 = _mm_unpackhi_epi8(under, xmzero);
over_a1 = _mm_extract_epi16(xmo16, 3);
over_a2 = _mm_extract_epi16(xmo16, 7);
under_a1 = _mm_extract_epi16(xmu16, 3);
under_a2 = _mm_extract_epi16(xmu16, 7);
if (((over_a1 == 255 && over_a2 == 255)) || (under_a1 == 0 && under_a2 == 0)) {
xmhi = xmo16;
} else if (under_a1 == 0 && under_a2 == 0) {
xmhi = xmu16;
} else if (under_a1 == 255 && under_a2 == 255) {
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmu16 = _mm_mullo_epi16(xmu16, xmover_ralpha); // under_pix * over_ralpha
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // (under_pix * over_ralpha + over_pix * over_alpha)
xmhi = _mm_srli_epi16(xmo16, 8); // (under_pix * over_ralpha + over_pix * over_alpha) >> 8
xmhi = _mm_insert_epi16(xmhi, 255, 3); //write alpha
xmhi = _mm_insert_epi16(xmhi, 255, 7); //write alpha
} else {
//
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmunder_alpha = _mm_shufflelo_epi16(xmu16, _MM_SHUFFLE(3,3,3,3));
xmunder_alpha = _mm_shufflehi_epi16(xmunder_alpha, _MM_SHUFFLE(3,3,3,3));
xmunder_ralpha = _mm_sub_epi16(xmff16, xmunder_alpha);
xmalpha = _mm_sub_epi16(xmff16, _mm_srli_epi16( _mm_mullo_epi16(xmover_ralpha, xmunder_ralpha), 8)); //255 - ((over_ralpha * under_ralpha) >> 8)
xmu16 = _mm_mullo_epi16( _mm_srli_epi16( _mm_mullo_epi16(xmunder_alpha, xmover_ralpha), 8), xmu16); //under_pix * ((under_alpha * over_ralpha)>>8)
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // over_pix * over_alpha + (under_pix * ((under_alpha * over_ralpha)>>8)
//calc pixel1
xmf_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmo16, xmzero) );
xmf_alpha_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmalpha, xmzero) );
xmf_lo = _mm_mul_ps(xmf_lo, _mm_rcp_ps(xmf_alpha_lo));
//calc pixel 2
xmf_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmo16, xmzero) );
xmf_alpha_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmalpha, xmzero) );
xmf_hi = _mm_mul_ps(xmf_hi, _mm_rcp_ps(xmf_alpha_hi));
//pack to 16bit data
xmhi = _mm_packs_epi32(_mm_cvtps_epi32(xmf_lo), _mm_cvtps_epi32(xmf_hi));
int a3 = _mm_extract_epi16(xmalpha, 3);
xmhi = _mm_insert_epi16(xmhi, a3, 3); //write alpha
int a4 = _mm_extract_epi16(xmalpha, 7);
xmhi = _mm_insert_epi16(xmhi, a4, 7); //write alpha
}
//store to memory
_mm_store_ps((float*)dst_pix,
_mm_castsi128_ps( _mm_packus_epi16(xmlo, xmhi)) );
//
over_pix += 4;
under_pix += 4;
dst_pix += 4;
}
over_line = GetNextLineAddress32(img1, over_line);
under_line = GetNextLineAddress32(img2, under_line);
dst_line = GetNextLineAddress32(img5 , dst_line);
}
}
std::cout << timer.elapsed() << " msec" << std::endl;
cvSaveImage("test5.bmp", img5);
//SSE: _mm_load_ps + _mm_stream_si128
timer.restart();
{
over_line = GetPixelAddress32(img1, 0, 0);
under_line = GetPixelAddress32(img2, 0, 0);
dst_line = GetPixelAddress32(img6, 0, 0);
__m128i xmzero = _mm_setzero_si128();
__m128i xmff16 = _mm_set1_epi16(0xff);
for(y=0; y<img1->height; y++){
over_pix = over_line;
under_pix = under_line;
dst_pix = dst_line;
for(x=0; x<img1->width; x+=4){
__m128i under = _mm_castps_si128( _mm_load_ps((float*)under_pix) );
__m128i over = _mm_castps_si128( _mm_load_ps((float*)over_pix) );
__m128i xmo16, xmu16;
__m128i xmlo, xmhi;
__m128i xmover_alpha, xmover_ralpha, xmunder_alpha, xmunder_ralpha, xmalpha;
__m128 xmf_lo, xmf_hi, xmf_alpha_lo, xmf_alpha_hi;
//load lo double ward
xmo16 = _mm_unpacklo_epi8(over, xmzero);
xmu16 = _mm_unpacklo_epi8(under, xmzero);
int over_a1 = _mm_extract_epi16(xmo16, 3);
int over_a2 = _mm_extract_epi16(xmo16, 7);
int under_a1 = _mm_extract_epi16(xmu16, 3);
int under_a2 = _mm_extract_epi16(xmu16, 7);
if (((over_a1 == 255 && over_a2 == 255)) || (under_a1 == 0 && under_a2 == 0)) {
xmlo = xmo16;
} else if (under_a1 == 0 && under_a2 == 0) {
xmlo = xmu16;
} else if (under_a1 == 255 && under_a2 == 255) {
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmu16 = _mm_mullo_epi16(xmu16, xmover_ralpha); // under_pix * over_ralpha
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // (under_pix * over_ralpha + over_pix * over_alpha)
xmlo = _mm_srli_epi16(xmo16, 8); // (under_pix * over_ralpha + over_pix * over_alpha) >> 8
xmlo = _mm_insert_epi16(xmlo, 255, 3); //write alpha
xmlo = _mm_insert_epi16(xmlo, 255, 7); //write alpha
} else {
//
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmunder_alpha = _mm_shufflelo_epi16(xmu16, _MM_SHUFFLE(3,3,3,3));
xmunder_alpha = _mm_shufflehi_epi16(xmunder_alpha, _MM_SHUFFLE(3,3,3,3));
xmunder_ralpha = _mm_sub_epi16(xmff16, xmunder_alpha);
xmalpha = _mm_sub_epi16(xmff16, _mm_srli_epi16( _mm_mullo_epi16(xmover_ralpha, xmunder_ralpha), 8)); //255 - ((over_ralpha * under_ralpha) >> 8)
xmu16 = _mm_mullo_epi16( _mm_srli_epi16( _mm_mullo_epi16(xmunder_alpha, xmover_ralpha), 8), xmu16); //under_pix * ((under_alpha * over_ralpha)>>8)
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // over_pix * over_alpha + (under_pix * under_alpha * over_ralpha)>>8
//calc pixel 1
xmf_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmo16, xmzero) );
xmf_alpha_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmalpha, xmzero) );
xmf_lo = _mm_mul_ps(xmf_lo, _mm_rcp_ps(xmf_alpha_lo));
//calc pixel 2
xmf_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmo16, xmzero) );
xmf_alpha_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmalpha, xmzero) );
xmf_hi = _mm_mul_ps(xmf_hi, _mm_rcp_ps(xmf_alpha_hi));
//pack to 16bit data
xmlo = _mm_packs_epi32(_mm_cvtps_epi32(xmf_lo), _mm_cvtps_epi32(xmf_hi));
//write alpha
int a1 = _mm_extract_epi16(xmalpha, 3);
xmlo = _mm_insert_epi16(xmlo, a1, 3); //write alpha
int a2 = _mm_extract_epi16(xmalpha, 7);
xmlo = _mm_insert_epi16(xmlo, a2, 7); //write alpha
}
//load hi double ward
xmo16 = _mm_unpackhi_epi8(over, xmzero);
xmu16 = _mm_unpackhi_epi8(under, xmzero);
over_a1 = _mm_extract_epi16(xmo16, 3);
over_a2 = _mm_extract_epi16(xmo16, 7);
under_a1 = _mm_extract_epi16(xmu16, 3);
under_a2 = _mm_extract_epi16(xmu16, 7);
if (((over_a1 == 255 && over_a2 == 255)) || (under_a1 == 0 && under_a2 == 0)) {
xmhi = xmo16;
} else if (under_a1 == 0 && under_a2 == 0) {
xmhi = xmu16;
} else if (under_a1 == 255 && under_a2 == 255) {
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmu16 = _mm_mullo_epi16(xmu16, xmover_ralpha); // under_pix * over_ralpha
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // (under_pix * over_ralpha + over_pix * over_alpha)
xmhi = _mm_srli_epi16(xmo16, 8); // (under_pix * over_ralpha + over_pix * over_alpha) >> 8
xmhi = _mm_insert_epi16(xmhi, 255, 3); //write alpha
xmhi = _mm_insert_epi16(xmhi, 255, 7); //write alpha
} else {
//
xmover_alpha = _mm_shufflelo_epi16(xmo16, _MM_SHUFFLE(3,3,3,3));
xmover_alpha = _mm_shufflehi_epi16(xmover_alpha, _MM_SHUFFLE(3,3,3,3));
xmover_ralpha = _mm_sub_epi16(xmff16, xmover_alpha);
xmunder_alpha = _mm_shufflelo_epi16(xmu16, _MM_SHUFFLE(3,3,3,3));
xmunder_alpha = _mm_shufflehi_epi16(xmunder_alpha, _MM_SHUFFLE(3,3,3,3));
xmunder_ralpha = _mm_sub_epi16(xmff16, xmunder_alpha);
xmalpha = _mm_sub_epi16(xmff16, _mm_srli_epi16( _mm_mullo_epi16(xmover_ralpha, xmunder_ralpha), 8)); //255 - ((over_ralpha * under_ralpha) >> 8)
xmu16 = _mm_mullo_epi16( _mm_srli_epi16( _mm_mullo_epi16(xmunder_alpha, xmover_ralpha), 8), xmu16); //under_pix * ((under_alpha * over_ralpha)>>8)
xmo16 = _mm_mullo_epi16(xmo16, xmover_alpha); // over_pix * over_alpha
xmo16 = _mm_add_epi16(xmu16, xmo16); // over_pix * over_alpha + (under_pix * ((under_alpha * over_ralpha)>>8)
//calc pixel1
xmf_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmo16, xmzero) );
xmf_alpha_lo = _mm_cvtepi32_ps( _mm_unpacklo_epi16(xmalpha, xmzero) );
xmf_lo = _mm_mul_ps(xmf_lo, _mm_rcp_ps(xmf_alpha_lo));
//calc pixel 2
xmf_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmo16, xmzero) );
xmf_alpha_hi = _mm_cvtepi32_ps( _mm_unpackhi_epi16(xmalpha, xmzero) );
xmf_hi = _mm_mul_ps(xmf_hi, _mm_rcp_ps(xmf_alpha_hi));
//pack to 16bit data
xmhi = _mm_packs_epi32(_mm_cvtps_epi32(xmf_lo), _mm_cvtps_epi32(xmf_hi));
int a3 = _mm_extract_epi16(xmalpha, 3);
xmhi = _mm_insert_epi16(xmhi, a3, 3); //write alpha
int a4 = _mm_extract_epi16(xmalpha, 7);
xmhi = _mm_insert_epi16(xmhi, a4, 7); //write alpha
}
//store to memory
_mm_stream_si128((__m128i*)dst_pix, _mm_packus_epi16(xmlo, xmhi));
//
over_pix += 4;
under_pix += 4;
dst_pix += 4;
}
over_line = GetNextLineAddress32(img1, over_line);
under_line = GetNextLineAddress32(img2, under_line);
dst_line = GetNextLineAddress32(img6, dst_line);
}
}
std::cout << timer.elapsed() << " msec" << std::endl;
cvSaveImage("test6.bmp", img6);
cvReleaseImage(&b_ch);
cvReleaseImage(&g_ch);
cvReleaseImage(&r_ch);
cvReleaseImage(&a_ch1);
cvReleaseImage(&a_ch2);
cvReleaseImage(&img1);
cvReleaseImage(&img2);
cvReleaseImage(&img3);
cvReleaseImage(&img4);
cvReleaseImage(&img5);
cvReleaseImage(&img6);
}