rydotの呟''

プログラミングとかCGとかDTMとか適当にいろいろのことを適度にやる気なく綴るはず。

楕円を投影

どっかに浮かんでいる楕円をどっかの平面に投影するとどうなるの？
というのを聞かれたので概略を。

楕円の表現

楕円の乗っている平面の法線向き \bf{n}, \: |\bf{n}|=1
楕円の中心点 \bf{p}
長径の向き \bf{x}, \: |\bf{x}|=1, \: \bf{x}\cdot\bf{n}=0
長径a
短径b
として、楕円E(\bf{n},\bf{p},\bf{x},a,b)を \{\bf{p}+a\bf{x}\cos t+b\bf{y}\sin t \: | \: 0 \leq t < \pi\}とパラメトリック形式で表す。
ここで \bf{y}=\bf{n} \times \bf{x}, \: |\bf{y}|=1, \: \bf{x}\cdot\bf{y}=0
「\cdot」は内積、「\times」は外積を表す。

平面の表現

平面の法線向き \bf{m}, \: |\bf{m}|=1
原点からの距離h
として、平面P(\bf{m},h)を \{\bf{r} \: | \: \bf{r}\cdot\bf{m}-h=0\}と陰関数形式で表す。

楕円を平面に投影

点qを平面Pに投影した点q'は \bf{q}-(\bf{q}\cdot\bf{m})\bf{m}+h\bf{m}となる。
楕円の各点を投影すると (\bf{p}+a\bf{x}\cos t+b\bf{y}\sin t) - \{(\bf{p}+a\bf{x}\cos t+b\bf{y}\sin t)\cdot\bf{m}\}\bf{m}+h\bf{m}となる。
内積の線形性より
\bf{p}-(\bf{p}\cdot\bf{m})\bf{m}+h\bf{m} + a\{\bf{x}-(\bf{x}\cdot\bf{m})\bf{m}\}\cos t + b\{\bf{y}-(\bf{y}\cdot\bf{m})\bf{m}\}\sin tである。
ここで、
\bf{p}'=\bf{p}-(\bf{p}\cdot\bf{m})\bf{m}+h\bf{m}
\bf{x}'=\bf{x}-(\bf{x}\cdot\bf{m})\bf{m}
\bf{y}'=\bf{y}-(\bf{y}\cdot\bf{m})\bf{m}
とおくと、投影した楕円E'は新たなp', x', y'を用いて \{\bf{p}'+a\bf{x}'\cos t+b\bf{y}'\sin t \: | \: 0 \leq t < \pi\}と書ける。
投影楕円E'の法線向きは投影された平面Pの法線向きmとなる。
p'は元の楕円の中心点pを単に平面Pに投影したもの等しい。
しかし、ここでの\bf{x}', \bf{y}'は必ずしも長さ1でなく、必ずしも直交しないため適切な長径向き・長径・短径を求め直す必要がある。

長径・短径向きの直交化

ここでは楕円中心点から楕円上の点までのベクトルを \bf{p}=\bf{u}\cos t+\bf{v}\sin tと表す。
uはax'に、vはby'に相当する。
長径・短径はpの長さの最大値・最小値であるから、極値となるtを求めることで長径または短径を見つけることができる。長さのままでは扱いづらいため二乗を用いる。
pの長さの二乗 |\bf{p}|^2=|\bf{u}|^2\cos^2t+2(\bf{u}\cdot\bf{v})\sin t \cos t+|\bf{v}|^2\sin^2t
をパラメータtで微分し0と等しいときのtを求める。
\frac{d|\bf{p}|^2}{dt}=|\bf{u}|^2(-2\sin t\cos t)+2(\bf{u}\cdot\bf{v})(\cos^2 t - \sin^2 t)+|\bf{v}|^2(2\sin t\cos t) = 0
変形して
(|\bf{v}|^2-|\bf{u}|^2)(\sin 2t) + 2(\bf{u}\cdot\bf{v})(\cos 2t) = 0
tについてまとめると
\tan 2t = \frac{2\bf{u}\cdot\bf{v}}{|\bf{u}|^2-|\bf{v}|^2}
t=\frac{1}{2}\tan^{-1}(\frac{2\bf{u}\cdot\bf{v}}{|\bf{u}|^2-|\bf{v}|^2})
となる。
tまたはt+\frac{\pi}{2}がそれぞれ長径または短径におけるパラメータとなる。
ここから長径向き・長径・短径を求めることで最初に示した形式に一致する表現を得ることができる。

確認

tとt+\frac{\pi}{2}でのpが直交することを確認する。
それぞれを代入して内積をとる。
\bf{p}=\bf{u}\cos t+\bf{v}\sin t
\bf{p}'=\bf{u}\cos (t+\frac{\pi}{2})+\bf{v}\sin (t+\frac{\pi}{2})=-\bf{u}\sin t+\bf{v}\cos t
\bf{p}\cdot\bf{p}'=|\bf{u}|^2(-2\sin t\cos t)+2(\bf{u}\cdot\bf{v})(\cos^2 t - \sin^2 t)+|\bf{v}|^2(2\sin t\cos t)=\frac{d|\bf{p}|^2}{dt}
ここで、 \frac{d|\bf{p}|^2}{dt}=0となるようにtを求めたのであるから、内積は0となり直交することが確認できた。

TODO

図を書く

ARC004_C

AtCoder Regular Contest #004 C問題
http://arc004.contest.atcoder.jp/tasks/arc004_3

ケアレスミス記念＆はてダ記念。

僕の解法(?)はこんなのだった。

1からNまでの平均値= \frac{N+1}{2}

足し忘れ数M= \frac{(N+1)N}{2} - N\frac{X}{Y}

0 < M \leq N, \: M \in \mathcal{N}となるNを見つける。

このままではNの範囲がわからないので不等式を変形して範囲を求める。

下限

0 < \frac{N}{2} + (\frac{1}{2} - \frac{X}{Y})

2\frac{X}{Y} - 1 < N

上限

\frac{N}{2} + (\frac{1}{2} - \frac{X}{Y}) \leq N

N \leq 2\frac{X}{Y} + 1

というわけで下限と上限が求まった。
よくみると範囲はめちゃくちゃ狭い。
ところが下限の計算をミスって O(\frac{X}{Y})の線形探索になってしまい、どこで計算を削ろうかと考えてしまいタイムアップだった。もったいない。

適応的しきい値

借りた本をスキャンしたんだけれど 真ん中のへこんだ部分がどうしても暗くなる。

本の印刷のコントラストが適度に低くて しかも本が適度に分厚いと 真ん中のへこんだ部分がより大きく暗くなるので 全体の文字がつぶれないようにがんばると どうしてもこうなる。

threshold0.png threshold0.png

そこで適応的にしきい値かけるソフトを探したんだけれど openCVに標準で関数があったので テキトーに書いてみた。 ていうか「適応的しきい値」とか「適応的二値化」ってマイナーなのかしら あんまりソフトがひっかからんかったような気がする。

適応的しきい値をかけるとこういう感じになってくれてとてもうれしい。

threshold0.png threshold1.png

参考 http://opencv.jp/sample/filter_and_color_conversion.html#threshold

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

int main (int argc, char **argv)
{
  IplImage *src_img = 0, *dst_img;
  IplImage *src_img_gray = 0;

  //args
  // inputfile outputfile blocksize thresholdparameter
  char *outfile="out.jpg";
  int block=5;
  double param=5;

  if (argc >= 2)
  {
    src_img = cvLoadImage(argv[1], CV_LOAD_IMAGE_COLOR);
  }
  if (src_img == 0)
  {
    return -1;
  }
  if (argc >= 3)
  {
    outfile=argv[2];
  }
  if(argc>=4)
  {
    int b=atoi(argv[3]);
    if(b<3)b=3;
    b/=2;
    b*=2;
    b++;
    block=b;
  }
  if(argc>=5)
  {
    double p=atof(argv[4]);
    param=p;
  }

  dst_img = cvCreateImage(cvGetSize (src_img), IPL_DEPTH_8U, 1);
  src_img_gray = cvCreateImage(cvGetSize (src_img), IPL_DEPTH_8U, 1);

  cvCvtColor(src_img, src_img_gray, CV_BGR2GRAY);

  cvAdaptiveThreshold(src_img_gray, dst_img, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY, block, param);

  cvSaveImage(outfile,dst_img);

  cvReleaseImage(&src_img);
  cvReleaseImage(&src_img_gray);
  cvReleaseImage(&dst_img);

  return 0;
}