← Back to branch summary

~hrg/hrg-packaging/zeo

« back to all changes in this revision

Viewing changes to src/zeo_vec3d_pca.c

  • Committer: Yosuke Matsusaka
  • Date: 2016-01-12 04:55:51 UTC
  • Revision ID: yosuke.matsusaka@gmail.com-20160112045551-9xbrsuokuwqdfnc9
initial import

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/zeo_vec3d_pca.c
 
1
/* Zeo - Z/Geometry and optics computation library.
 
2
 * Copyright (C) 2005 Tomomichi Sugihara (Zhidao)
 
3
 *
 
4
 * zeo_vec3d_pca - principal component analysis for 3D vectors.
 
5
 */
 
6
 
 
7
#include <zeo/zeo_mat3d.h>
 
8
 
 
9
/* zVec3DBarycenterPL
 
10
 * - barycenter of vector cloud given by a list.
 
11
 */
 
12
zVec3D *zVec3DBarycenterPL(zVec3DList *vl, zVec3D *c)
 
13
{
 
14
  zVec3DListCell *vc;
 
15
 
 
16
  zVec3DClear( c );
 
17
  zListForEach( vl, vc )
 
18
    zVec3DAddDRC( c, vc->data );
 
19
  return zVec3DDivDRC( c, zListNum(vl) );
 
20
}
 
21
 
 
22
/* zVec3DBarycenter
 
23
 * - barycenter of vector cloud given by an array.
 
24
 */
 
25
zVec3D *zVec3DBarycenter(zVec3D v[], int num, zVec3D *c)
 
26
{
 
27
  register int i;
 
28
 
 
29
  zVec3DClear( c );
 
30
  for( i=0; i<num; i++ )
 
31
    zVec3DAddDRC( c, &v[i] );
 
32
  return zVec3DDivDRC( c, num );
 
33
}
 
34
 
 
35
/* zVec3DPCA_PL
 
36
 * - PCA against vector cloud given by a list.
 
37
 */
 
38
zVec3D *zVec3DPCA_PL(zVec3DList *vl, zVec3D evec[])
 
39
{
 
40
  zMat3D pv, vm;
 
41
  double eval[3];
 
42
  zVec3DListCell *vc;
 
43
 
 
44
  zMat3DClear( &vm );
 
45
  zListForEach( vl, vc ){
 
46
    zMat3DDyad( vc->data, vc->data, &pv );
 
47
    zMat3DAddDRC( &vm, &pv );
 
48
  }
 
49
  zMat3DSymEig( &vm, eval, evec );
 
50
  return evec;
 
51
}
 
52
 
 
53
/* zVec3DPCA
 
54
 * - PCA against vector cloud given by an array.
 
55
 */
 
56
zVec3D *zVec3DPCA(zVec3D v[], int num, zVec3D evec[])
 
57
{
 
58
  zMat3D pv, vm;
 
59
  double eval[3];
 
60
  register int i;
 
61
 
 
62
  zMat3DClear( &vm );
 
63
  for( i=0; i<num; i++ ){
 
64
    zMat3DDyad( &v[i], &v[i], &pv );
 
65
    zMat3DAddDRC( &vm, &pv );
 
66
  }
 
67
  zMat3DSymEig( &vm, eval, evec );
 
68
  return evec;
 
69
}
 
70
 
 
71
/* zVec3DBaryPCA_PL
 
72
 * - barycenter of and PCA against vector cloud given by a list.
 
73
 */
 
74
zVec3D *zVec3DBaryPCA_PL(zVec3DList *vl, zVec3D *c, zVec3D evec[])
 
75
{
 
76
  zMat3D pv, vm;
 
77
  double eval[3];
 
78
  zVec3DListCell *vc;
 
79
  zVec3D dp;
 
80
 
 
81
  zVec3DBarycenterPL( vl, c );
 
82
  zMat3DClear( &vm );
 
83
  zListForEach( vl, vc ){
 
84
    zVec3DSub( vc->data, c, &dp );
 
85
    zMat3DDyad( &dp, &dp, &pv );
 
86
    zMat3DAddDRC( &vm, &pv );
 
87
  }
 
88
  zMat3DSymEig( &vm, eval, evec );
 
89
  return c;
 
90
}
 
91
 
 
92
/* zVec3DBaryPCA
 
93
 * - barycenter of and PCA against vector cloud given by an array.
 
94
 */
 
95
zVec3D *zVec3DBaryPCA(zVec3D v[], int num, zVec3D *c, zVec3D evec[])
 
96
{
 
97
  zMat3D pv, vm;
 
98
  double eval[3];
 
99
  zVec3D dp;
 
100
  register int i;
 
101
 
 
102
  zVec3DBarycenter( v, num, c );
 
103
  zMat3DClear( &vm );
 
104
  for( i=0; i<num; i++ ){
 
105
    zVec3DSub( &v[i], c, &dp );
 
106
    zMat3DDyad( &dp, &dp, &pv );
 
107
    zMat3DAddDRC( &vm, &pv );
 
108
  }
 
109
  zMat3DSymEig( &vm, eval, evec );
 
110
  return c;
 
111
}