↑            
  | yc         e の符号と dx,dyの符号でxc,ycの上下左右の関係が判定可能
  |＼          この原理は、内外判定等に応用出来ます。
  |φ＼        
  |    Ｘ      
  |rl／  ＼    
  |／φ    ＼xc 
  0+----------+-->x

求め方
     rl= xc * cosφ  ここで cosφ= yc/√(xc^2+yc^2)
       = xc * yc / √(xc^2+yc^2)
       =e^2 / (dy*dx)/√(e^2/dx^2+e^2/dy^2)
       =e   /  √(dx^2+dy^2)

応用

線分が交わるか？

片方の線分が水平で始点が原点になるよう回転変換＋平行移動し、も
う一方の線分とｘ軸がどこで交わるかを求めて判定する方法
まず全体の平行移動によって
片方の線分の始点を原点として 原点からsx,sy迄、 もう一方が
　　　　点rx,ry nx,nyとします。 sx,syをベクトル積(回転x拡大)すると sx,syは
			sx'=sx*sx+sy*sy 
			sy'=sx*sy-sy*sx =0 つまり sは原点からsx'迄の水平線
もう一方のy座標は
			ry'=rx*sy-ry*sx
			ny'=nx*sy-ny*sx
		ここで ry',ny'の符号が違うという条件でフルイ別けします
		   [例:(ry'*ny')>0 なら交差しないとする]
		この条件を満たした場合(つまりx軸と交差するなら)
もう一方のx座標も同じように変形して  ( dx = nx-rx ;  dy= ny-ry )
 	  rx'=rx*sx+ry*sy 
 	  nx'=nx*sx+ny*sy 
          dx' = nx'-rx'=dx*sx+dy*sy;
          dy' = ny'-ry'=dx*sy-dy*sx; 
x軸 と交わる点 xc= (rx'*dy'-ry'*dx')/(dy') =rx'-ry'*dx'/dy'が
0<xc<sx'なら交わると判定されます 
         xc=(rx*sx+ry*sy)-(rx*sy-ry*sx)*(dx*sx+dy*sy)/(dx*sy-dy*sx)
         xc=(sx^2+sy^2)(ry*dx-rx*dy)/(dx*sy-dy*sx)
         xc=sx'(ry*dx-rx*dy)/(dx*sy-dy*sx) 

         結局、e0=(ry*dx-rx*dy)とe2=(dx*sy-dy*sx)  を求めて
         0<e0/e2で e0/e2 < 1なら 交差という判断をする事になります 
          0 < e2 なら   e0<e2 および 0<e0
          0 > e2 なら e0>e2   および 0>e0

と符号だけの判断で済みます(サンプル)
	このe0/e2で交点の位置の線分上の比を得られますから、 e0.e2から
	sx*e0/e2 sy*e0/e2 と元座標での交点の位置が得られます

これも色々あるのですが多いのは
　　　　１）線分の始点終点で出来る長方形が交わるか？
　　　　２）線分が収まる円を描いて円同士が交わるか？
　　　　３）片方の端点から線分との距離２組みを求めて短い方
の３つくらいかな？
	１）は判断を並べれば良いだけ しかし、分岐が入るとCPUは遅いので、
	   アセンブラレベルで大小比較結果(CF)を集めてまとめて比較するような工夫が必要です。
	
	２）円同士が交わるかは、中心同士の距離を求めて双方の円の半径の
		和(線分の長さ合計/2)との大小比較します。 hypotenuse
		を荒く、誤差分を安全側に判定れば(分岐予測のミスでペナルティを払うCPUでは)
		1)より高速です。
	３）は線分の交差判定の初段と実は同じ計算をします。
		これ使うなら素直に交差判定した方がマシ

	

GUIの為に、必要な処理です。

マウス座標をp0 線分が点p1,p2を通るとして 
rx:=p1.x-p0.x , ry:=p1.y-p0.y
dx:=p2.x-p1.x , dy:=p2.y-p1.y
まず p1,p2で指定される図形の内側かどうかを判定し、それに
　　　合格すれば 直線迄の距離が十分小さい時に線上と判定する

１）４角形の箱内に入るかどうかを判定する方法は
if dx>0 then begin {箱になければ次の判定に飛ぶ}
     if rx<0  then continue
     if rx>dx then continue;
 end else begin
     if rx>0 then continue;
     if rx<dx then continue;
 end;
if dy>0 then begin
     if ry<0  then continue;
     if ry>dy then continue;
 end else begin
     if ry>0 then continue;
     if ry<dy then continue;
 end;
   ２）中心を通る円内に入るかどうかで判定する方法は、
       (dx2+dy2)/4> | (p2+p1)/2-p02|2
       (dx2 + dy2) > ( (p2.x+p1.x-2*p0.x)2+(p2.y+p1.y-p0.y*2)2 )
       (dx2 + dy2) > ( (rx*2-dx)2+(ry*2-dy)2 )
    計算量は多いが分岐が少ないので　最近のCPUでは検討に値する
    ただし２乗を使っているので、整数演算では桁溢れに注意

直線迄の距離は (rx*dy-ry*dx)/√(dx2+dy2) この絶対値が1〜2
なら線上であるとする。hypotenuseの高速近似
　　　 dlen=hypot(dx,dy)
if abs(rx*dy-ry*dx)<dlen*2 then 線上

 その他に直接線分と点との距離　を求める方法もあります

　　点 q を指定し 点p配列で指定した多角形の内側かどうか？
     方式１ 交差回数法: 点qから直線を引き、多角形と何度交わるか
     方式２ 巻付判定法: 多角形をゴムに換えたら点qに巻付くか？
                        ＝点qから見て多角形の線分を順に追いかけてゆくと１回転以上するかどうか
           

　　全く違う種類に見えるこの２つの方法は、効率を追いかけるとプログラムは
　　殆ど相似になる。

    どちらも 多角形pを-qだけ平行移動し、原点を基準に考える
     方式１ 交差回数法:y軸と交わった時の符号が正ならインクリメント
          線分が y軸と交わるかどうかは 線分の始点終点のx座標の符号が違うかどうか
         その座標は   yc=-(rx*dy-ry*dx)/dx
          符号だけ見れば良いので e=(rx*dy-ry*dx)と dxの符号だけで判断
          つまりdxとeの符号が同じならinc

     方式２ 巻付判定法:
         y軸を左から右に移動して交わった時、上ならインクリメント下なら
　　　　 デクリメントとする。右->左ならその反対とすれば１回転すれば±２
         になり 巻き付かない時はゼロとなる。

           y軸と交わるかどうかは 線分の始点終点のx座標の符号が違うかどうか
　　　　  で上と同じ、
          その座標は yc=-(rx*dy-ry*dx)/dx
          左から右に移動してるかどうかは dxの符号を見る
          上か下かは ycの符号を見る dxの符号で反対にするという事は
          e=rx*dy-ry*dx の符号だけでインクリメント・デクリメントす
　　　　　れば良い事が判る。

 方式２ 巻付判定法には色々な実現方法が考えられる
　　　　　１）素直に頂点の角度差を求めて(atan2で)それを加算する
               ニフテイ HGH02363 渡辺 義則 さんのアイデア
 
              浮動小数点で座標を持っていれば atan2は比較的高速に
　　　　　　　計算出来るのでそう遅くはない。atan2に精度は要らない
　　　　　　　ので高速なのを自分で作れば良い (整数版の高速なのはここ)

　　　　　２）平面をqを中心として上下左右で４事象に分割し a b c dとする 
               a->b b->c c->d d->aと移動すれば+1
               a->d b->a c->b d->cと移動すれば-1
               a<->c b<->dと移動すれば どっち側を移動したかを交点計算
　　　　　　　　して+2 または-2する
              余計に複雑にしたようだが　４角形が多い場合には +1 -1され
　　　　　　　るだけで交点計算されないので　高速になるメリットがある。


交差回数法は Ray-Intersection Method と呼ぶらしい

cY+sX+d=0       s*s+c*c=1;と表現した時

1)原点から直線迄の距離は |d|　その位置は(-s*d,-c*d)
2)点(X1,Y1)からの距離rは r=|cY1+sX1+d|
3)点(X1,Y1)に一番近い点 {c(cX1-sY1)-s*d,-s(cX1-sY1)-c*d}
4)媒介変数表示すると
        X= ct-s*d ;
        Y=-st-c*d ;
5)X'=x-dx,Y'=y-dyと座標変換した時 cY'+sX'+d+(cdy+sdx)=0
6)直線がX 方向に平行になるよう 回転させたい。
   Y'+d=0;にしたい場合
       X'= cx-sy      x= cX'+sY'
       Y'= sx+cy      y=-sX'+cY'
   Y'-d=0;にしたい場合
       X'=-cx+sy      x=-cX'-sY'
       Y'=-sx-cy      y= sX'-cY'

  ２直線の交点

c1*y+s1*x+d1=0
c2*y+s2*x+d2=0

x=( d1*c2-d2*c1)/q
y=(-d1*s2+d2*s1)/q
q=  c1*s2-c2*s1

s=(a+b+c)/2とおいて
S=sqrt(s(s-a)(s-b)(s-c)) または
=sqrt(((aabb+aacc+bbcc)-(aaaa+bbbb+cccc)/2)

サンプルリスト

----------------------------------------------------------------------
多角形の内と外の判定 

 ３角形を含む任意の多角形の内外判定

 点 x,y と n多角形 ptを与えて 点が内側か外側かを判定
function PolygonExterior(x,y,n:integer;pt:array of TPoint):boolean;
var i,ct,nx,dx,dy,rx,ry:integer;
begin
pt[n]:=pt[0]; {始点は終点と等しいとします}
ct :=0;
 for i:=0 to n-1 do begin
  rx:=x-pt[i].x;
  nx:=x-pt[i+1].x;
   if ((rx<0)and(nx>=0)) or ((rx>=0) and (nx<0))  then begin
    ry:=y-pt[i].y;
    dx:= pt[i+1].x- pt[i].x;
    dy:= pt[i+1].y- pt[i].y;
    if longint(rx)*longint(dy)<longint(ry)*longint(dx)
       then Inc (ct) else Dec(ct);
   end;
  end;
  Result:=ct=0;
end;

アルゴリズム解説:
  点を通る直線を引いて何度線分が横切るかを調べる方法もありますが、今回
  は、点の位置に立って多角形データを順に眺めていって自分の体が１回転す
  るかどうかを調べます。自分が内側に立っていれば体は1回転しますね？
    詳しく角度を調べると遅いので、180度単位に横切るかどうかを調べます
   自分の体を回す代わりに、多角形の線分が自分より上下右左にあるかで
   上にあって左から->右に変化すれば +1 右から左なら-1とします。
   下にあって左から->右に変化すれば -1 右から左なら+1とします。
                              A -> B -> C -> D -> E -> F -> G -> H ->A  計
  A-------------------B   P0    +1   0    0    0    0    0   +1    0    2
  |                   |   P1    +1   0    0    0    -1   0   0     0    0
  |   F-----------E   |   P2    +1   0   +1    0    0    0   0     0    2
  |P0 | P1        |P2 |
  |   |           |   |
  H---G           |   |
                  D---C
  左右の変化を判定してる部分が
    if ((rx<0)and(nx>=0)) or ((rx>=0) and (nx<0))  then  ....
  上か下かを判定してる部分が
    if longint(rx)*longint(dy)<longint(ry)*longint(dx) then ...

コメント:
 整数演算だけで判定していますから、TPointのbit数*2<=LongIntのビット数 
 である必要があります。
 多角形上に完全に一致した場合どうなるかはご自分で調べてみて下さい。

使用例:
  マウスを動かしながらクリックすると多角形が描かれ、クリックせずに
  動かすとその内側が赤く塗られます。

type
  TForm1 = class(TForm)
    procedure FormPaint(Sender: TObject);
    procedure FormMouseDown(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure FormMouseMove(Sender: TObject; Shift: TShiftState; X,
      Y: Integer);
  private
    { Private 宣言 }
  public
    { Public 宣言 }
  end;

var
  Form1: TForm1;
  n:integer;
  pt :array [0..101] of TPoint;

implementation

{$R *.DFM}

procedure TForm1.FormPaint(Sender: TObject);
var i:integer;
begin
if n>0 then   begin
 Form1.Canvas.MoveTo (pt[n-1].x,pt[n-1].y);
for i:=0 to n-1 do
 with pt[i] do Form1.Canvas.LineTo (x,y);
 end;
end;

procedure TForm1.FormMouseDown(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
if n <99 then begin
  pt[n].x:=X;
  pt[n].y:=Y;
  inc(n);
  end;
 Invalidate;
end;

procedure TForm1.FormMouseMove(Sender: TObject; Shift: TShiftState; X,
  Y: Integer);
begin
if not PolygonExterior(x,y,n,pt) then  Canvas.Pixels[x,y]:=clRed;
end;
end.


-----------------------------------------------------------------------
線分の上かどうか（多角形)
  点 x,y と n多角形 ptを与えて 
  多角形の内側かどうかと同時に使う事が多い、指定した線分の上かどうか
  の判定方法です。帰り値が -1 なら線分の上ではないです。でなければ、
  帰り値の配列番号上の線分の上に居ます。

function PolygonNear(x,y,n:integer;pt:array of TPoint):Integer;
var i,x0,y0,x1,y1:integer;
var rx,ry,dx,dy,dlen,slen:LongInt; //D4ならInt64として下さい
begin
pt[n]:=pt[0];
 for i:=0 to n-1 do begin
  with pt[i  ] do begin  x0:=x;  y0:=y;  end;
    rx:= x -x0;
    ry:= y -y0;
  with pt[i+1] do begin  x1:=x;  y1:=y;  end;
    dx:= x1-x0;
    dy:= y1-y0;
     if dx>0 then  begin   {クリッピン判断。線分の入る箱にあるかどうか}
         if rx<0  then continue; {箱になければ次の判定に飛ぶ}
         if rx>dx then continue;
       end else begin
         if rx>0  then continue;
         if rx<dx then continue;
       end;
      if dy>0 then begin
         if ry<0  then continue;
         if ry>dy then continue;
       end else begin
         if ry>0  then continue;
         if ry<dy then continue;
       end;
    dlen:=abs(rx*dy-ry*dx)    ;
     dx:=abs(dx);  dy:=abs(dy);
     if dx>dy then slen:=dx*2+dy
            else   slen:=dy*2+dx;  {距離は約２ピクセルに調整しています}
      if dlen < slen then    begin  Result:=i;  exit; end; {発見した}
   end;
  Result:= -1;                                             {なかった}
 end;

「解説」 クリッピング判定を先にしてから線分迄の距離を計算しています。
   割算やルートを実際に計算すると遅いので
   dlen=線分迄の距離*線分の長さ として 
   slen=線分の長さx2=2*√(dx*dx+dy*dy);の近似値
   を求めて比較しています。
   データ精度とピクセルとが違う場合は、slenをその比率倍する必要がある
   のは当然です。

「使用例」他は前回と同じで MouseMoveだけ変更

var oldcc:integer; //
procedure TForm1.FormMouseMove(Sender: TObject; Shift: TShiftState; X,
  Y: Integer);
  var cc,typ:integer;
begin
  if n>0 then
   begin
    cc :=  PolygonNear(x,y,n,pt);
   if oldcc<>cc then begin      {判定結果が前回と違ってたら}
     Canvas.Pen.Color:=clRed;   {線を赤で塗る}
     if cc>=0 then begin
      with pt[cc  ] do Canvas.MoveTo (x,y);
      with pt[cc+1] do Canvas.LineTo (x,y);
     end;
     Canvas.Pen.Color:=clBlack; {前の線を黒で塗って戻す}
     if Oldcc>=0 then begin
      with pt[Oldcc  ] do Canvas.MoveTo (x,y);
      with pt[Oldcc+1] do Canvas.LineTo (x,y);
     end;
    Oldcc:=cc;
   end;
  end;
end;

------------------------------------------------------------
Int32x32To64のインラインasm版
 32bitの掛算をして64bitの結果が欲しい場合、DELPHI4 で拡張された
 INT64に先に変換してから掛算させると64bitの掛算をするために非常に
 効率が悪く速度が落ちるようです。
 そこで WinAPIの同名の関数をインラインアセンブラで書いてみました。


function Int32x32To64(X, Y: Integer): int64;
asm
  MOV     EAX,X  {register呼び出し規約の場合はXはEAXなのでホントは不要}
  IMUL    Y
end;


 また、D3を使っていると前回掲示した内外判定や線上判定はTPointを16bit以
 上で使っている場合には失敗します。といってINT64はありませんから 以下の
 ような関数を書いてみました。
 以下はD1では使えません D3でテストしています。

{掛算(32x32=64)して引算(64)しその結果が正かゼロか負なら1,0,-1を返す}
function MulSubSGN(a,b,c,d:integer):Integer;register;{-a*b + c*d}
asm
          mov eax,a     {a=eaxなのでホントは不要}
          IMUL b         {b=edx}

          xchg c,eax     {c=ecx}
          xchg d,edx     {dはregisterでもメモリ}
          IMUL edx
          sub  eax,c     {a*b-c*dに見えるが交換してるので注意}
          sbb  edx,d

          je   @pltest  {edxがゼロなら 下位テスト}
          jge  @set1    {edxが正なら 1}
          mov  eax,-1   {edxが負なら-1}
          jmp @exit
@SET1:    mov  eax,1
          jmp @exit
@pltest:
          cmp eax,0
          jne @SET1
@exit:
end;

{ (-a*b+c*d)/√(b^2+d^2) を計算するのだが √(b^2+d^2)をま
  ともに計算すると遅いので 概算して求めている }

function RoughDistance(a,b,c,d:integer):Integer;stdcall;{ (-a*b+c*d)/√(b^2+d^2)}
asm
 mov  eax,a
 IMUL b
 mov  a,eax
 mov  ecx,edx

 mov  eax,c
 IMUL     d
 sub eax,a
 sbb edx,ecx
 mov a,eax
 mov c,edx

 {b^2+d^2}
 mov eax,b
 mov edx,d
 test eax,eax
 jns @absskp1
 neg  eax
@absskp1:
 test edx,edx
 jns @absskp2
 neg  edx
@absskp2:
 cmp eax,edx
 jge @swapSkp
 xchg eax,edx
@swapSkp:
 shr edx,1
 add eax,edx
 je @exit  {ゼロなら割算が出来ない}
 mov ecx,eax
 {割算}
 mov eax,a
 mov edx,c
 idiv ecx
@exit:
end;
----------------------------------------------------------
MulSubSGN は 内外判定 #786
    if rx*dy>ry*rdx  then inc (ct) else dec(ct);
 を
     if MulSubSGN(rx,dy,ry,rdx)< 0  then inc (ct) else dec(ct);
  と置き換えて使う事が出来ます。
  

RoughDistance は rx,ry 先の dx,dyの線分までの距離を概算で求める演算
 線の上の判定 #788 の
    dlen:=abs(rx*dy-ry*dx)    ;
     dx:=abs(dx);  dy:=abs(dy);
     if dx>dy then slen:=dx*2+dy
            else   slen:=dy*2+dx;  {距離は約２ピクセルに調整しています}
      if dlen < slen then    begin  Result:=i;  exit; end; {発見した}

  以上の５行を
   if abs(RoughDistance(rx,dy,ry,dx))<2 then
     begin  Result:=i;  exit; end; {発見した}
 と 置き換えて使います。

線同士が交わるかの判定

function xchk(x0,y0, sx,sy,rx,ry,nx,ny:double):boolean ;
var dx,dy,e0,e1,e2:double;
var ryq,nyq:double;
begin
 sx:=sx-x0; rx:=rx-x0;nx:=nx-x0;
 sy:=sy-y0; ry:=ry-y0;ny:=ny-y0;
  ryq:=rx*sy-ry*sx;
  nyq:=nx*sy-ny*sx;
 if (ryq*nyq)>0  then {同じ符号なら交差しない}
  begin
    result:=False;
    exit;
  end;

 dx:=nx-rx; dy:=ny-ry;
 e0:=(ry*dx-rx*dy);e2:=(sy*dx-sx*dy);

 if 0 < e2 then result:=( e0<e2) and (0<e0)
           else result:=( e0>e2) and (0>e0);
end;

http://www.and.or.jp/~ikeuchi/makedemo/effect/math/hypot/index.htm

ABS(a)は、多くのコンパイラで組込みマクロによって
   mov   eax,a 
   CDQ           //eaxが負数ならedx=-1 
   xor   eax,edx //eaxが負数なら全部のビットを反転 
   sub   eax,edx //さらに+1すれば完全な負数になる 
 と展開されるので、内部で分岐が無く高速です。
 
 この方式を利用して インラインアセンブラを使うと分岐無しで書くことが出来ます。
最後の２行(sarとadd)を取り除くと、絶対値の大きな方を選ぶ事になる。
 --------------- 誤差 最悪 12% の整数 高速 hypot(a,b) ------
function ihypot(a,b:integer):Cardinal;
asm
   push      b   //delphiではレジスタ渡しなので
   mov   eax,a   // この行はたぶん不要
   CDQ           //eaxが負数ならedx=-1 
   xor   eax,edx //eaxが負数なら全部のビットを反転 
   sub   eax,edx //さらに+1すれば完全な負数になる 
   mov   ecx,eax

   pop   eax  
   CDQ           //eaxが負数ならedx=-1 
   xor   eax,edx //eaxが負数なら全部のビットを反転
   sub   eax,edx //さらに+1すれば完全な負数になる 
   mov   edx,eax

   sub   eax,ecx
   rcr   eax,1
   sar   eax,31  //eaxが負数ならedx=-1
   xor   edx,ecx
   and   eax,edx
   xor   edx,ecx //edxを元に戻す 
   xor   ecx,eax //小さい方がecxになる
   xor   eax,edx //大きい方がeaxになる
   
   shr   ecx,1   //小さい方を半分にして
   add   eax,ecx
end;

---------- これを Cで書くと ----
ihypot(int a,int b)
{
 a=abs(a);
 b=abs(b);
   { register w=( ( a-b ) >> 31 ) & ( a^b ); /* a>bなら w=0  */
     return   ( a^w )+( b^w )/2;
   } 
}

//精度がもう少し必要な場合用に0.2%精度のルーチン
// 計算中の桁溢れに注意 24bit以下なら安全
function Int32x32Div(X, Y, Z: Integer): Integer;
asm
  MOV     EAX,X  {register呼び出し規約の場合はXはEAXなのでホントは不要}
  IMUL    Y
  IDIV    Z
end;

function iHypot1(a,b:Integer):Integer;
var w:integer;
begin
 a:=abs(a);
 b:=abs(b);
 Result:=0;
 if a<b then  begin w:=a;a:=b;b:=w;end;
 if b*2< a then
  begin
  Result:=a+Int32x32Div(b,b,a)*121 div 256;
  exit;
  end;
  w:=(b-a div 2);
  Result:= ( 143*a + 62*w ) div 128 + (Int32x32Div(w,w,a)*7) div 32 ; 
end;

//同じく 1% 程度の精度 (GUI 処理ならこれで十分)
function iHypot2(a,b:Integer):Integer;
var w:integer;
begin
  a:=abs(a);
  b:=abs(b);
  if a<b then  begin w:=a;a:=b;b:=w;end;
  Result:=a+ Int32x32Div(b,b,a)*55 div 128;
end;

　線分と点の距離    
  直線迄の距離は (rx*dy-ry*dx)/√(dx2+dy2) 　この時の直線上に投影した点の線分始点からの距離は 
　　　√( (rx2+ry2)-(rx*dy-ry*dx)2/(dx2+dy2) ) 　この値が0〜線分の長さ√(dx2+dy2)に収まっていれば直線迄の距離が答え 
　負なら線分始点迄の距離　正なら線分終点迄の距離が答えとなります    

 場合分けを整理してpascalで書くと    
function Point2SegmentDistance(x0,y0,x1,y1,px,py:double):double;    
var rx,ry,dx,dy,qx,qy,L2,W:double;    
begin    
   dx:=x1-x0;    
   dy:=y1-y0;    
   L2:=dx*dx + dy*dy;    
   rx:=px-x0;    
   ry:=py-y0;    
   W :=rx*dx + ry*dy;    
   if L2 < W  then    
    begin    
     rx:=x1-px;    
     ry:=y1-py;    
     W :=rx*dx + ry*dy;    
    end;    
     if W<0  then result := sqrt(rx*rx+ry*ry)    
             else result := abs (rx*dy-ry*dx)/sqrt(L2);    
end;    
// 最善時の計算量に注目した場合    
function Point2SegmentDistance(x0,y0,x1,y1,px,py:double):double;    
var ly,cx:double;    
var pxm,pym,x0m,y0m:double;    
begin    
 pxm:=px-x1; pym:=py-y1;    
 x0m:=x0-x1; y0m:=y0-y1;    
    cx:= x0m*pxm+y0m*pym ; //線分がX軸に水平になるよう回転    
 if cx<0  then begin    
    Result:=hypot(pxm,pym);    
  end else begin    
   pxm:=px-x0; pym:=py-y0;    
   x0m:=x1-x0; y0m:=y1-y0;    
    cx:= x0m*pxm+y0m*pym ; //線分がX軸に水平になるよう回転    
 if cx<0  then begin    
   Result:=hypot(pxm,pym);    
  end else    
  Result:= abs(y0m*pxm-x0m*pym)/hypot(x0m,y0m) ;    
 end;    
end;

//ミッチェナー(Michener)のアルゴリズムを利用し
//円を水平線で塗り潰す。
//Michenerのアルゴリズムについては検索して下さい
//　水平線で塗る為の関数呼び出し場所を確認の事
//
procedure CirFill(Xc,Yc,r:Integer);
procedure xline(y,x0,x1:Integer);
begin
    Canvas.MoveTo(x0+Xc,y+Yc);//動作確認用の水平線描画
    Canvas.LineTo(x1+Xc,y+Yc);
end;
 var x,y,df:Integer;
begin
   x :=r;
   y :=0;
   df:=-2*r+3;  // 判定値の初期値 誤差は２倍して考える
 while  x >= y  do begin
   xline(-y,-x,x);
   xline( y,-x,x);
     if  df > 0  then begin  // (x-1,y+1) の半径誤差(x変化)の方が小さい
     xline( x,-y,y);
     xline(-x,-y,y);
       df := df+ 4*(-x+y) +10;  
       dec(x);
     end else begin  // (x,y+1) の半径誤差(x変化無) の方が小さい
       df := df + 4*y +6;
     end;
   inc(y);
 end;
end;

★楕円
媒介変数r: x=Asin(r),y=Bcos(r) これを回転/平行移動したもの

★最小２乗法による直線推定
　【回帰直線】
Ｙ＝ａＸ＋ｂのａ，ｂを求めるのに
ａ＝(XY -X*Y/N) /(XX-X*X/N)
ｂ＝(XX*Y-X*XY)/N/(XX-X*X/N)
XX=Σ(Xn*Xn) XY=Σ(Yn*Xn)
X =ΣXn Y =ΣYn
　【問題点】
通常知られている直線の推定方法は回帰直線による方法だが，
実際には誤差を最小にする方法でなく，Ｘ方向かＹの誤差の２
乗和なので，計測には向かない．回帰直線を求める方法では
回転変換したデータで計算して同じ結果が獲られない．

　【直線の方程式】
cY+sX+d=0 s*s+c*c=1;{または Y=st-d/c ;X=ct-d/s;}

とした時，点(x1,y1)から直線迄の距離z1が誤差になる．
z1=cy1+sx1+d

最小２乗法では誤差の２乗和を最小にするので

ΣZn*Zn を最小にするs,dを求める
ΣZn*Zn =Σ(cYn+sXn+d)^2

最小値は微分してゼロになる点なのでdで微分すると
2Σ(cYn+sXn+d)=0
cΣYn+sΣXn+Σd=0 ------------ a)

cで微分する為にsをcで微分する事を考える
s=(1-c^2)^(1/2) ---> s'=-c/cqrt(1-c*c)=-c/s
cで微分すると
Σ2(c*Yn+s*Xn+d)*( Yn - Xn*c/s)=0; 
Σ2(c*Yn+s*Xn+d)*(s*Yn - c*Xn)=0; 
Σ (c*Yn+s*Xn )*(s*Yn - c*Xn) +Σd*( s*Yn - c*Xn)=0;
Σ{c*Yn*s*Yn -c*Yn*c*Xn+s*Xn*s*Yn -s*Xn*c*Xn } =-Σd*( s*Yn - Xn*c)
Σ{s*c*Yn*Yn -c*c*Yn*Xn+s*s*Xn*Yn -s*c*Xn*Xn } =-Σd*( s*Yn - Xn*c)
Σ{s*c*(Yn*Yn-Xn*Xn) -(c*c-s*s)Yn*Xn } =-Σd*( s*Yn - Xn*c);
s*c*Σ(Yn*Yn-Xn*Xn) -(c*c-s*s)ΣYn*Xn =-d*Σ( s*Yn - Xn*c)

★３つの方程式をまとめると
YY=Σ(Yn*Yn) XX=Σ(Xn*Xn) XY=Σ(Yn*Xn)
X =ΣXn Y =ΣYn N=Σ1
c*c+s*s=1 ----------- 1)
c*Y+s*X+Nd=0 ------------ 2)
s*c(YY-XX)-(c*c-s*s)XY =-d(s*Y-c*X) 3)
3)に2)式を代入する
s*c(YY-XX)-(c*c-s*s)XY =( c*Y+s*X)(s*Y-c*X)/N
=( c*s*Y*Y-c*c*Y*X+s*s*X*Y-c*s*X*X)/N
=(-s*c*(X*X-Y*Y)-(c*c-s*s)*X*Y)/N
よって s*c(YY-XX+(X*X-Y*Y)/N) =(c*c-s*s)(XY-X*Y/N)
s*c/(c*c-s*s)= (XY-X*Y/N )/(YY-XX+(X*X-Y*Y)/N)
s*c/(c*c-s*s)は s=sin(a)c=cos(a)とすれば
sin(2a)/cos(2a)/2=tan(2a)/2なので

a=arctan{2*(XY-X*Y/N )/(YY-XX+(X*X-Y*Y)/N)}/2
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
と求められる
aからsin,cosを求めて
d=-(cY+sX)/N;
Σ(cYn+sXn+d)*(cYn+sXn+d)
=c*c*YY + s*s*XX + d*d*N + 2*c*s*XY + 2*d*(c*Y+s*X)
=c*c*YY + s*s*XX + 2*c*s*XY -d*d*N

★最小２乗法による円弧推定
　【円弧の方程式】
(X-a)^2+(Y-b)^2=r^2

とした時，点(x1,y1)から円弧迄の距離z1が誤差になる．２乗誤差は
Σz1=Σ{sqrt((X1-a)^2+(Y1-b)^2)-r}^2
=Σ{sqrt((X1-a)^2+(Y1-b)^2)^2 -2r*sqrt((X1-a)^2+(Y1-b)^2)+r*r
=Σ(X1-a)^2+(Y1-b)^2) -2r*sqrt((X1-a)^2+(Y1-b)^2)+r*r
=XX-2aX+aaN+YY-2bY+bbN+rrN-2rΣ{sqrt((X1-a)^2+(Y1-b)^2)}

右辺のsqrt項がある為 このような距離誤差の解を得るのは困難．
その代用として中心からの距離の２乗の誤差の２乗を最小にする．
いわば面積誤差を最小にする．
そうする事で，円の内側と外側とで誤差の重みが違う事になるが，
それは良い方向に働くようなので可とする．

Σ zn= Σ((Xn-a)^2+(Yn-b)^2-r^2)^2
a,b,rで微分すると
Σ-2((Xn-a)^2+(Yn-b)^2-r^2)2(Xn-a)=0 ---a)
Σ-2((Xn-a)^2+(Yn-b)^2-r^2)2(Yn-b)=0 ---b)
Σ-4((Xn-a)^2+(Yn-b)^2-r^2)r=0 ---r)
分解して
Σ(Xn*Xn+a*a-2aXn+Yn*Yn+b*b-2bYn-r^2)(Xn-a)=0 ---a1)
Σ(Xn*Xn+a*a-2aXn+Yn*Yn+b*b-2bYn-r^2)(Yn-b)=0 ---b1)
Σ(Xn*Xn+a*a-2aXn+Yn*Yn+b*b-2bYn-r^2)=0 ---r1)
a1),b1)の中で(Xn-a)の項の内，-aはr1)式より相手がゼロなので
Σ(Xn*Xn+a*a-2aXn+Yn*Yn+b*b-2bYn-r^2)Xn=0 ---a2)
Σ(Xn*Xn+a*a-2aXn+Yn*Yn+b*b-2bYn-r^2)Yn=0 ---b2)
XXX=Σ(Xn*Xn*Xn)のように表示する．
XXX+aaX-2aXX+YYX+bbX-2bYX=rrX ---a3)
XXY+aaY-2aXY+YYY+bbY-2bYY=rrY ---b3)
XX +aaN-2aX +YY +bbN-2bY =rrN ---r3)

XXX-2aXX+YXY-2bXY =(rr-aa-bb)X ---a4)
XXY-2aXY+YYY-2bYY =(rr-aa-bb)Y ---b4)
XX -2aX +YY -2bY =(rr-aa-bb)N ---r4)
r4)の右辺ををa4),b4)に代入
XXX+XYY -2bXY-2aXX =X(XX+YY -2aX -2bY)/N ---(a2
YYY+XXY -2bYY-2aXY =Y(XX+YY -2aX -2bY)/N ---(b2

XXX+XYY -2bXY-2aXX =X*(XX+YY)/N -2aX*X/N -2bY*X/N ---(a2
YYY+XXY -2bYY-2aXY =Y*(XX+YY)/N -2aX*Y/N -2bY*Y/N ---(b2

2a( X*X/N -XX)+ 2b(X*Y/N -XY) =X(XX+YY)/N-XXX-XYY ----1)
2a( X*Y/N -XY)+ 2b(Y*Y/N -YY) =Y(YY+XX)/N-YYY-XXY ----2)
r*r=(XX+YY-2.0(a*X+b*Y) )/N+a*a+b*b; ----3)

この３つの方程式を解けば良い．

  重心というのは、そこを持つと左右のバランスが取れる位置

  面積を細かく分割したものをdSとするとそれを集めた 
    ∫dS が面積
  座標Pとその微少な面積を積和し、それを全体の面積で割る
    ∫PdS/∫dS  これが重心
    xn座標に重みgnがかかっていれば (Σxn*gn)/(Σgn)という形式になります


  面積を求める方法は
  1) 「原点を頂点とする３角形」で分割して積算  と
  2) 「X軸/Y軸への台形」で分割して積算
  の２つが良く使われています。
  2)の方法のサンプル
  精度が要求される場合は、重心付近を原点に移動して計算します。
  この時の重心付近は、単に多角形の頂点の平均でもOK

  どちらも、線分の進行方向の右/左側に閉曲面があるかどうかで面積の符号
  を反転させれば良いのです。 
  重心だけが必要なら、X軸とY軸で符号を反対にして計算し 面積が負数になっ
  てもモーメントも負数になるから右回/左回りは無視して計算すればok

  多角形を構成する M個の線分 Pn,Pn+1 が ある場合、
    |     /pn+1      この線分とx軸との 面積を求めてゆきます。
    |    / |         dx=x(n+1)-x(n) としてその符号そのまま
    |   /  |         台形の面積の公式を使い積算してゆきます
    |   |pn|         つまり dxが負なら減算されます
    +------------>X
 この方法で計算する時は輪郭線の接続順(回転方向)が問題になります。
 方向が逆になれば負数になります。

 PASCAL=DELPHI
program test;
uses
 SysUtils;
type
  TPointF = record
    x,y: Extended;
  end;
var
 Sx,Sy:Extended;
 Mx,My:Extended;
 dx,dy:Extended;
 px,py:Extended;
 DT   :array [0..100]of TPointF;
 i,N:integer;
function PointF(x,y:Extended):TPointF;
 begin
  Result.x:=x;
  Result.y:=y;
 end; 
begin
 N:=0;
 DT[N]:=PointF(31,37);inc(N);
 DT[N]:=PointF(43,53);inc(N);
 DT[N]:=PointF(27,47);inc(N);
 DT[N]:=PointF(31,62);inc(N);
 DT[N]:=PointF(10,51);inc(N);
 DT[N]:=DT[0]; //終点は始点と同じ

 Sy:=0 ; Sx:=0 ;  //面積累積レジスタをクリア
 My:=0 ; Mx:=0 ;  //モーメント累積レジスタをクリア
 for i:=0 to N-1 do begin
  dy:=(DT[i+1].y-DT[i].y);
  dx:=(DT[i+1].x-DT[i].x);
   {      三角形                 +  四角形  }
  Sx:=Sx+dx*dy/2                 +dx*DT[i].y ;               //台形の面積
  Sy:=Sy+dx*dy/2                 +dy*DT[i].x ;               //台形の面積
  Mx:=Mx+dx*dy/2*(DT[i+1].x-dx/3)+dx*DT[i].y*(DT[i+1].x-dx/2);//右回モメント
  My:=My+dx*dy/2*(DT[i+1].y-dy/3)+dy*DT[i].x*(DT[i+1].y-dy/2);//右回モメント
 end;
//計算結果は Sx=-Syになりますから、ホントはSx,Syの一方を計算すればok
  px:=Mx/Sx; //重心X座標
  py:=My/Sy; //重心Y座標
end.

  いわゆる３次スプラインは行列演算が入るので結構厄介ですし、
  全部の点列が得られないと計算に入れません。
  比べてこの方式は、点列が入手出来る度に計算出来る利点があります。
  その上、混合スプラインや雲形定規スプラインに比べて接続感が良好。
  また、滑らか度合いを係数 k で調整可能です。
  さらに、整数化すれば、Win32APIにあるBezier系の関数を使って
  高速描画出来る特徴もあります。
   注)この方法は一応裏目小僧オリジナルの方法と信じていますが、
    　　　オリジナル性の調査はしていません。


DELPHI
type TDPoint = record
x: double;
y: double;
end;
type TDp4 = record
x1,y1: double;
x2,y2: double;
x3,y3: double;
x4,y4: double;
end;
///////////////////////////////////////////////////
// ４点を与えてp0〜p1の範囲のベジェ点化して返す
function BezierCvt(pr,p0,p1,p2:TDPoint):TDp4;
var dx,dy,k,r,b:double;
begin
 k:=0.3;  //kは 0.1〜0.3の範囲  小さい程折線に近い
 with Result do
  begin
    x1:=p0.x;
    y1:=p0.y;
    x4:=p1.x;
    y4:=p1.y;
    dx:=x4-x1;
    dy:=y4-y1;
    r:=k*sqrt(dx*dx+dy*dy);
    //  x2,y2 は p0点からpr-p1直線延長上距離rの点
    dx:=p1.x-pr.x;
    dy:=p1.y-pr.y;
    b:=sqrt(dx*dx+dy*dy);
    x2:=x1+r*dx/b;
    y2:=y1+r*dy/b;

    //  x2,y2 は p1点からp2-p0直線延長上距離rの点
    dx:=p0.x-p2.x;
    dy:=p0.y-p2.y;
    b:=sqrt(dx*dx+dy*dy);
    x3:=x4+r*dx/b;
    y3:=y4+r*dy/b;
  end;
end;
///////////////////////////////////////////////////
//ベジェ点データqを与えてその区間 t=0〜1 で点を返す
function BezierData(q:TDp4;t:double):TDpoint;
var   v:double;
begin
 v:=1-t;
 with Result do with q do begin
 x:= v*v*v*x1+3*v*v*t*x2+3*v*t*t*x3+t*t*t*x4;
 y:= v*v*v*y1+3*v*v*t*y2+3*v*t*t*y3+t*t*t*y4;
 end;
end;
///////////////////////////////////////////////////
//呼び出し方のサンプル
// 始点、終点では　始点終点を重ねるだけで、まあまあの結果が得られます。
// なお、 閉曲線なら始点終点を互いに重ねるだけの事です。

procedure TForm1.FormPaint(Sender: TObject);
var i,j:integer;
t:double;
p:TDpoint;
q:TDp4;
var dd:array [0..10] of TDpoint;
begin
 i:=0;
 dd[i].x:=100; dd[i].y:=10; inc(i); //始点を重ねる
 dd[i].x:=100; dd[i].y:=10; inc(i);
 dd[i].x:=150; dd[i].y:=200; inc(i);
 dd[i].x:=200; dd[i].y:=150; inc(i);
 dd[i].x:=250; dd[i].y:=300; inc(i);
 dd[i].x:=300; dd[i].y:=400; inc(i);
 dd[i].x:=350; dd[i].y:=300; inc(i);
 dd[i].x:=400; dd[i].y:=300; inc(i);
 dd[i].x:=400; dd[i].y:=300; inc(i);//終点も同じ点を重ねる
 with Canvas do begin
 Pen.color:=clRed;
 MoveTo(round(dd[1].x),round(dd[1].y));
 for i:=2 to 7 do LineTo(round(dd[i].x),round(dd[i].y));
 Pen.color:=clBlack;
 for i:=1 to 6 do
 begin
  q:= BezierCvt(dd[i-1],dd[i],dd[i+1],dd[i+2]) ;
  MoveTo(round(q.x1),round(q.y1));
  for j:=1 to 100 do 
   with BezierData(q,j/100) do LineTo(round(x),round(y));
   end;
 //なお このqの値から直接WIN32APIのベジェ系の関数を呼ぶのは簡単でしょう
 end;
end;

Ｂスプラインの解説

http://www.sra.co.jp/people/aoki/Kjo/Manuals/Kjo/Curve/Bspline/Bspline.htm

曲線・曲面入門ベジェ曲線 スプライン曲線 その他の曲線と比較

http://www.ccad.sccs.chukyo-u.ac.jp/manual/prgrm/curve/index.htm

http://ww3.tiki.ne.jp/~y-kawano/daigaku/cg-6.html

 混合系のスプライン
  x(t)= (1-t)3*x1+3(1-t)^2*t*x2+3(1-t)*t^2*x3+t3*x4
  y(t)= (1-t)3*y1+3(1-t)^2*t*y2+3(1-t)*t^2*y3+t3*y4

雲型定規式は 混合スプラインと同じで補間式だけが違い、
x(t) = ( -x0 * t   *(t-1)*(t-1)
         +x1 *(t-1)*(3*t*t-2*t-2)
         -x2 * t   *(3*t*t-4*t-1)
         +x3 * t   *  t *(t-1)    )/2