この関数は、2つの半透明画像の混合を実装します。 このタスクは、不透明な背景上に画像をアルファチャネルでオーバーレイする一般的なケースであり、おそらく多くの人によく知られています。 ウィキペディアで半透明の画像の混合について読むことができます。 私はこれについて長い間記事を書きました 。
それ以来、私はすでに、その記事の最後にリストされているものよりも、自分のニーズに最適な実装を書くことができました。 そして、この実装は、Cで書かれた新しいバージョンのPillowのalpha_composite（）よりも高速であることが判明しました。 もちろん、私はこれに満足していましたが、Pillowから実装を改善しようとすることにしました。
最初にテストベンチが必要です。
bench.py
from PIL import Image from timeit import repeat from image import paste_composite im1 = Image.open('in1.png') im1.load() im2 = Image.open('in2.png') im2.load() print repeat(lambda: paste_composite(im1.copy(), im2), number=100) print repeat(lambda: Image.alpha_composite(im1, im2), number=100) out1 = im1.copy() paste_composite(out1, im2) out1.save('out1.png') out2 = Image.alpha_composite(im1, im2) out2.save('out2.png')
スクリプトは、現在のフォルダーにあるin1.pngおよびin2.pngファイルを取得し、最初にpaste_composite（）関数、次にalpha_composite（）で100回混合します。 前の記事（ one 、 two ）からファイルを取りました。 paste_composite（）の平均実行時間は235ミリ秒でした。 元のalpha_composite（）の動作時間は400ミリ秒でした。
なぜalpha_composite（）が非常に遅いのか、長い間推測する必要はありませんでした。 実装を見ると、関数が浮動小数点数ですべての計算を実行していることが明らかです。 グラフィックスを扱う場合、浮動小数点数を使用するとアルゴリズムを簡単に記述できますが、そのパフォーマンスは常に整数のパフォーマンスよりも低くなります。 私が最初にしたことは、すべてを整数で動作するように変換することでした：
typedef struct { UINT8 r; UINT8 g; UINT8 b; UINT8 a; } rgba8; Imaging ImagingAlphaComposite(Imaging imDst, Imaging imSrc) { Imaging imOut = ImagingNew(imDst->mode, imDst->xsize, imDst->ysize); for (int y = 0; y < imDst->ysize; y++) { rgba8* dst = (rgba8*) imDst->image[y]; rgba8* src = (rgba8*) imSrc->image[y]; rgba8* out = (rgba8*) imOut->image[y]; for (int x = 0; x < imDst->xsize; x ++) { if (src->a == 0) { *out = *dst; } else { UINT16 blend = dst->a * (255 - src->a); UINT8 outa = src->a + blend / 255; UINT8 coef1 = src->a * 255 / outa; UINT8 coef2 = blend / outa; out->r = (src->r * coef1 + dst->r * coef2) / 255; out->g = (src->g * coef1 + dst->g * coef2) / 255; out->b = (src->b * coef1 + dst->b * coef2) / 255; out->a = outa; } dst++; src++; out++; } } return imOut; }
そしてすぐに5.5倍の増加を受けました。 テストは75ミリ秒で完了しました。
良い結果が得られましたが、写真がどれだけ高品質であるかを確認する必要がありました。 Photoshopの結果を参考にすることほど良いものはありませんでした。 外観上、結果の2つの画像はPhotoshopで得られた画像と見分けがつかなかったため、違いをより正確に視覚化する方法を見つけ出す必要がありました。
diff.py
#!/usr/bin/env python import sys from PIL import Image, ImageMath def chanel_diff(c1, c2): c = ImageMath.eval('127 + c1 - c2', c1=c1, c2=c2).convert('L') return c.point(lambda c: c if 126 <= c <= 128 else 127 + (c - 127) * 10) im1 = Image.open(sys.argv[1]) im2 = Image.open(sys.argv[2]) diff = map(chanel_diff, im1.split(), im2.split()) Image.merge('RGB', diff[:-1]).save('diff.png', optimie=True) diff[-1].convert('RGB').save('diff.alpha.png', optimie=True)
入力時に、スクリプトは写真付きの2つのファイル名を受け入れます。 これらの画像のチャンネルの各ペアには違いがあります。最初に、違いが単純に検索され、灰色の背景（値127）に追加されるため、両方向の偏差を確認できます。 次に、1つの値よりも大きいすべての偏差が10倍に増幅され、視覚的に見えるようになります。 比較結果は2つのファイルに書き込まれます：diff.pngのrgbチャンネル、diff.alpha.pngのグレーの濃淡の形式のアルファチャンネル。
alpha_composite（）とPhotoshopの結果の違いは非常に重要であることがわかりました。
もう一度コードを注意深く見て、丸めを忘れていることに気付きました。整数は、残りを破棄することで除算されます。 丸められた結果を取得する場合、除数の半分を被除数に追加する必要があります。
UINT16 blend = dst->a * (255 - src->a); UINT8 outa = src->a + (blend + 127) / 255; UINT8 coef1 = src->a * 255 / outa; UINT8 coef2 = blend / outa; out->r = (src->r * coef1 + dst->r * coef2 + 127) / 255; out->g = (src->g * coef1 + dst->g * coef2 + 127) / 255; out->b = (src->b * coef1 + dst->b * coef2 + 127) / 255; out->a = outa;
ランタイムは94ミリ秒に増加しました。 しかし、リファレンスとの違いははるかに少ないです。 しかし、それらはまったく消えませんでした。
原則として、ここでやめることができます。この時点で存在する差異は、主に透明度の高いピクセルに集中しているため、カラーチャンネルの値はそれほど重要ではありません。 ただし、念のため、alpha_composite（）の元のバージョンを確認すると、実際には標準から逸脱していないことがわかりました。 そこで止めることは不可能でした。
明らかに、得られた不正確さは乗算中の係数の不十分な精度の結果です。 各変数に追加の有効ビットを格納でき、除算の精度が向上します。 ブレンド変数は除算なしで計算され、何もする必要はありません。 outa変数は、ブレンドを255で除算することによって取得されます。計算時には、blendとsrc-> aの値を4ビット上にシフトします。 その結果、outaには12の有効ビットがあります。 両方の係数で、outaによる除算が発生し、現在は4ビット大きくなっています。 さらに、係数自体はさらに4ビットを取得したいと考えています。 その結果、被除数はすでに8ビットシフトされています。
UINT16 blend = dst->a * (255 - src->a); // 16 bit max UINT16 outa = (src->a << 4) + ((blend << 4) + 127) / 255; // 12 UINT16 coef1 = ((src->a * 255) << 8) / outa; // 12 UINT16 coef2 = (blend << 8) / outa; // 12 out->r = ((src->r * coef1 + dst->r * coef2 + 0x7ff) / 255) >> 4; out->g = ((src->g * coef1 + dst->g * coef2 + 0x7ff) / 255) >> 4; out->b = ((src->b * coef1 + dst->b * coef2 + 0x7ff) / 255) >> 4; out->a = (outa + 0x7) >> 4;
もちろん、同じ4ビットですべての計算結果をシフトしてから画像ピクセルに配置する必要があります。 このオプションの動作は少し遅くなります：108ミリ秒。ただし、標準から複数の値が異なるピクセルはほとんどありません。 品質に関する結果が達成され、パフォーマンスについて再度考えることができます。
明らかに、このコードはプロセッサの観点からは最適ではありません。 分割は最も困難なタスクの1つであり、ここでは各ピクセルに対して最大6分割が実行されます。 それらのほとんどを取り除くことは非常に簡単であることが判明しました。 Paste.cファイルの同じPILで、丸めによる高速除算の方法が説明されています。
a + 127 / 255 ≈ ((a + 128) + ((a + 128) >> 8)) >> 8
その結果、1つの部門が2つのシフトと1つの追加に置き換えられ、ほとんどのプロセッサーで高速になります。 シフトの1つを既存のシフト4にグループ化すると、次のコードが得られました。
UINT16 blend = dst->a * (255 - src->a); UINT16 outa = (src->a << 4) + (((blend << 4) + (blend >> 4) + 0x80) >> 8); UINT16 coef1 = (((src->a << 8) - src->a) << 8) / outa; // 12 UINT16 coef2 = (blend << 8) / outa; // 12 UINT32 tmpr = src->r * coef1 + dst->r * coef2 + 0x800; out->r = ((tmpr >> 8) + tmpr) >> 12; UINT32 tmpg = src->g * coef1 + dst->g * coef2 + 0x800; out->g = ((tmpg >> 8) + tmpg) >> 12; UINT32 tmpb = src->b * coef1 + dst->b * coef2 + 0x800; out->b = ((tmpb >> 8) + tmpb) >> 12; out->a = (outa + 0x7) >> 4;
すでに65ミリ秒で完了しており、以前のバージョンと比較して結果は変わりません。
合計で、作業速度が6倍に改善され、プールリクエストがリポジトリに送信されました。 標準とのより大きな類似性の達成を試みることができます。 たとえば、Photoshopで生成されたアルファチャンネルを完全に正確に繰り返すことができましたが、ピクセルにさらに歪みが生じました。
大きなアップデート
もう一度、私は何が行われたのかを検討し、どのタスクが行われ、どの実装が見られないかについて不明なPhotoshopアルゴリズムを思い付くのは愚かなことだと判断しました。 浮動小数点数のアルゴリズムの動作を参照する必要があります。 いくつかの統計を表示するように、違いを見つけるスクリプトを少し修正しました。
diff.py
#!/usr/bin/env python import sys from PIL import Image, ImageMath def chanel_diff(c1, c2): return ImageMath.eval('127 + c1 - c2', c1=c1, c2=c2).convert('L') def highlight(c): return c.point(lambda c: c if 126 <= c <= 128 else 127 + (c - 127) * 10) im1 = Image.open(sys.argv[1]) im2 = Image.open(sys.argv[2]) diff = map(chanel_diff, im1.split(), im2.split()) if len(sys.argv) >= 4: highlight(Image.merge('RGB', diff[:-1])).save('%s.png' % sys.argv[3]) highlight(diff[-1]).convert('RGB').save('%s.alpha.png' % sys.argv[3]) def stats(ch): return sorted((c, n) for n, c in ch.getcolors()) for ch, stat in zip(['red ', 'grn ', 'blu ', 'alp '], map(stats, diff)): print ch, ' '.join('{}: {:>5}'.format(c, n) for c, n in stat)
また、各ピクセルの組み合わせが一意である16メガピクセルの2つの画像の非常に深刻なテストケースを準備しました。
make_testcase.py
def prepare_test_images(dim): """Plese, be careful with dim > 32. Result image is have dim ** 4 pixels (ie 1Mpx for 32 dim or 4Gpx for 256 dim). """ i1 = bytearray(dim ** 4 * 2) i2 = bytearray(dim ** 4 * 2) res = 255.0 / (dim - 1) rangedim = range(dim) pos = 0 for l1 in rangedim: for l2 in rangedim: for a1 in rangedim: for a2 in rangedim: i1[pos] = int(res * l1) i1[pos + 1] = int(res * a1) i2[pos] = int(res * l2) i2[pos + 1] = int(res * a2) pos += 2 print '%s of %s' % (l1, dim) i1 = Image.frombytes('LA', (dim ** 2, dim ** 2), bytes(i1)) i2 = Image.frombytes('LA', (dim ** 2, dim ** 2), bytes(i2)) return i1.convert('RGBA'), i2.convert('RGBA') im1, im2 = prepare_test_images(63) im1.save('im1.png') im2.save('im2.png')
その後、もう一方の端から始めることにしました。整数アルゴリズムの精度を上げるのではなく、浮動小数点アルゴリズムを整数で動作させようとします。 これは私が始めたところです：
double dsta = dst->a / 255.0; double srca = src->a / 255.0; double blend = dsta * (1.0 - srca); double outa = srca + blend; double coef1 = srca / outa; double coef2 = 1 - coef1; double tmpr = src->r * coef1 + dst->r * coef2; out->r = (UINT8) (tmpr + 0.5); double tmpg = src->g * coef1 + dst->g * coef2; out->g = (UINT8) (tmpg + 0.5); double tmpb = src->b * coef1 + dst->b * coef2; out->b = (UINT8) (tmpb + 0.5); out->a = (UINT8) (outa * 255.0 + 0.5);
そして、これは私が来たものです：
UINT16 blend = dst->a * (255 - src->a); UINT16 outa255 = src->a * 255 + blend; // There we use 7 bits for precision. // We could use more, but we go beyond 32 bits. UINT16 coef1 = src->a * 255 * 255 * 128 / outa255; UINT16 coef2 = 255 * 128 - coef1; #define SHIFTFORDIV255(a)\ ((a >> 8) + a >> 8) UINT32 tmpr = src->r * coef1 + dst->r * coef2 + (0x80 << 7); out->r = SHIFTFORDIV255(tmpr) >> 7; UINT32 tmpg = src->g * coef1 + dst->g * coef2 + (0x80 << 7); out->g = SHIFTFORDIV255(tmpg) >> 7; UINT32 tmpb = src->b * coef1 + dst->b * coef2 + (0x80 << 7); out->b = SHIFTFORDIV255(tmpb) >> 7; out->a = SHIFTFORDIV255(outa255 + 0x80);