Bilinear(バイリニア法) ～3.アセンブラ＆MMX編～

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

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

MMXって何だっけ？

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

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

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

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

MMX - Wikipedia

ja.wikipedia.org

で、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という補間アルゴリズムのサンプルを作ってみます！