More servicesWindows Live
HomeHotmailSpacesOneCare
 
MSN
Sign in
 
 
Spaces home  ※ 画面は開発中のものですPhotosProfileFriendsBlog Tools Explore the Spaces community

Blog

    • View next 20 entriesView last 20 entries
    September 23

    レーシック

    手術日の朝は7:30に起床することにしていたが、早めに目覚めた。
    が、全く寝られないほどではない、その程度の緊張感だった。
     
    ギリギリまで 「やめとけ、レーシック」みたいなページばかり見ていた。
    当日キャンセルは原則無理と書いてあった。
     
    • 病院に着いた
    • 呼ばれた
    • 視力検査をした
    • 高いお金を払った
    • 待機
    • 手術の説明
       術後にかけるメガネと眼帯と、説明が書いてある紙を貰う。
    • 手術フロアに行った
       今までにつけてためがねと手荷物をロッカーに入れる。
    • 問診
       「質問ある?」に対して「今は無い」(後からされても困るだろうが)
       しかし、挙げるとすれば先知れぬ不安。
    • 待機

    最初の部屋
     フラップというものを作成するらしい部屋。
     イスに座らされると機械が頭上に現れる。
     目をぐりっとあけられる、今までに無いくらい。
     そして、光が目に当たる。
     「光をみててね~」といわれる。
     ものすごく押さえつけられる感じがするとすぐに機械がうなる。
     目には光しか見えない。
     耳には秒数を数えるナースの声と機械の音。

     終わると、目の前が白くなったが、全く見えないわけではない。
     ナースに付き添われて次の部屋の前へ移動する。

    2番目の部屋
     部屋があくまで目を閉じて待たされる。
     白い世界の状態なので、放置されているような気になる。
     10分くらい待って部屋に入ると同じような手術イスがある。

     またレーザーを照射される。
     つい力が入り、「力ぬいてね~、そうそう、うごかないよ~」 といわれる。
     穴の開いたシートをかぶせられ、レーザーが数十秒当たり、洗浄されてはがされる。
     道具を使って目をなぞられて(貼り付けるように)

     「終わりです」と言われて立ち上がり、回復室に連れて行かれる。

    回復室
     イスに座って眠らずに待つ。
     作業は終わったが緊張感が残っていたので眠気とは無縁。

    最後のチェック・解散
     特に問題ないので目薬を貰ってさようなら

    帰りの状況
     涙がひどかった。おまけに鼻水も出てきた。
     ただし、痛み止めは即効性があり、半分くらい帰ったところで点したらすぐにとまった。
     一人でも帰れそうだが辛そう、というレベル。

    これから
     いちばん怖いのは近視の戻り。
     せっかく1や1.5を越えたりしても、数年後に0.7などというざんねんなこともあるらしい。
     それは体質とも言われるが、なるべく目は酷使しないようにしたい。

    すべきか
     値段: 待てば少しは良心的になるかもしれない。紹介制度など不透明な部分がなくなって。
     安全性: 安心だと思うまで待てばよい。待ちすぎると老眼のためにメリットが減る。
     勇気: 手術は怖かった。痛くないが、怖いことに変わりは無い。
          術後の痛みはたいしたこと無かった。

    August 15

    レーシックの事前検査

    手術が可能かどうか調べにいってきた。(ここまでは無料)
    結論は、可能。 ただ、消耗品であるらしい角膜の厚さが1回分+a しかないため、
    極わずかな確率での失敗時に再手術ができないかもしれないとのこと。
    視力1.0獲得率は98%で、だいたい0.6くらいあればいいと思っているので成功率は 98%+少し。
     
    この確率は、信頼すべきものだと思うが、勇気は要る

    眼鏡屋でフレームを選ぶためにも。(まだかけるつもり)
    高い特殊レンズ代を省略して安く買うためにも。
    成功すればプールサイドで躓かない。
    成功すれば免許も裸眼で取れる。
    視野が広がる(最近のメガネはレンズ小さいから。)
    シャンプーとリンスを間違えない。
    メガネをなくしても匍匐前進は不要。
    悲惨な事故(踏み潰されたり)で大事なメガネを失わずに済む。
    裸眼時の目は3でもεでもない。
     
    夢のような生活に戻れるらしい、18万くらい & 98+a% で。

    今は昔、最初メガネをかけるときは、
    それまでの数週間、寝る前にミドリンとかいう目薬を点眼したりしてメガネにならぬよう粘った。
    結局、周りの人より早くメガネをかけ、「メガネくん」とお決まりの名前をもらいながら、
    「どうしようもないのか」と、しばらくは我慢の日々だった。

    それが今日はなんという技術の進歩だ。このまま行けば不老不死も夢じゃないかと思えるほどだ。

    でもやはり少し怖いとは思う。
     

    検査は結構長かった。
    明るいときと暗いときの見え方の検査、視力検査をして、
    よくわからない目薬を点眼されて、麻酔を点眼されて、目玉をグリグリされる。
    それで最後に結果発表。

    帰りは目薬の効果でまぶしい。夕日がまぶしくてお化けになったかと思った。
    地下通路を駆使して何とか帰宅。それでもまぶしさで普段より疲れた。
    手術のときはタクシーが良いかもしれない。

    さて、残るは最後の決断である。

    August 06

    Nicovolume 0.252

    1. ちょっとした機能追加をしてみた。

    時間が来ると、下のような画面が表示され、スライダーが左右に動きます。
    手動で動かすことによって、自動調整をキャンセルすることが可能です。

    0.252

    2. IE の自動起動が不要な場合のため、インストーラで ブラウザ拡張の有無を選択できるようにしました。
    その結果 Vista 環境で .msi 形式のインストーラがうまくいかなくなってしまったので、 exe 形式も入れておきました。

    3. Vista 環境で、 互換モード を使用する必要がなくなりました。
     Vista では 新しいAPI を使うようにしました。

    * 古いバージョンを消してから新しいバージョンをインストールする必要があります。
    * インストールフォルダ内にある time.conf を一時退避させておくと同じ設定で使えると思います。

    Nicovolume の記事
    http://prog-city.spaces.live.com/Blog/cns!E5379F11D9E8BB3F!918.entry


    Vista でのマスタボリュームコントロール についての 参考文献

    http://msdn.microsoft.com/en-us/library/ms679161(VS.85).aspx
    http://msdn.microsoft.com/en-us/library/bb331828(VS.85).aspx

    August 01

    フーリエ級数

    フーリエ何とかというものがさっぱりわからないので、教科書に書いてある次のような数式をもとに計算させてみた。

    周期Tの周期関数 f(t) が下のように展開される。

    image

    OnPaint に直接適当にコードを書く。
    f(t) は、連続してない関数にした。

        public int n;
        public double func(double d) {
            return ((int)Math.Abs(d) % 50 < 25) ? 0 : 50;
        }
    
        public double func_a(double d) {
            return func(d) * Math.Cos(n * 2 * Math.PI / 50 * d);
        }
    
        public double func_b(double d)
        {
            return func(d) * Math.Sin(n * 2 * Math.PI / 50 * d);
        }
    
        private void Form1_Paint(object sender, PaintEventArgs e)
        {
            int x;
            double y = 0;
    
            Point prev = new Point(0, 0);
            for (x = 0; x < 300; x++) {
                double a0 = 1.0 / 50.0 * mymath.simpson(-25, 25, 0.001, func);
                y = a0;
                
                for (n = 0; n < 10; n++) {
                    double an = 2.0 / 50.0 * mymath.simpson(-25, 25, 0.01, func_a);
                    double bn = 2.0 / 50.0 * mymath.simpson(-25, 25, 0.01, func_b);
    
                    y += an * Math.Cos(n * 2 * Math.PI / 50 * x);
                    y += bn * Math.Sin(n * 2 * Math.PI / 50 * x);
                }
                e.Graphics.DrawLine(Pens.Black, new Point(x, 200 - (int)y), prev);
                prev = new Point(x, 200 - (int)y);
            }
        }

    * simpson関数はこちら

    すると、こんな感じに表示された。

    25行目のループ条件: n < 2
    image

    n < 5
    image

    n < 10
    image

    本来のグラフ
    image

    ループすればするほど元の四角いグラフに近づいている。
    しかし、その分時間がかかった。

    シンプソンの公式

    シンプソンの公式のプログラムを汎用的にしておく。

    シンプソンの公式は、関数を二次曲線で近似させて積分する方法。
    台形公式より正確な値を出せる、とのこと。

    関数f の、st~ed の間を、刻み幅h以内で積分する。
    hは、刻み数が奇数になってしまった場合には変化する。

        //シンプソンの公式で積分する
        public delegate double simpsonFunc(double d);
        static public double simpson(double st, double ed, double h, simpsonFunc f){
            //刻み数を決める
            int nh = (int)(Math.Abs((ed - st) / h) + 1.0);
            if (nh % 2 != 0) nh++;
            h = (ed - st) / (nh-1);
    
            int cnt;
            double total = 0;
            for (cnt = 0; cnt < nh; cnt++ )
            {
                double tmpadd = f(h * cnt);
    
                if (cnt == nh - 1)
                { //最後:そのまま
                }
                else if (cnt == 0)
                { //最初:そのまま
                }
                else if (cnt % 2 != 0)
                { //奇数
                    tmpadd *= 4;
                }
                else
                { //偶数
                    tmpadd *= 2;
                }
                total += tmpadd;
            }
            return h / 3.0 * total;
        }
    

    この前のシンプソンの公式
    http://prog-city.spaces.live.com/blog/cns!E5379F11D9E8BB3F!940.entry

    July 11

    HSP でモンテカルロ

    image

    提出用ではない。
    やっぱりアルファベットが並んでいるとプログラミング言語らしい。

    乱数をきちんと発生させているつもりなのに、十字の模様がうっすらと見えてしまう。
    乱数生成まで自分で書かなければならないかも。


    //モンテカルロ
    title "モンテカルロ?"

    //使用変数(円)
    max_in = 1000000     ; 総ドット数(入力用
    max = 0;           
    in_dot = 0;            ; 内部ドット数
    rest = 10000         ; await をいれる間隔
    visible = 0         ; 1 で、描画過程を表示する
    vx = 390 : vy = 00    ; 表示座標
    r = 200                ; 半径
    pi = 0                ; 求めるpi

    //グラフ
    graph_x = 100 : graph_y = 450
    graph_w = 600 : graph_h = 150

    randomize            ; 乱数初期化
    rand_max = 32768    ; 乱数の最大値

    //メインウィンドウ初期化
    screen 0, 800, 600
    //コントロールの配置
    objsize 200, 25
    pos 0, 25    : combox dimensions, 100, "円\n球"

    pos 120, 0    : chkbox "過程描画", visible_chk

    pos 0, 55    : mes "ドット数"
    pos 80, 55    : input max_in, 120, 25
    pos 0, 80    : mes "内部ドット数"
    line 0, 100, 100, 100
    pos 0, 102    : mes "全ドット数"
    pos 100, 90    : mes "x volume="

    //数式を読む
    pos 25, 200 : picload "eq.jpg", 1

    objsize 200, 45
    pos 0, 350        : button "Start",*start
    obj_start = stat

    //軸の描画
    line vx-3, vy, vx-3, vy+2*r
    line vx, vy+2*r + 3, vx+2*r, vy+2*r + 3
    //目盛の描画
    pos vx-20, vy : mes "1"
    line vx-9, vy, vx-3, vy

    pos vx-20, vy+r : mes "0"
    line vx-9, vy+r, vx-3, vy+r

    pos vx-20, vy+r*2 : mes "-1"
    line vx-9, vy+r*2, vx-3, vy+r*2
    line vx, vy+r*2+3, vx, vy+r*2+9

    pos vx+r-9, vy+r*2+6 : mes "0"
    line vx+r, vy+r*2+3, vx+r, vy+r*2+9
    pos vx+r*2-9, vy+r*2+6 : mes "1"
    line vx+r*2, vy+r*2+3, vx+r*2, vy+r*2+9

    //グラフ軸の描画
    gosub *init_graph
    stop
    *start
        dimension = dimensions//これ以降dimensions はない
        max = max_in;
        visible = visible_chk
        //円描画用バッファ
        buffer 2, r*2, r*2;
        cls
        if(visible == 1): gsel 0;
        pos vx, vy
        gcopy 2, 0, 0, r*2, r*2
        //円を作る繰り返し
        gosub *init_graph;
        in_dot = 0
        for i, 0, max, 1
            //必要な値の計算
            rx = rnd(rand_max) - rand_max/2;
            ry = rnd(rand_max) - rand_max/2;
            if(dimension == 1) : rz = rnd(rand_max) - rand_max/2:else:rz = 0;
            len = sqrt(rx*rx + ry*ry + rz*rz)
            px = rx * r / (rand_max / 2);
            py = ry * r / (rand_max / 2);
            //直接描画の場合、描画位置を座標に加算
            if(visible == 1) : px += vx : py += vy ;
            //色の指定
            if (len/(32768/2) < 1){
                in_dot ++
                color 100, 100, 255
            }else{
                color 0, 0, 0
            }
            //描画
            pset px + r, py+ r
            //一定間隔で待機を行い、進行状況・円周率を表示
            if(dimension == 0): vol = 4:else : vol = 8;
            pi = (((1.0*in_dot)/(i+1))*vol)
            if( i \ rest == 0){
                prev = ginfo(3)
                gsel 0
                color 100, 100, 255
                boxf 0, 0, 100, 10
                color 0, 0, 0;
                boxf 0, 0, (i * 100) / max, 10

                color 255, 255, 255
                boxf 180, 90, 300, 115
                color 0, 0, 0
                pos 180, 90
                mes ""+pi//円周率
                gsel prev
                await 10
            }
            //グラフに打点すべきタイミング
            if(max / graph_w < 1){
                if(dimension == 0):def = 3.1416:else:def = 4.1888
                prev = ginfo(3)
                gsel 0
                color 255,0,0
                pset graph_x + i, graph_y + graph_h/2 - pi*500 + def*500
                gsel prev
            }else{
                if(i\(max / graph_w) == 0){
                    if(dimension == 0):def = 3.1416:else: def = 4.1888
                    prev = ginfo(3)
                    gsel 0
                    color 255,0,0
                    pset graph_x + (1.0*i/max)*(graph_w), graph_y + graph_h/2 - pi*500 + def*500
                    gsel prev
                }
            }
        next
        gsel 0
        color 0, 0, 0
        boxf 0, 0, 100, 10
        //非表示モードではまとめて円を描画
        if(visible == 0) :     pos vx, vy:  gcopy 2, 0, 0, r*2, r*2
        stop
    *init_graph
        if(dimensions == 0):def = 3.1416:else:def = 4.1888
        prev = ginfo(3)
        gsel 0
        color 255, 255, 255
        boxf 0, graph_y, graph_x + graph_w, graph_y + graph_h
        color 0, 0, 0
        pos 30, graph_y + graph_h/2 : mes ""+def
        line 99, graph_y, 99, graph_y + graph_h
        line 99, graph_y + graph_h/2, 700, graph_y + graph_h/2

        pos 30, graph_y
        mes "" + (def + (1.0 * graph_h / 500))
        line 99, graph_y, graph_x + graph_w, graph_y
        gsel prev
        return ;

    July 09

    C# でベジエ曲線(4)

    ベジエ曲線が四角形の中にあるのか、無いのかを知るには、どうすればよいのだろうか。

    1. ベジエ曲線は、制御点の凸包の外側には飛び出さないということなので、四角形の中に凸包の領域が無ければ曲線は含まれず、有れば曲線が含まれるかもしれない。
    2. 最初か最後の制御点が四角形に含まれていれば、四角形内に曲線が含まれる。
    3. 曲線と四角形の辺に交点があれば、四角形内に曲線が含まれる。

    1の判別:
     凸包を求める必要がある。
     グラハムスキャンという方法があるらしい。

    2の判別:
     if 文を使う

    3の判別:
     辺の直線の式にベジエ曲線の式を代入し、t についての方程式を解く。
     が、次数が高くなるかもしれないことになっているので、工夫がいる。

     ニュートン法という方法で 解が求められるというので、実装してみたのだが、どんなケースでも解が求まるわけではないらしい。 難しそうなのでスルーしていた bezier clipping に手を出してみる。

    こちらのページとにらめっこして理解を試みる。
    http://nis-lab.is.s.u-tokyo.ac.jp/nis/bezclip.html

    Bezier Clipping は、最近のアルゴリズムらしい。
    曲線の制御点による凸包の  y= 0 を通る辺 の x軸との交点2つをそれぞれmin, max として、曲線の当区間を切り出し、それを再帰的に行うというもの。

    まず、 曲線の式を用意する
    bezier

    次に、直線の式を用意する
    image

    そして、解くべき式は次のようになる。

    image

    この式の各項 の f i が制御点 i の値になる。 各項での t は、 1 を 項数で均等に割ったものらしい。

    この式に対して、次のように処理する。参考にしたのはこちらのページ
    http://www-ui.is.s.u-tokyo.ac.jp/~kenshi/assignment/cg07/bezierclip.html

    1. 凸包を得る
    2. 各辺のうち、f(t)=0 を通る辺の、f(t) = 0 での t を求め、min, max とする
    3. そのような辺がなければこの部分には解がない
    4. min, max の範囲が、もともとの半分以上なら、解が複数あるとみて、曲線を半分に分割する。
      その両方の曲線を再帰し、両方の結果を戻す
    5. そうでない場合、収束しているか調べ、していれば値を返す
    6. していなければmin, max の範囲で曲線を切り出し、再帰し、結果を戻す。

    自分の場合、切り取った曲線の t の変域は 0 <= t <= 1 だったので、 再帰した処理が戻ってくるとき、(max – min)倍に縮めて最小値を加算する処理を付け加えた。

    //このオブジェクトをベジエクリッピングする
    public float[] bezierclipping() {
        //Pointf にする t は項数で均等に分けたもの
        PointF[] pts = new PointF[points.Length];
        PointF[] convex;
        for (int i = 0; i < points.Length; i++) {
            pts[i].X = (1.0f / (points.Length-1)) * i;
            pts[i].Y = points[i];
        }
    
        //凸包を得る
        convex = mymath.convex_hull(pts);
    
        //凸包のうち、ゼロを通る辺を最小・最大値とする
       float min = float.MinValue, max = float.MaxValue;
       for (int i = 0; i < convex.Length; i++)
       {
           int j = (i + 1) % convex.Length;
           if (convex[i].Y * convex[j].Y < 0)
           {
               if (min == float.MinValue)
               {
                   min = (convex[i].X * (float)Math.Abs(convex[j].Y)
                       + convex[j].X * (float)Math.Abs(convex[i].Y))
                       / ((float)Math.Abs(convex[i].Y) + (float)Math.Abs(convex[j].Y));
    
               }
               else if (max == float.MaxValue)
               {
                   max = (convex[i].X * (float)Math.Abs(convex[j].Y)
                       + convex[j].X * (float)Math.Abs(convex[i].Y))
                       / ((float)Math.Abs(convex[i].Y) + (float)Math.Abs(convex[j].Y));
                   break;
               }
           }
           else if (convex[i].Y == 0) {
               if (min == float.MinValue) min = convex[i].X;
               else if (max == float.MaxValue) { max = convex[i].X; break; }
           }
       }
        if (min == float.MinValue || max == float.MaxValue)
            return new float[] { };
    
        if(min > max){
            float tmp;
            tmp = min; min = max; max = tmp;
        }
        if (max - min < (float)Math.Abs(pts[pts.Length - 1].X - pts[0].X) * 0.5f)
        {
            //0.5以下に小さくなるのならば、解はひとつ
            if (max - min < 0.0001) //収束判別
            {
                return new float[] { (max - min) / 2 };
            }
    
            //もう一度縮めたい
            float[] tmp = cut(min, max).bezierclipping();
            for (int cnt = 0; cnt < tmp.Length; cnt++) {
                tmp[cnt] *= (max - min)/1;
                tmp[cnt] += min;
            }
            return tmp;
        }
        else 
        { 
            //0.5以上ならば解が複数ある
            bezierfunc[] cut_half = cut(0.5f);
            float[] a, b, c;
            a = cut_half[0].bezierclipping();
            b = cut_half[1].bezierclipping();
    
            c = new float[a.Length + b.Length];
            for (int i = 0; i < c.Length; i++) {
                c[i] = (i < a.Length) ? a[i]/2 : (b[i - a.Length]/2 + 0.5f);
            }
            return c;
        }
    }

    image

    July 04

    C# でベジエ曲線(3)

    計算が複雑になってきてきた。そこで、多項式の計算をクラスにしてみた。
    四則演算と微分・積分を実装した。
    ただし、 1/x の積分は出来ない。 とか、速度の問題とかいろいろとある。

    カーブ半径を取得したい。
    曲線の曲がり具合は、曲率というもので表されるらしく、
    一般に用いる半径は曲率半径といって曲率の逆数らしい。

    したがって、曲率を求めたい。が、
    http://ja.wikipedia.org/wiki/%E6%9B%B2%E7%8E%87
    わからない。

    検索して調べていると

    (x’y’’ – x’’y’ ) / √{({x’}^2 + {y’}^2)^3}

    x(t), y(t) : 曲線上の点
    x’(t), y’(t) : 接線の傾き

    こんな式を見つけたのでそのまま計算して逆数にすると、正しいらしい結果が出てきた。image

    カーブ半径に合わせて円を描いてみた。


    下がコードの一部。 tmp1 が上式の 分子、 tmp2 が分母。

           public float curvature( float t) {
                float tmp1 = -f_tangental_x.differentiate().GetValue(t) * f_tangental_y.GetValue(t) +
                                f_tangental_y.differentiate().GetValue(t) * f_tangental_x.GetValue(t);
                float tmp2 = (float)Math.Sqrt((double)powFI((powFI(f_tangental_x.GetValue(t), 2) 
                                + powFI(f_tangental_y.GetValue(t), 2)), 3));
                if (tmp2 == 0) return 0;
                return  (tmp1 / tmp2);
            }

    differentiate は多項式クラスで微分をするメソッドで、 (x^a)’ = ax^(a-1) をすべての項に対して行う。
    GetValue は多項式の変数に代入し、値を得るもの。

            //微分する
            public polynomial differentiate() {
                polynomial result = new polynomial();
    
                foreach (KeyValuePair<float, float> kpy in this.values)
                {
                    if (kpy.Key != 0)
                        result.values[kpy.Key - 1.0f] = kpy.Key * kpy.Value;
                }
    
                return result;
            }

    多項式クラスでは、各項を、キーを指数、値を係数とするリストに保存している。

     

    正確さや速度に問題が出てくるかもしれないが、機能としてはこれくらいかと思う。
    次はカーブをつなげたレールの設計とか。

    July 01

    夏が近づく / C#

    ビワを食べた。
    家族には不人気のビワ、好きなのは自分だけらしい。
    おいしかったからこれも植えてしまおうか。

    DSC00791

    トカゲを見た。 トカゲと呼んでいるがこれはカナヘビのような気がする。
    http://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%8A%E3%83%98%E3%83%93
    DSC00780 
    DSC00781 
    飼いたいが餌は虫。しかも手ごろなデュビアというムシがなんと ゴキの仲間。
    無理だ。

    C# の機能を一通り勉強。 しかし ・・・ 情報が古い。
    とりあえず気にせずに コーディング。

    June 29

    C# でベジエ曲線(2)

    無事、ベジエ曲線とその並行曲線がかけたのが前回まで。

    曲線の式が分かることで、列車を滑らかにカーブさせることができそうだが、
    まだ曲線の上を一定速度で動かすことができない。
    前回のベジエ曲線の式は、パラメータ t で表されていて、
    t を一定に増加させても 等速度で動いてくれないように見える。

    だから、 適当な距離 l において、パラメータ t はいくつなのか調べないといけない。
    また、曲線全体の距離も知りたい。

    そこで、高校の数学教科書の最後のほうのページに、「曲線の長さ」という項目があったことを思い出し、開いてみた。

    (b~a の定積分)∫√[ { f’(t) }^2  +  { g’(t) }^2 ] dt

    次数の高い多項式がルートの中に入ったものを積分するのは非常に難しそう。
    つい最近、大学で積分の数値計算を習ったことから、それで計算することにする。

    とはいっても、積分は基本的に、地道に面積を計上していくただのループ。
    なんとかならないものかという気持ちが残る。

    --------------

    積分のアルゴリズムには、台形公式とシンプソンの公式と、ほかに各点の重みや積分区間を不等間隔にする何とかという公式があるらしい。

    あるグラフを区間ごとに区切って面積を計算する時、その区間のグラフが曲線でないこととして計算すれば、台形公式になる。 それでは誤差が大きくなってしまうので、 区間上の3点を通る2次曲線の積分として考えれば、少しは誤差は小さくなる。これがシンプソンの公式である。

    他にもいい方法があると本には書いてあるが、難しくなると面倒だし、ピクセル単位であっていればいいのでシンプソン法で計算してみることにした。

    ------------

    とりあえず、 曲線の パラメータ t での長さ l が何ピクセルなのかを知りたい。

    計算のやり方は教科書のまる写しでなんとかする。
    グラフを、 (t0, l0), (t1, l1), (t2, l2),,,,,(tN, lN) と区切って、
    (h/3)(l0 + 4l1 + l2) + (h/3)(l2 + 4l3 + l4) + … と計算していく。

    これを簡単にするとこんな感じ。
    (h/3) * {
    l0 + lN
    + l偶数番 * 2
    + l奇数番 * 4
    }

    ちなみに、区間は偶数個に区切らなければならないらしい。

            //パラメータ t を 開始点からの距離 l に変換する h は積分刻み幅
            //シンプソンの公式を用いた
            static public float t2l(float t, PointF[] pt, float h)
            {
                float tcnt, prev =0.0f, tmplen, tmpadd, total = 0.0f;
                int i;
                for (tcnt = 0.0f, i = 0; tcnt <= t; i++, tcnt += h) {
                    tmplen = length(bezier.tangental_func(pt, tcnt));
     
                    if (tcnt + h > t) 
                    { //最後
                        if (i % 2 == 1)
                        { //奇数個で終わる
                            tmpadd = (prev + tmplen * 4 + length(bezier.tangental_func(pt, tcnt + h))) / 2;
                        }
                        else
                        { //偶数個で終わる
                            tmpadd = tmplen;
                        }
                    }
                    else if (tcnt < h)
                    { //最初:そのまま
                        tmpadd = tmplen;
                    }
                    else if (i % 2 != 0)
                    { //奇数
                        tmpadd = tmplen * 4;
                    }
                    else
                    { //偶数
                        tmpadd = tmplen * 2;
                    }
                    total += tmpadd;
                    prev = tmplen;
                }
                return h/3.0f * total;
            }

    積分刻み幅を任意に定められることにしたかったので、項数が奇数になってしまった場合、区間を超えて計算し、超えた区間を2で割るようにしてみた。

    そして、その逆変換

            //開始点からの距離 l を パラメータ t に変換する h は積分刻み幅
            //シンプソンの公式
            static public float l2t(float l, PointF[] pt, float h)
            {
                float tcnt, prev =0.0f, tmplen, tmpadd, total = 0.0f, lcnt, endlen = 0.0f;
                int i;
                for (tcnt = 0.0f, i = 0; tcnt <= 1.0f && endlen < l; i++, tcnt += h) {
                    tmplen = length(bezier.tangental_func(pt, tcnt));
     
                    //ここで終了した場合の長さ
                    if (i % 2 == 1)
                    { //奇数個で終わる
                        endlen = (total + (prev + tmplen * 4 + length(bezier.tangental_func(pt, tcnt + h))) / 2)
                            * (h / 3.0f);
                    }
                    else
                    { //偶数個で終わる
                        endlen = (total + tmplen)
                             * (h / 3.0f);
                    }
    
                    //加算処理
                    if (tcnt < h)
                    { //最初
                        tmpadd = tmplen;
                    }
                    else if (i % 2 != 0)
                    { //奇数
                        tmpadd = tmplen * 4;
                    }
                    else
                    { //偶数
                        tmpadd = tmplen * 2;
                    }
    
                    total += tmpadd;
                    prev = tmplen;
                }
                return tcnt - h/2.0f;
            }

    距離から t を求める場合、ループするたびにその地点での距離を知りたいので、そのつど終了した場合の距離を計算し、終了判断しなければならない。 あと、戻り値は必ず 本来の値以上になることから、 刻み幅の半分を引いてみた。

    image

    t = 0.5 の位置を距離に変換し、戻してみた。
    大体あっているが、実用的なのか微妙なところ。

     

    それより、コードがなんか汚い。 言い訳としては、ご丁寧に中かっこをエディタが改行してくれるということ。
    IntelliSense も C++ に比べるとおせっかいに働いてくれる。

    C# でベジエ曲線

    C# を始めた。
    Visual studio にはいくつかの言語が入っていて、 VB とかもできるのだが 手は伸びていなかった。
    C を使っていると 効率が悪いということにいい加減気づくべきではないかと思い、 C# を始めてみようと思ったわけです。

    #include の順番とか、 定義の順番だとか、ヘッダファイルとか、悩まされることが少ない。
    メモリの確保とかも考えるべきことがない。

    C++はC++で深い。 Windows の深いところまで勉強するには C++ が必須。
    でも、普通のソフトを作るのに手間をかけすぎていては困る。器用に使い分けてみたい。

    -----------------

    C# で作りたいものは、ちょっとしたゲームのようなもの。
    鉄道系のシミュレーションが好きなので、それの簡易版っぽいものを考えている。
    それで、そのレールをベジエ曲線にしてその上をガタゴト走らせようと考えている。

    ただ、単純にベジエ曲線を 便利なライブラリにおそらく入っているであろう DrawBezier などという関数でビーッと引くだけではすまないであろう。その線に沿って、列車が動かなければならない。

    となると、数式でその座標を表さなければならない。
    なので、数式とにらめっこする必要が生じてくる。

    ----------------

    調べてみると、
    http://ja.wikipedia.org/wiki/%E3%83%99%E3%82%B8%E3%82%A7%E6%9B%B2%E7%B7%9A
    こんな感じの数式であることがわかった。

    ∑ Bi Jni(t)

    Jni(t) = (n i) t^i (1-t)^(n-i)

    (n i) という部分をベクトルかと思ってずっと考えていたらこれは nCi (組み合わせ)ということが判明。
    http://ja.wikipedia.org/wiki/%E4%BA%8C%E9%A0%85%E5%AE%9A%E7%90%86
    (n i) というのは二項係数というらしい。最初からそう書いてあればいいのに。

    とりあえず線の式はわかった。
    パラメータ t での座標を返す関数を書いてみた。制御点の数は可変にしておいた。

            static public PointF curve(PointF[] pt, float t) {
                PointF total = new PointF(0.0f, 0.0f);
                int n = pt.Length - 1;
                int i;
    
                for (i = 0; i < pt.Length; i++) {
                    total.X += pt[i].X * mCn(n, i) * (float)Math.Pow((double)t, i) * (float)Math.Pow(1.0 - t, (double)n - i);
                    total.Y += pt[i].Y * mCn(n, i) * (float)Math.Pow((double)t, i) * (float)Math.Pow(1.0 - t, (double)n - i);
                }
    
                return total;
            }

    でも1本だけでは足りない。平行にもう1本ほしい。レールだから。

    曲線の接線に対して直角に、一定距離進んだところをつないで線にすれば平行になりそう。
    接線を求めるには式を微分すればよさそう。

    Bi と(n i) を定数として、 { (t^i) * (1-t)^(n-i)} ’  の部分を考える。
    単純に考えると
    (t^i)’ * (1-t)^(n-i)  +  t^i * {(1-t)^(n-i)}’ 積の微分
    it^(i-1) * (1-t)^(n-i)  +  t^i * (n-i)(1-t)^(n-i-1)(-1) 合成関数

    よって ∑ Bi (n i) { it^(i-1) * (1-t)^(n-i)  +  t^i * (n-i)(1-t)^(n-i-1)(-1) }

    という感じになりそうなのだが、このままコードにすると失敗した。
    0乗のときは1なので微分するとそのまま消えないといけないから。

            //ベジエ曲線の、接線を表す関数
            public static SizeF tangental_func(PointF[] pt, float t)
            {
                PointF total = new PointF(0.0f, 0.0f);
                int i;
                int n = pt.Length - 1;
    
                for (i = 0; i < pt.Length; i++)
                {
                    float d1 = (i == 0) ? 0.0f : (float)i * powFI(t, i - 1);
                    float d2 = (n - i == 0) ? 0.0f : (float)(n - i) * powFI(1.0f - t, n - i - 1) * (-1);
    
                    total.X += pt[i].X * mCn(n, i) *
                    (
                        d1 * powFI(1.0f - t, n - i) + powFI(t, i) * d2
                    );
    
    
                    total.Y += (pt[i].Y * mCn(n, i) *
                    (
                        d1 * powFI(1.0f - t, n - i) + powFI(t, i) * d2
                    ));
                }
                return new SizeF(total.X, total.Y);
            }

    * powFI は float 型をキャストせずに べき乗できるようにした単純な関数。
    * d1 と d2 がそれぞれ t^i と (1-t)^(n-i) を微分したものをいれる変数。

    この戻り値をそのまま描画したらうまく接線が描画された。
    しかし、場所によって長さが変わってしまった。 曲線の パラメータ t の 1あたりの長さになってしまっているかららしい。

    よって、このx, y に同じ値 l をかけて 長さを1にする必要がある。
    そうしないと並行曲線が描けない。

    (lx)^2 + (ly)^2 = 1^2
    l^2 = 1 / (x^2 + y^2)
    l = 1 / √(x^2 + y^2)

    となり、(x / √(x^2 + y^2) , y / √(x^2 + y^2) ) が長さ1ピクセル分のベクトルになる。
    英語のwikipedia に Parallel curve という項目があって、これを見ると何となく正しいような気がする。
    http://en.wikipedia.org/wiki/Parallel_curve

    この値に必要な長さをかけて使えばよい。

    ----------

    接線を直角に曲げて、一定距離の地点を線でつないで行けば平行な線が引けそうである。
    接線のベクトルが (x, y) なら、 その直角は ( y, –x )。

    で、レール(道路)の幅だけ掛け算する。
    これで平行線のグラフもわかった。

    これらを使って適当に線を引く

    image

    はい、完成。

    ところで、 bezier 曲線は、 ベジエと発音するのがどちらかというと原語に近いとのこと。
    少し早くしゃべるだけで区別がつかなくなるような気もするが。

    June 17

    Visual studio コマンドウィンドウ

    先生は 「統合開発環境を使えば簡単ですが、 何かと応用が利かなくなったりするのでコマンドも覚えましょう。」 といいます。
     
    今日、 Visual Studio の画面上に 「コマンドウィンドウ」 を発見したので、適当に dir と打ち込んで見たところ、
     
    "コマンド "dir" は有効ではありません。"
    との反応。
     
    調べてみたところ、こんなコマンドを入力するらしいです。
     
    Visual Studio コマンドの定義済みのエイリアス
     
    ショートカットキーで十分そうなコマンドが多いですが、
    使いそうなものをメモ。
     
    • ? x : xの内容を表示する
    • gotoln : 指定行に移動(ステータスバーの行数をダブルクリックしてもOK)
    • lcase : エディタの選択範囲を小文字にする
    • ucase : エディタの選択範囲を大文字にする
    • nav : ブラウザタブを開く
    • shell x : シェルで x を実行
    • wordwrap : 右端で折り返す(折り返さない)

    エイリアスの作成方法
    http://msdn.microsoft.com/ja-jp/library/xasxzd71(VS.80).aspx

    June 04

    simutrans200806

    Simutrans 、フリーのシミュレーションゲーム。

    町にバスを走らせ、町外れに駅を作り、町同士をつないでネットワークを広げていく。
    最初はコストのかからない交通機関を多く使用してネットワークの拡大に資金を注ぎ、
    拡大した後は儲けすぎて楽になるので、想像力を働かせて遊ぶ。

    車両自体、小さくて識別しにくいが、アドオンが多く公開されている。
    東武の車両を一式入れてみた。
    地名は変えられるが、デフォルトではランダムに割り当てられる。

    下図 : 中央の駅で通過待ちをしている。
    信号を駆使して、追い越させるようにする。ちょっとしたパズルっぽい。
    simscr08


    常磐線の中距離電車は10両編成が多い。
    反面、緑色の快速は15両編成が多い。
    現状、明らかに乗車率が違う。

    短距離の電車と同じ停車駅にとまる上に、中距離を運転するのだから、同じ定員なら乗車率が高くなることは想像できる。しかし、なぜ短距離電車の編成が長く、中距離が短いのか。

    ・・・グリーン車に乗ってもらいたいからか。

    May 28

    Nicovolume

    Nicovolume は、あらかじめ指定した時刻に、マスターボリュームを自動的に調整するソフトウェアです。 
    特定のURL を Internet Explorer で閲覧していた場合 自動的に起動するので、起動し忘れることがありません。
    特に、時報を聞いたり聞かなかったりするのに便利だと思います。
     
    個人的な使用感(必ずしも得られるわけではありません。)
    1. 音量がゆっくり下がる
    2. 何事か!?と思う
    3. 再生が止まる
    4. あっ、そうかと思う
    5. 時報が流れる

    ありがたいことに紹介されました。
    http://www.forest.impress.co.jp/article/2008/06/09/okiniiri.html

    ----------
    バージョン 0.252

    * Vista でもそのまま使えるようになりました。
    * 音量調整時、手動調整も可能なダイアログでお知らせします。

    0.252

     

    -------------
    バージョン 0.25
    [対応OS]
     Vista : インストールディレクトリ\nicovolume.exe を互換モード(XP sp2) に変更する必要があります
     XP : そのまま動きます。
     
    [対応ブラウザ]
     IE7
     IE6
     
    May 27

    携帯電話不要?

    小中学生に携帯電話持たせないで…教育再生懇が第1次報告(Yahooから)
    http://headlines.yahoo.co.jp/hl?a=20080526-00000038-yom-pol

    政府の教育再生懇談会(座長・安西祐一郎慶応義塾長)は26日夕、首相官邸で会合を開き、第1次報告を福田首相に提出した。
    子供を有害情報から守るため、小中学生が携帯電話を持つことがないように関係者に協力を促している。また、英語教育の強化を掲げ、国に小学校3年から英語を必修化するように求めた。

    持たなければいいという考えになぜ行き着くのかと疑問に思う。

    いずれは自分の身は自分で守らねばならない時が来るわけで、それを先延ばしにしても何の解決にもならない。
    など、いろいろ思うことはあるのですが、
     
    "携帯電話会社に対して通話やGPS(全地球測位システム)機能に限定した機種"
     
    これについて、 
    自分なら、頻繁に家に置き忘れると思う。
    電話しかできない上に監視までされるという非常につまらない端末、忘れてもきっと気づかない。
     
    臭いものにふたをして、監視だけしようとしても、そううまくはいかないのではないかと思う。
     
    楽しいものだからちゃんと持っていてくれるわけで、
    携帯会社には安全で楽しい携帯を作ってもらいたい。
     
     
    May 15

    アクエリアスゼロ

    今日は体調がよく、途中で死にそうになることなくゆっくりと走りきることが出来た。
    帰ってからアクエリアス ゼロという新製品を飲んでみた。
     
    感想は、
    • 透明
    • さわやかな甘酸っぱさ

    普通のアクエリアスとは別の味と考えてよい。
    これならがぶ飲みしても気持ち悪くないかもしれない。

    製品サイト
    http://www.aquarius-sports.jp/zero/index.html

    ちゃんとしたサイトのレビュー
    「見た目はまるで水、カロリーゼロの「アクエリアス・ゼロ」試飲レビュー」(GIGAZINE)
     
     
    May 07

    ゴランノス・ポンサー株式会社

    勝手にググって笑ってる。
     
    ゴランノスポンサー
     
    ついでに小売店 「字幕スーパー」
     
    この世界はゴランノスに 支配されているといっても過言ではないかも。
     
    株式会社 提供
     株式会社 提供 は、日本における民間放送のほぼ全ての番組で社名表示を行っている。その際、音声で「ゴランノス・ポンサーの提供」という表現がなされ、提供が ゴランノス・ポンサーグループに属していることを強調する形となっている。
     提供のロゴは必ず、画面の最も上に、全てのページで表示される。