エッシャーっぽい絵を生成する "escherize" を OpenCV で作ってみた

インスパイア元 --> 映像奮闘記: エッシャーっぽい絵を生成する「エッシャーくん」を作ってみた。
画像を変換してエッシャーっぽく無限渦巻きにする，というプログラム．面白そうだったので，OpenCV を使って作ってみました．というか，上のサイトの python で書かれていたコードをほとんどそのまま C++ に落としただけです．コードは記事の最後に載っけておきます．
変換のアルゴリズムの詳細はこちら --> "A logarithmic image transformation"．要は，入力画像からドーナツ型の領域を切り取って，一ヶ所ハサミを入れて，渦巻き一周分になるようにぐねっと歪めて，拡縮を繰り替えしながら貼り付けていく，と．これを複素平面での単純な変換でやっているのがすごい．
変換結果は下のようになります．補間をしていないので，よく見るとちょっとジャギーな仕上がりです．

ちなみに，のコードはこちら．

#include <math.h>
#include <libgen.h>
#include <stdio.h>

#include <cxcore.h>
#include <cv.h>
#include <highgui.h>

/**
 * エッシャー風の画像に変換する
 * src と dst は同じサイズでなくても，同じタイプでなくても大丈夫
 * @param src 入力画像
 * @param dst 出力画像
 * @param r1  変換に使用する領域の内径
 * @param r2  変換に使用する領域の外径
 */
void escherize(CvMat *src, CvMat *dst, double r1, double r2) {
  CV_FUNCNAME("escherize");
  __BEGIN__

  double tw, th, alpha, f, ca, sa;

  // 引数を検査する
  if (r1 >= r2)
    CV_ERROR(CV_StsBadArg, "inner radius is too long or outer raidus is too short");
  if (r1 <= 0.0)
    CV_ERROR(CV_StsOutOfRange, "inner radius is too short or negative");
  if (r2 >= src->rows/2.0 || r2 >= src->cols/2.0)
    CV_ERROR(CV_StsOutOfRange, "outer radius is too long");

  // 計算に必要な定数を求める
  tw    = log(r2/r1);  // タイルの幅
  th    = M_PI * 2.0;  // タイルの高さ
  alpha = atan(tw/th); // タイルに変換する際の回転角度
  f     = cos(alpha);  // タイルに変換する際のスケール
  ca    = cos(alpha);
  sa    = sin(alpha);

  // 出力画像を生成する
  for (int j = 0; j < dst->height; ++j) {
    double dy = j - dst->height / 2.0;
    for (int i = 0; i < dst->width; ++i) {
      double dx = i - dst->width / 2.0;
      // 極座標に変換する
      double r = sqrt(dx*dx + dy*dy);
      if (r < 0.5) {
        cvSet2D(dst, j, i, cvScalarAll(0.0));
        continue;
      }
      double theta = atan2(dy, dx);
      // タイル空間に変換する
      double lr = log(r);
      double re = (lr *  ca + theta * sa) / f;
      double im = (lr * -sa + theta * ca) / f;
      // 中心のタイル範囲に収める
      while (re > tw ) re -= tw;
      while (re < 0.0) re += tw;
      while (im > th ) im -= th;
      while (im < 0.0) im += th;
      // 入力画像の座標に変換する
      double sx = r1 * exp(re) * cos(im) + src->width  / 2.0;
      double sy = r1 * exp(re) * sin(im) + src->height / 2.0;
      // 出力画像の画素値を求める
      // TODO: 補間したい
      cvSet2D(dst, j, i, cvGet2D(src, (int)sy, (int)sx));
    }
  }

  __END__
}

/**
 * 使い方を表示して，プログラムを終了する
 * @param status  終了時に親プロセスに返す値
 * @param command コマンド名
 */
void usage(int status, const char* command) {
  printf("Usage: %s src dst r1 r2\n", command);
  printf("Arguments:\n"
         "  src      source image file\n"
         "  dst      dest image file\n"
         "  r1, r2   transform parameters (r1 < r2 < image_size/2)\n");
  exit(status);
}

/**
 * メイン関数
 * @param argc コマンドライン引数の数
 * @param argv コマンドライン引数の配列
 * @retval 終了ステータス
 */
int main(int argc, char** argv) {
  CvMat *src = 0, *dst = 0;
  __BEGIN__

  double r1, r2;
  const char *src_file, *dst_file;

  // コマンドライン引数をパースする
  if (argc != 5)
    usage(EXIT_FAILURE, basename(argv[0]));
  src_file = argv[1];
  dst_file = argv[2];
  r1 = atof(argv[3]);
  r2 = atof(argv[4]);

  // 画像を読み込む
  src = cvLoadImageM(src_file, CV_LOAD_IMAGE_UNCHANGED);

  // 出力領域を確保する
  dst = cvCreateMat(src->rows, src->cols, CV_MAT_TYPE(src->type));

  // エッシャー風に変換する
  escherize(src, dst, r1, r2);

  // 出力画像を保存する
  cvSaveImage(dst_file, dst);

  __END__
  cvReleaseMat(&src);
  cvReleaseMat(&dst);
}

コンパイルと実行は下のような感じで．

$ g++ `pkg-config opencv --cflags` -o escherize escherize.cc `pkg-config opencv --libs` -lm
$ ./escherize lena.png lena_escherize.png 50 150