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

[numeron]

　[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]

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 );　として変数に代入すれば１回で済む

[usao]

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

コード:

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

のようにそれ自体を１個の関数にまとめてしまえば良いかと．

[ym114]

シュミットの直行化はΣを使って一般化できます．
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]

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

[numeron]

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

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