最小二乗法による球の中心・半径のC言語による導出

質問

工学部の学生です。

最小二乗法による球の中心・半径のC言語による導出についてのご質問です。

ある物体の表面座標群を取得し、点iの位置（x座標、y座標、z座標）を
zahyou.x［i］
zahyou.y［i］
zahyou.z［i］　（iは0から300程度）
として保存している状況です。

この座標群にC言語で最小二乗法を適用し、中心座標と球半径を導出する場合、
どのようにすればよろしいでしょうか？

座標群が歪な物体であった場合、むりやりにでも導出することは可能でしょうか？
形だけでも点iの存在する空間の中心・半径っぽいものを求められると助かります。
（たとえば、最大のｘ座標と最小のｘ座標÷2≦導出半径に収まるなど、ありえない結果は除外できるでしょうか…）

実に他力本願な質問事項で心苦しいのですが、切羽詰っております。
恥を忍んで、どうか皆様方のご助力お願いします。

参考
http://questionbox.jp.msn.com/qa2652396.html

回答No.6

ur2c

> 切羽詰っております

なのに 3 日も応答がないのは「お前の考えはどうでも良い．C で最小 2 乗法をさっさと教えろ」ということかな，と Numerical Recipes in C の旧版を示しときます．

まあ，無駄だと思うわけですけど．これで，私にできることは，やり尽くしました．

参考URL：

http://www.nrbook.com/a/bookcpdf.php

回答No.5

ur2c

> 円の最小二乗法のサンプルソースがネット上に多くあったので、球も・・・と安易に考えました。

なら，最小2乗法の制約はない，と思って良いですね？ どうも以下のようなことらしい，と推測します．

O := 対象物体
R := robot
R から O を見て左手の座標系をとる．
Sensor は R に付いているものとする．
x := 中指，y := 親指，z := 人差し指
v := [v_x, v_y, v_z]

やりたいこと:
O を R のかごに入れる

そのためにできること:
R の進行方向と速度を制御する．

そのために必要な情報:
O の左端と右端の座標 x_l, x_r
R から O までの距離 z_O

そのために取れる情報:
V := {v}
card V = n ~= 300

O を外れた測定点は infty の値を持つとして
W := {v | v_z < infty} subset V

私なりの答です:

hat{x_l} := min {v_x | v in W}
hat{x_r} := max {v_x | v in W}
hat{z_O} := min {v_z | v in W}

ここで不安なのは hat{.} は n 個ある data points の 1 つの値だけに依存して決まっており，したがって不安定なこと．まずは上の実験をしてみて，実際に不安定なら，たとえば hat{x_l} は左端の数個の data を見て決めるようにすれば良い．その他も同様．

data update があるなら，指数平滑法くらいをやってみる．平滑 parameter は実験で決める．

かごを上下に調整できないなら y は不要．data に球面をあてはめるのも不要．中心座標も不要．目的と無関係な部分で不必要に話をややこしくしてます．手法を決める前に仕様を詰めるべきだという意見です．

回答No.4

FT56F001

>座標群が歪な物体であった場合、むりやりにでも導出することは可能でしょうか？
可能です。むりやり球形に近似する解が出ます。

>（たとえば、最大のｘ座標と最小のｘ座標÷2≦導出半径に収まるなど、ありえない結果は除外できるでしょうか…）
点の分布によっては，この不等式は崩れると思います。

基本的には#1さんの回答どおり，参考の式をCで書くだけの問題です。
4元連立方程式を解くのに，クラメールの公式かガウスの消去法かは悩みますが，簡単に書ける方で書けばよいでしょう。
Cの書き方なら「[技術者向] コンピューター」カテの方が適切でしょうが，
「参考の式の由来が分からない」なら，本カテで聞けばよいきっと説明が得られます。

なお，工学部のプログラミングあるいは数値計算の教育用の例題として，この問題は妥当かつ自然だと思います。

>既存のプログラム（たくさんある）を探して持って来ちゃうのが早いです。ネットで探せば、すぐみつかります。
実際の問題を解決する時，そうする方が早くてより高度なソフトが手に入る場合も多いです。

でも，中身も知らず，ブラックボックスとして使うだけなら，シロウトさんでもできます。
「一度は，簡単でよいから自分で書いてみる」という姿勢は，
学ぶ立場，次に来るべきものを創り出す立場としては重要なことです。それをすっとばすと，
砂上の楼閣の上になります。

回答No.3

alice_44

＞ Cで具体的に

カテゴリー違いですね。
具体的なコードが知りたいなら、
それなりの場所で訊きましょう。
ここは数学カテです。
最小二乗法をこの問題にあてはめた式は、
御自身のリンク先に解説済みのようです。

回答No.2

ur2c

> 工学部の学生です。

この問題は工学としては、へんな所が多いです。

> ある物体の表面座標群を取得し
> （iは0から300程度）

まず、測定方法の説明がありません。ですから、測定値をランダムサンプルと見なせるかどうかがわかりません。精度もわかりません。

> 中心座標と球半径

測定値に球面をあてはめる目的が書いてありません。そのため、どの程度の精度の出力が望まれているのか、わかりません。

> 最小二乗法による
> C言語による

なぜこのような制約が付いているのか、わかりません。

球面をあてはめるだけなら算数です。わざわざ「最小二乗法による」とするのはなぜでしょう？ 正規方程式を解かないといけないという意味でしょうか？ それとも計算結果が最小二乗法によるものと結果的に同等ならば良いのでしょうか？

どうしても最小二乗法によらねばならないなら、C より書きやすい言語がいくらもあります。なぜ C に限るのでしょう？ C でなければならないのは通常、マイコンとかに組込む必要性からです。けれど、それならプログラムをうんと小さくまとめる必要があるはずなので、ますます「最小二乗法による」といいう制約がかかっている理由が不明です。

先生が教育上の理由から「C言語で最小二乗法」という指定をしたのかもしれませんけど。

> この座標群にC言語で最小二乗法を適用し、中心座標と球半径を導出する場合、
どのようにすれば

どうしても「C言語で最小二乗法」なら、既存のプログラム（たくさんある）を探して持って来ちゃうのが早いです。ネットで探せば、すぐみつかります。ただし、英語で検索しないと出てこないかもしれません。

> 座標群が歪な物体であった場合、むりやりにでも導出することは可能でしょうか？

当然です。初めから球面だけがデータのモデルなのですから、立方体でも球とみなしちゃうわけで。

> 最大のｘ座標と最小のｘ座標÷2≦導出半径に収まるなど、ありえない結果は除外できるでしょうか…

歪みや測定誤差があるなら、これはありえる結果ですよね？

お礼

ご丁寧な回答ありがとうございます。

>> 工学部の学生です。
>この問題は工学としては、へんな所が多いです。
>> ある物体の表面座標群を取得し
>> （iは0から300程度）
>まず、測定方法の説明がありません。ですから、測定値をランダムサンプルと見なせるかどうかが
>わかりません。精度もわかりません。

測定方法はLRFを利用した3次元レーザスキャンです。
計測誤差はmmオーダとなっています。
また、測定対象は真球ではなく、お手玉のようなひしゃげた形状をしています。
さらに、一方向からのみの測定が条件となっているため、測定対象の半面のみ座標が明らかです。
測定対象の半径は30cm前後です。

>> 中心座標と球半径
>測定値に球面をあてはめる目的が書いてありません。そのため、どの程度の精度の出力が望ま
>れているのか、わかりません。

測定対象の中心に移動ロボットを向かわせます。
移動ロボットには測定対象捕獲用の籠がついていて大きさは全幅40cmです。
そのため中心座標の許容する誤差は5cmまでが望ましいです。

>> 最小二乗法による
>> C言語による
>なぜこのような制約が付いているのか、わかりません。

他の制御部がC言語によって構築され、そのソースへの書き加えが求められています。
また先生からの指導により、C以外の言語は残念ながら認められていません。

>球面をあてはめるだけなら算数です。わざわざ「最小二乗法による」とするのはなぜでしょう？ 正
>規方程式を解かないといけないという意味でしょうか？ それとも計算結果が最小二乗法によるも
>のと結果的に同等ならば良いのでしょうか？

球面をあてはめるだけ、という処理が不勉強ゆえ理解できません。
円の最小二乗法のサンプルソースがネット上に多くあったので、球も・・・と安易に考えました。
残念ながら球の最小二乗法や、球のカーブフィッティングのサンプルソースが見あたらなかったため、ご質問させて頂きました。

>> この座標群にC言語で最小二乗法を適用し、中心座標と球半径を導出する場合、
どのようにすれば
>どうしても「C言語で最小二乗法」なら、既存のプログラム（たくさんある）を探して持って来ちゃうのが早いです。ネットで探せば、すぐみつかります。ただし、英語で検索しないと出てこないかもしれません。

探したのですが・・・、すいません・・・

>> 最大のｘ座標と最小のｘ座標÷2≦導出半径に収まるなど、ありえない結果は除外できるでしょうか…
>歪みや測定誤差があるなら、これはありえる結果ですよね？

LRFによる計測座標の範囲外に中心が存在する、という結果がでないように、計測範囲外の解を除外し、中心座標が計測座標の系の中に内包するなかでの誤差ε＾２が最小になるような解を求めさせる。こういった処理はどのようにすればよろしいでしょうか。

実に厚かましいお願いなのですが、兎にも角にも、それっぽい中心が求まるサンプルソースなど御座いましたらお教え頂けると幸いです。

回答No.1

Tacosan

その「参考」にある式を C で書けばいいだけじゃないかなぁ.

お礼

おっしゃることごもっともです。
しかし残念ながらCで具体的にどうかくのか分かりません。

よろしければC言語を利用してどう行列式を解くのか、ソースを用いてご指南頂けないでしょうか？