グラム・シュミットの正規直交化法

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
numeron

グラム・シュミットの正規直交化法

#1

投稿記事 by numeron » 6年前

 [1.1] 自分が今行いたい事は何か
=>シュミットの正規直交化を行いたい
 [1.2] どのように取り組んだか
=>下記参照
 [1.3] どのようなエラーやトラブルで困っているか
=>プログラムは完成したがもっとsmartにしたい
 [1.4] 今何がわからないのか、知りたいか
=>力技が多い(printf etc...の多用など)ので
それらをfor文や関数で表記し
より分かりやすいプログラムにしたい

コード:

#include<stdio.h>
#include<float.h>
#include<math.h>
typedef struct vector{
  float x;
  float y;
  float z;
} VECTOR;

float naiseki(VECTOR *p,VECTOR *q);
float schmidt(VECTOR *p,VECTOR *q);
float normalize(VECTOR *p);
float main(void)
{
  VECTOR x1,x2,x3,a,b;
  float i,j,k;

  printf("\n");
  printf("3次ベクトルx1,x2,x3を入力してください.\n\n");
    
  printf("x1_x=");
  scanf("%f",&x1.x);
  printf("x1_y=");
  scanf("%f",&x1.y);
  printf("x1_z=");
  scanf("%f",&x1.z);
  printf("\n");
  
  printf("x2_x=");
  scanf("%f",&x2.x);
  printf("x2_y=");
  scanf("%f",&x2.y);
  printf("x2_z=");
  scanf("%f",&x2.z);
  printf("\n");
  
  printf("x3_x=");
  scanf("%f",&x3.x);
  printf("x3_y=");  
  scanf("%f",&x3.y);
  printf("x3_z="); 
  scanf("%f",&x3.z);
  printf("\n");

  printf("x1=(%f,%f,%f)\n",x1.x,x1.y,x1.z);
  printf("x2=(%f,%f,%f)\n",x2.x,x2.y,x2.z);
  printf("x3=(%f,%f,%f)\n",x3.x,x3.y,x3.z);
  printf("\n");
  
  printf("内積x1*x2=%f\n",naiseki(&x1,&x2));
  printf("内積x2*x3=%f\n",naiseki(&x2,&x3));
  printf("内積x3*x1=%f\n",naiseki(&x3,&x1));

  a.x=x2.x-schmidt(&x1,&x2)*x1.x;
  a.y=x2.y-schmidt(&x1,&x2)*x1.y;
  a.z=x2.z-schmidt(&x1,&x2)*x1.z;

  b.x=x3.x-(schmidt(&x1,&x3)*x1.x)-(schmidt(&a,&x3)*a.x);
  b.y=x3.y-(schmidt(&x1,&x3)*x1.y)-(schmidt(&a,&x3)*a.y);
  b.z=x3.z-(schmidt(&x1,&x3)*x1.z)-(schmidt(&a,&x3)*a.z);

  printf("\n");
  printf("3次ベクトルx1,x2,x3に対する直交基底v1,v2,v3は以下です.\n");
  printf("v1=(%f,%f,%f)\n",x1.x,x1.y,x1.z);
  printf("v2=(%f,%f,%f)\n",a.x,a.y,a.z);
  printf("v3=(%f,%f,%f)\n",b.x,b.y,b.z);
  printf("\n");

  i=normalize(&x1);
  j=normalize(&a);
  k=normalize(&b);

  x1.x=x1.x/i;
  x1.y=x1.y/i;
  x1.z=x1.z/i;
  a.x=a.x/j;
  a.y=a.y/j;
  a.z=a.z/j;
  b.x=b.x/k;
  b.y=b.y/k;
  b.z=b.z/k;

  printf("3次ベクトルx1,x2,x3に対する正規直交基底u1,u2,u3は以下です.\n");
  printf("u1=(%f,%f,%f)\n",x1.x,x1.y,x1.z);
  printf("u2=(%f,%f,%f)\n",a.x,a.y,a.z);
  printf("u3=(%f,%f,%f)\n",b.x,b.y,b.z);
  printf("\n");
  return(0);
}

float naiseki(VECTOR *p,VECTOR *q)
{
  float naiseki;
  naiseki= p->x*q->x+p->y*q->y+p->z*q->z;
  return(naiseki);
}

float schmidt (VECTOR *p,VECTOR *q)
{
  float schmidt;
  schmidt=((p->x*q->x)+(p->y*q->y)+(p->z*q->z))/((p->x*p->x)+(p->y*p->y)+(p->z*p->z));
  return(schmidt);
}

float normalize(VECTOR *p)
{
  float normalize;
  normalize=sqrt((p->x*p->x)+(p->y*p->y)+(p->z*p->z));
  return(normalize);
}

アバター
usao
記事: 1569
登録日時: 6年前

Re: グラム・シュミットの正規直交化法

#2

投稿記事 by usao » 6年前

VECTOR x1,x2,x3; → VECTOR Vec[3]; として,
for( int i=0; i<3; i++ ){ Vecの入力処理 } みたいな話でしょうか?

あとは,せっかく内積計算の関数があるので schmidt() や normalize() で使うとか.

気になる細かい点は
・0での除算チェックを入れた方がいい
・normalize()が正規化処理をしていないので名前はnorm()とかの方が合ってるかな? とか.

[追記]
・同じ引数で関数schmidt()を複数回呼んでいる → const float schmidt_result = schmidt( xxx, yyy ); として変数に代入すれば1回で済む
最後に編集したユーザー usao on 2013年7月09日(火) 17:41 [ 編集 1 回目 ]

アバター
usao
記事: 1569
登録日時: 6年前

Re: グラム・シュミットの正規直交化法

#3

投稿記事 by usao » 6年前

アルゴリズムを実装したものを今後使うのであれば,例えば,

コード:

//[グラムシュミットの正規直行化]
//nVecs本の入力ベクトルpSrcVecsに対する正規直行化処理結果をpDstResultVecsに入れる.
//(pSrcVecsとpDstResultVecsが同じだったり重複領域を持つような指定をしたらどうなるのか,とかも書いておくとよい)
//成功時は正の値を,失敗時は0を返す.
int GramSchmidt_Orthonormalization(
    const VECTOR *pSrcVecs,  //※要素数nVecs個以上の配列の先頭
    VECTOR *pDstResultVecs,  //※同上
    size_t nVecs  //ベクトル数
    )
{
    ....
}
のようにそれ自体を1個の関数にまとめてしまえば良いかと.

ym114

Re: グラム・シュミットの正規直交化法

#4

投稿記事 by ym114 » 6年前

シュミットの直行化はΣを使って一般化できます.
Wn = Vn - Σ_(i=0 to n-1){ <Wi, Vn>*Wi}
Wi...Viを正規化したベクトル
<a,b> ... べクトルとbの内積
* ... 乗算

◆1. 入力をsmartに
for文を使うためにも,配列は必須です.
VECTOR v[3];
また,scanfをまとめて一行にもできます.
3つの数字をスペースで挟んで入力してみてください.

コード:

for(int i=0; i<3; i++)
{
    printf("Enter V%d",i);
    scanf("%f%f%f", &v[i].x, &v[i].y, &v[i].z);
}
◆2. 計算をsmartに
そもそも,シュミットの直行化はベクトルを得るものではないでしょうか?
VECTOR W2 = Schmidt(v, 2);  ※vはベクトルの集合
このように呼び出せたらsmartだな思いました.
ポイントは,関数の中で自分の関数を呼び出している点です

コード:

/* Wiを求める関数*/
VECTOR Schmidt(VECTOR v[], int n)
{
    // コピーを作成
    VECTOR v_n = v[n];
    for(int i=0; i<=n-1; i++)
    {
        VECTOR w_i = Schmidt(v, i);         // Wiを求める(x,y,zのコピー)
        Multiply(&w_i, Naiseki(&v_n, &w_i));// <Wi,Vn>*Wiを計算
        Subtract(&v_n, &w_i);               // v_nから引く
    }
    return Normalize(&v_n);
}
void Multiply(VECTOR* v, float x); ..... 渡したベクトルをx倍する関数
void Subtract(VECTOR* a, VECTOR* b); ... 最初のベクトルから後のベクトルを引く関数
VECTOR Normalize(VECTOR* v) ... 渡したベクトルを正規化する関数

全容の一部ですが,これだけで済みます。

コード:

int main()
{
    VECTOR v[3];
    for(int i=0; i<3; i++)
    {
        printf("Enter V%d\n",i);
        scanf("%f%f%f", &v[i].x, &v[i].y, &v[i].z);
    }

    int n;
    printf("Enter Wn(n=?) to calc\n");
    scanf("%d", &n);
    VECTOR a = Schmidt(v, n);
   printf("W%d = %f,%f,%f", n, a.x, a.y, a.z);
}
Cの文法に慣れていないので,怒られたら変更してください.
あと何かあればまた誰かに聞いてみてください.

※変数命名はこの考え方にのっとりました
3次元,n=2のとき{V0 V1 V2}に対して
W0 = Normalize(W0)
W1 = V1 - <W0,V1>*W0
W2 = V2 - <W0,V2>*W0 - <W1,V2>*W1


(お断り)
自身を呼び出す関数を再帰関数と呼びますが,上の例ではW0を何度も計算するので少し無駄があります。
例えば,初めて計算したW_iを外部の配列に保存できれば,二回目以降の無駄な計算は減らせます.多少smartな書き方ではなくなりますが。
また,W2を計算したあとにW0やW1の値も配列から取り出すことができます.

アバター
usao
記事: 1569
登録日時: 6年前

Re: グラム・シュミットの正規直交化法

#5

投稿記事 by usao » 6年前

その書き方だと,せっかく計算した正規直交基底のn番目の基底ベクトルしかもらえないように見受けます.
通常,得たいのは正規直交基底の全ての基底ベクトルだと思うので
再帰処理の中で,各Wiの計算結果を結果格納域に保存するようにすれば(無駄な計算も自然と省けますし)良いように思います.

numeron

Re: グラム・シュミットの正規直交化法

#6

投稿記事 by numeron » 6年前

返信が遅くなってしまい申し訳ありません、質問者のnumeronです。

皆様のおかげで、より見やすいプログラムを作ることができました(ベクトルの配列化、関数norm,normalizeの定義など)。本当にありがとうございました。また質問することもあると思いますので、よろしくお願いします。

閉鎖

“C言語何でも質問掲示板” へ戻る