[Delphi(Win32)] Bilinear(バイリニア法) ~3.アセンブラ&MMX編~

  • このエントリーをはてなブックマークに追加

久しぶりのDelphiネタ。そしてBilinear補間処理の続き。
何年ぶりだろ?ようやくDelphi復帰宣言です。

今回はかなり時代遅れな感じが否めないMMXによる最適化です。
本当はSSEとかマルチスレッドとか使うべきなんだろうけど・・・

MMXって何だっけ?

もう世間では当たり前になっているので手短に済まします。

複数のデータを一挙に処理する事ができるCPUの拡張命令セットです。

整数演算を同時処理する為のもので、古くからあるので機能的には結構物足りないところも多いです。

まあ何だか分からなければ下記のWiki記事読んだりYahoo!でググってくださいw
もうそこらのPentium互換CPUは必ず搭載してて当たり前になっているので、特別説明することも無いでしょう。

MMX – Wikipedia

で、MMXを利用するにはIntelのお高いコンパイラとかを駆使すれば高級言語からでもいい感じになるのでしょうが、金かけたくないですし、手動チューニングの面白さを尊重したいと思います。
そこでアセンブラでのコード記述が必要となってきます。流石にプログラム全体をアセンブラで書くのは無謀過ぎるので、インラインアセンブラを利用します。

コードサンプル

インラインアセンブラはDelphiでも当然サポートされているので使ってみました。VisualStudioでC++などで同様にインラインアセンブラを使うときもほとんど同じコードを流用できるので、他の言語使いの方でも多少は参考になるかと思います。
ただ、x86アセンブラ初挑戦な上に、4年も前に組んだシロモノ・・・大分至らぬところはあると思います。
今回も前回の固定小数点数化を踏襲して整数演算を主体に処理しています。

function ResizeBilinearMMX(SrcBmp, DstBmp: TBitmap; Width, Height: LongWord): Boolean;
const
  FIX_SHIFT      = 7;
  FIX_SHIFT32    = 15;
  FIX_SHIFT12    = 11;
  FIX_DEC_MASK32 = $00007FFF;
var
  FIX_ONE         : LongWord;// 8bitの1
  FIX_DEC_MASK12  : LongWord;// 12bitの小数部
  FIX_ONE12       : LongWord;// 12bitの1

  wfactor, hfactor: LongWord;
  pSrcBase        : Pointer;
  pDstBase        : Pointer;
  srcOffset       : LongInt;
  dstOffset       : LongInt;
  wMax, hMax      : LongWord;
  xPos, xCoef     : array of LongWord;
  ix, src_ix      : LongWord;
begin
  Result := False;

  if (Width < 1) or (Height < 1) then Exit;

  if (DstBmp.PixelFormat <> pf32bit)
  or (SrcBmp.PixelFormat <> pf32bit) then Exit;

  FIX_ONE         := $00000000;// MMX用定数
  FIX_DEC_MASK12  := $000007FF;//
  FIX_ONE12       := $00000800;//

  DstBmp.Width  := Width;                                   // 出力BMPサイズ設定
  DstBmp.Height := Height;                                  //

  SetLength(xPos,  Width + 1);                              //
  SetLength(xCoef, Width);                                  //

  hfactor := LongWord(SrcBmp.Height shl FIX_SHIFT12) div Height; // 入出力サイズ比率を求める(Q16)
  wfactor := LongWord(SrcBmp.Width shl FIX_SHIFT32) div Width;   //

  wMax := SrcBmp.Width - 1;                                 // 横最大インデックス
  hMax := SrcBmp.Height - 1;                                // 縦最大インデックス

  pSrcBase := SrcBmp.ScanLine[0];
  srcOffset := Integer(SrcBmp.ScanLine[1]) - Integer(pSrcBase);

  pDstBase := DstBmp.ScanLine[0];
  dstOffset := Integer(DstBmp.ScanLine[1]) - Integer(pDstBase);

  // 横のピクセル位置
  src_ix := 0;
  for ix := 0 to Width - 1 do
  begin
    xPos[ix]  := src_ix shr FIX_SHIFT32;                    // x座標
    xCoef[ix] := (src_ix and FIX_DEC_MASK32) shr 8;         // x座標における係数
    src_ix := src_ix + wfactor;                             //
  end;

  // 補間処理
  asm
    push eax; push ebx;         // レジスタ退避
    push ecx; push edx;         //
    push edi; push esi;         //

    // iy loop -----------------------------------------------------------------
    xor ecx, ecx;               // iy
    xor ebx, ebx;               // src_iy

    // src_iyはhfactorを足していく
    @Loop_iy:

      push ecx;                 // dst_iy退避
      push ebx;                 // src_iy退避

      // yCoef(8bit)
      movd mm7, ebx;            //
      movd mm0, FIX_DEC_MASK12; //
      movd mm2, FIX_ONE12;      //
      pand mm7, mm0;            // mm7 = yCoef1
      psubd mm2, mm7;           // mm2 = yCoef2

      punpcklwd mm7, mm7;       // mm7 = 0 0 yc1 yc1
      punpcklwd mm2, mm2;       // mm2 = 0 0 yc2 yc2
      punpckldq mm7, mm2;       // mm7 = yc2 yc2 yc1 yc1

      shr ebx, FIX_SHIFT12;     // yPos0 = src_iy >> 11
      mov edx, pSrcBase;        //
      mov edi, srcOffset;       //
      mov eax, ebx;             // ebx退避
      imul ebx, edi;            // ebx = srcOffset*yPos0
      lea esi, [ebx + edx];     // pSrcY0 = pSrcBase + ebx

      cmp eax, hMax;            // yPos1 = Min(hMax, yPos0+1 >> 15)
      jae @below_hMax;          //
      add ebx, edi;             //
      @below_hMax:              //
      lea edi, [ebx + edx];     // pSrcY1 = pSrcBase - 4*yPos1

      // ix loop ---------------------------------------------------------------
      xor ecx, ecx;             //
      mov edx, pDstBase;        //

      @Loop_ix:

        mov ebx, xPos;          //
        mov eax, [ebx + 4*ecx]; //
        movd mm2, [esi + 4*eax];// (xPos0, yPos0)
        movd mm4, [edi + 4*eax];// (xPos0, yPos1)

        cmp eax, wMax;          // xPos1 = Min(wMax, xPos0+1)
        jae @below_wMax;        //
        inc eax;                //
        @below_wMax:            //
        movd mm3, [esi + 4*eax];// (xPos1, yPos0)
        movd mm5, [edi + 4*eax];// (xPos1, yPos1)

        // pixel data unpack
        pxor mm0, mm0;          //
        punpcklbw mm2, mm0;     // mm2 = 0 A1 0 B1 0 G1 0 R1
        punpcklbw mm3, mm0;     // mm3 = 0 A2 0 B2 0 G2 0 R2
        movq mm1, mm2;          //
        punpcklwd mm2, mm3;     // mm2 = 0 G2 0 G1 0 R2 0 R1
        punpckhwd mm1, mm3;     // mm1 = 0 A2 0 A1 0 B2 0 B1

        punpcklbw mm4, mm0;     // mm4 = 0 A3 0 B3 0 G3 0 R3
        punpcklbw mm5, mm0;     // mm5 = 0 A4 0 B4 0 G4 0 R4
        movq mm3, mm4;          //
        punpcklwd mm4, mm5;     // mm4 = 0 G4 0 G3 0 R4 0 R3
        punpckhwd mm3, mm5;     // mm3 = 0 A4 0 A3 0 B4 0 B3

        // xCoef
        mov eax, xCoef;         //
        movd mm6, [eax + 4*ecx];// xCoef1
        movd mm5, FIX_ONE;      //
        psubw mm5, mm6;         // xCoef0
        punpcklwd mm5, mm5;     // mm5 = 0 0 x1 x1
        punpcklwd mm6, mm6;     // mm6 = 0 0 x0 x0
        punpcklwd mm5, mm6;     // mm5 = x0 x1 x0 x1

        // multi add xCoef (7bit)
        pmaddwd mm1, mm5;       // A2x0 + A1x1 B2x0 + B1x1
        pmaddwd mm3, mm5;       // A4x0 + A3x1 B4x0 + B3x1
        pmaddwd mm2, mm5;       // G2x0 + G1x1 R2x0 + R1x1
        pmaddwd mm4, mm5;       // G4x0 + G3x1 R4x0 + R3x1

        psrld mm1, FIX_SHIFT + FIX_SHIFT12 - 16;   // pmulhwでレンジ0になるように調整
        psrld mm3, FIX_SHIFT + FIX_SHIFT12 - 16;   //
        psrld mm2, FIX_SHIFT + FIX_SHIFT12 - 16;   //
        psrld mm4, FIX_SHIFT + FIX_SHIFT12 - 16;   //

        // multi yCoef(かけたあと上位のみ取り出しでQ0)
        movq mm5, mm7;          //
        movq mm6, mm7;          //
        punpckhwd mm5, mm0;     // yCoef0(12bit)
        punpcklwd mm6, mm0;     // yCoef1(12bit)

        pmulhw mm1, mm5;        // Ay0*yc0 By0*yc0
        pmulhw mm3, mm6;        // Ay1*yc1 By1*yc1
        pmulhw mm2, mm5;        // Gy0*yc0 Ry0*yc0
        pmulhw mm4, mm6;        // Gy1*yc1 Ry1*yc1

        // add (byte単位ラップアラウンド付き加算)
        paddw mm1, mm3;         // A B
        paddw mm2, mm4;         // G R

        // to 4bytes data
        packuswb mm2, mm1;      // A B G R
        packuswb mm2, mm2;      //

        // store to memory
        movd [edx + 4*ecx], mm2;// [pDstBase + 4*dst_iy]

        // end of ix loop
        inc ecx;                // ix++
        cmp ecx, Width;         // if (ix != Width) goto LOOP_IX

      jne @Loop_ix;             //
      // -----------------------------------------------------------------------
      pop ebx;                  //
      pop ecx;                  //

      add ebx, hfactor;         // src_iy += hfactor

      mov eax, dstOffset;       // dstBase += dstOffset
      add pDstBase, eax;        //

      inc ecx;                  // dst_iy++
      cmp ecx, Height;          //

    jne @Loop_iy;               //

    pop esi; pop edi;           // レジスタ復活
    pop edx; pop ecx;           //
    pop ebx; pop eax;           //

    emms;                       // MMX後始末
  end; // end of for of iy
  // ---------------------------------------------------------------------------

  Result := True;
end;

多分これチラっと見ただけじゃあなんだかさっぱりだと思います。
アセンブラはデータの流れを追うのは結構大変です。デバッガでステップ実行しながらCPUのレジスタを見ると何をやっているのかわかると思われます。
しかしながらパック、アンパックという作業は複数のビットが一気に交換されるので混乱してしまうかもしれません。

結果的には普通にDelphiのコードを記述したときより圧倒的に早くなりました(10倍位)。PC環境が大幅に変わってしまって、以前の結果と比較が出来なくなってしまったので、具体的な数値は出せません。

なお、MMXはFPU(浮動小数点演算ユニット)とレジスタを共有するので、後片付け処理(emms命令)をしなければなりませんし(初期化は自動的に行われる)、MMX使用中は浮動小数点演算ができません。
なのでMMXの利用は今日においては非推奨なんですよね…

で、そもそもアセンブラ分からない人は参考書読め!とかしか言い様が無いです。CPUアーキテクチャを理解しなければならないので、一言で片付くレベルではありません。

MMXを初めとしたSIMDで最適化したロジックは実に奥深いです。やっぱプロのアセンブラ使いはこんな付け焼刃レベルじゃないんでしょうね。

次回はBicubicという補間アルゴリズムのサンプルを作ってみます!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。




コメント

  1. たかおファン(surface0) より:

    ようやくですが、先日Delphi 10 Starterが無償入手できたので試してみよう。。。
    と思ったんですが、そもそもDelphiインストール直後の素の状態で
    http://edn.embarcadero.com/jp/article/39704
    ここのCode#002をまんま試した所、32bit透過bmpファイルを開いても何も表示されてませんでした。
    (しかもそのままSaveToFileすると透過情報が無くなる)
    もはやTBitmapの扱い自体が謎なので、上記コードの検証には至りませんでした。。。

  2. すー より:

    正直めんどくさいですよね。
    もし、お暇とやる気ができたらお願いします。
    気長にお待ちしてますので。

  3. たかおファン(surface) より:

    >すーさん
    他の方法で大丈夫であるとなれば、やはりアセンブラコードがなにか悪さしてるんですね。
    ここまで来るとステップ実行してレジスタやメモリの変化を見てくという細かい作業が必要そうです。

    はたまたコンパイルオプションを何かいじらないといけないとかあったりして…

    Delphi XE1は持ってるのでやろうとすれば検証可能ですが、ちょっと今入れる余裕が無くて…

  4. すー より:

    昔の記事に今更何度もすいません。
    以前のコード(基本/固定小数)を試しましたが両方とも問題ありませんでした。
    真っ黒になってしまった画像をPhotoshopで確認すると、RGBAとも全て黒でした。
    アルファチャンネルとはいっても他のRGBチャンネル同様にアルファ同士で補間処理すればいいはずなので、32bit画像だからといって特別違わないはずなんですが…
    ちなみにDelphiが32bit Bitmapに対応したのは2009からです。

  5. たかおファン(surface) より:

    >すーさん
    なるほど。。。
    実は当時TBitmapが透過をサポートしておらず、PNGも標準で使えませんでした。
    その為アルファチャンネル付きの画像で試してない気がします。

    おそらく透過を扱う場合は乗算済みアルファの考慮が必要になると思うので、そのあたりは再考が必要になると思います。

    ちなみにこれより以前のバイリニア法についての記事で掲載したコードは試されましたでしょうか?
    もしそちらが大丈夫であれば単に計算途中でオーバーフローしてるだけかもしれません…

  6. すー より:

    アルファチャンネルが原因だろうと思って、いくつかの画像で試しているのですが、同じ症状です。アルファチャンネルを黒一色または白一色にすると、結果は黒または透明になります。

  7. たかおファン(surface) より:

    >すーさん
    こんにちは。
    現在PCにDelphiがインストールされていないので、詳しく検証できないのですが、
    考えられる事として、このコードは32bitビットマップ画像専用の処理になっているので、ビット数が適合していないと正しく処理されませんが、その辺は大丈夫でしょうか?

  8. すー より:

    こんにちわ。いつも参考にさせてもらってます。
    今日サンプルコードを動かしてみようとしたんですが、
    どうやってもイメージが真っ黒(または透明)にしかなりません。
    お手間をとらせて申し訳ありませんが、
    なにかヒントを頂ければありがたいです。

  9. ねこ より:

    なるほどです、確かに最近は SSE やマルチスレッドなんかが使われていますね。
    個人的には Delphi より C# の方が読み易いので全然問題ないです。
    固定小数と言ってたのは Windows 以外の環境でも使えるので利便性が高かった為です。
    ともあれ楽しみにしております。

  10. たかおファン(surface) より:

    正直申しますとMMXはもう完全に時代遅れになってしまい、
    最近はMMXが進化し浮動小数点数も扱えるようになったSSE2、SSE3などが主流になっています。
    その為、今や固定小数演算を使う必要性が薄れて来ました。
    さらにマルチコアCPUの普及により、SIMDを使わずともマルチスレッドで十分な高速化が可能になったので、そちらに重点を置いた記事を執筆したいと思っています。
    また最近の志向からサンプルコードの言語はC#の利用を考えていますのでご了承下さいm(_ _)m

  11. ねこ より:

    お返事ありがとうございます。
    なるほど、そうだったんですね。
    アセンブラでなくて大丈夫ですー。
    Bilinear の時みたいな通常版、固定小数版
    みたいな感じであると嬉しいです。
    宜しくお願いします。

  12. たかおファン(surface) より:

    こんにちは、ねこさん!
    コメントありがとうございます。
    そして更新遅くてごめんなさい!
    アセンブラコードをご希望ですか?
    正直申しますと、現在職業が変わった影響でなかなかアセンブラをじっくり書く機会が無くなってしまいました…
    言語問わずのアルゴリズム紹介程度でしたらすぐかけます!

  13. ねこ より:

    Bicubicはいつ頃になりそうでしょうか…?
    大変参考にさせて頂いております。

コメントを残す