source: projectionDesigner/trunk/projdesigner/include/gmtl/Fit/GaussPointsFit.h @ 375

Last change on this file since 375 was 4, checked in by Torben Dannhauer, 15 years ago
File size: 9.3 KB
Line 
1/************************************************************** ggt-head beg
2 *
3 * GGT: Generic Graphics Toolkit
4 *
5 * Original Authors:
6 *   Allen Bierbaum
7 *
8 * -----------------------------------------------------------------
9 * File:          GaussPointsFit.h,v
10 * Date modified: 2002/01/31 01:16:22
11 * Version:       1.3
12 * -----------------------------------------------------------------
13 *
14 *********************************************************** ggt-head end */
15/*************************************************************** ggt-cpr beg
16*
17* GGT: The Generic Graphics Toolkit
18* Copyright (C) 2001,2002 Allen Bierbaum
19*
20* This library is free software; you can redistribute it and/or
21* modify it under the terms of the GNU Lesser General Public
22* License as published by the Free Software Foundation; either
23* version 2.1 of the License, or (at your option) any later version.
24*
25* This library is distributed in the hope that it will be useful,
26* but WITHOUT ANY WARRANTY; without even the implied warranty of
27* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28* Lesser General Public License for more details.
29*
30* You should have received a copy of the GNU Lesser General Public
31* License along with this library; if not, write to the Free Software
32* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
33*
34 ************************************************************ ggt-cpr end */
35// Based on code from:
36// Magic Software, Inc.
37// http://www.magic-software.com
38//
39#ifndef _GMTL_GAUSSPOINTSFIT_H
40#define _GMTL_GAUSSPOINTSFIT_H
41
42// Fit points with a Gaussian distribution.  The center is the mean of the
43// points, the axes are the eigenvectors of the covariance matrix, and the
44// extents are the eigenvalues of the covariance matrix and are returned in
45// increasing order.  The last two functions allow selection of valid
46// vertices from a pool.  The return value is 'true' if and only if at least
47// one vertex was valid.
48
49#include <gmtl/Vec3.h>
50#include <gmtl/Point3.h>
51#include <gmtl/Numerics/Eigen.h>
52
53namespace gmtl
54{
55
56/*
57void MgcGaussPointsFit (int iQuantity, const MgcVector2* akPoint,
58    MgcVector2& rkCenter, MgcVector2 akAxis[2], MgcReal afExtent[2]);
59*/
60
61void GaussPointsFit (int iQuantity, const Point3* akPoint,
62    Point3& rkCenter, Vec3 akAxis[3], float afExtent[3]);
63
64/*
65bool MgcGaussPointsFit (int iQuantity, const MgcVector2* akPoint,
66    const bool* abValid, MgcVector2& rkCenter, MgcVector2 akAxis[2],
67    MgcReal afExtent[2]);
68*/
69
70bool GaussPointsFit (int iQuantity, const Vec3* akPoint,
71    const bool* abValid, Vec3& rkCenter, Vec3 akAxis[3],
72    float afExtent[3]);
73       
74
75// --- Implementations ---- //
76void GaussPointsFit (int iQuantity, const Point3* akPoint,
77    Point3& rkCenter, Vec3 akAxis[3], float afExtent[3])
78{
79    // compute mean of points
80    rkCenter = akPoint[0];
81    unsigned i;
82    for (i = 1; i < iQuantity; i++)
83        rkCenter += akPoint[i];
84    float fInvQuantity = 1.0f/iQuantity;
85    rkCenter *= fInvQuantity;
86
87    // compute covariances of points
88    float fSumXX = 0.0, fSumXY = 0.0, fSumXZ = 0.0;
89    float fSumYY = 0.0, fSumYZ = 0.0, fSumZZ = 0.0;
90    for (i = 0; i < iQuantity; i++)
91    {
92        Vec3 kDiff = akPoint[i] - rkCenter;
93        fSumXX += kDiff[Xelt]*kDiff[Xelt];
94        fSumXY += kDiff[Xelt]*kDiff[Yelt];
95        fSumXZ += kDiff[Xelt]*kDiff[Zelt];
96        fSumYY += kDiff[Yelt]*kDiff[Yelt];
97        fSumYZ += kDiff[Yelt]*kDiff[Zelt];
98        fSumZZ += kDiff[Zelt]*kDiff[Zelt];
99    }
100    fSumXX *= fInvQuantity;
101    fSumXY *= fInvQuantity;
102    fSumXZ *= fInvQuantity;
103    fSumYY *= fInvQuantity;
104    fSumYZ *= fInvQuantity;
105    fSumZZ *= fInvQuantity;
106
107    // compute eigenvectors for covariance matrix
108    gmtl::Eigen kES(3);
109    kES.Matrix(0,0) = fSumXX;
110    kES.Matrix(0,1) = fSumXY;
111    kES.Matrix(0,2) = fSumXZ;
112    kES.Matrix(1,0) = fSumXY;
113    kES.Matrix(1,1) = fSumYY;
114    kES.Matrix(1,2) = fSumYZ;
115    kES.Matrix(2,0) = fSumXZ;
116    kES.Matrix(2,1) = fSumYZ;
117    kES.Matrix(2,2) = fSumZZ;
118    kES.IncrSortEigenStuff3();
119
120    akAxis[0][Xelt] = kES.GetEigenvector(0,0);
121    akAxis[0][Yelt] = kES.GetEigenvector(1,0);
122    akAxis[0][Zelt] = kES.GetEigenvector(2,0);
123
124    akAxis[1][Xelt] = kES.GetEigenvector(0,1);
125    akAxis[1][Yelt] = kES.GetEigenvector(1,1);
126    akAxis[1][Zelt] = kES.GetEigenvector(2,1);
127
128    akAxis[2][Xelt] = kES.GetEigenvector(0,2);
129    akAxis[2][Yelt] = kES.GetEigenvector(1,2);
130    akAxis[2][Zelt] = kES.GetEigenvector(2,2);
131
132    afExtent[0] = kES.GetEigenvalue(0);
133    afExtent[1] = kES.GetEigenvalue(1);
134    afExtent[2] = kES.GetEigenvalue(2);
135}
136
137
138//
139bool GaussPointsFit (int iQuantity, const Vec3* akPoint,
140    const bool* abValid, Vec3& rkCenter, Vec3 akAxis[3],
141    float afExtent[3])
142{
143    // compute mean of points
144    rkCenter = ZeroVec3;
145    int i, iValidQuantity = 0;
146    for (i = 0; i < iQuantity; i++)
147    {
148        if ( abValid[i] )
149        {
150            rkCenter += akPoint[i];
151            iValidQuantity++;
152        }
153    }
154    if ( iValidQuantity == 0 )
155        return false;
156
157    float fInvQuantity = 1.0/iValidQuantity;
158    rkCenter *= fInvQuantity;
159
160    // compute covariances of points
161    float fSumXX = 0.0, fSumXY = 0.0, fSumXZ = 0.0;
162    float fSumYY = 0.0, fSumYZ = 0.0, fSumZZ = 0.0;
163    for (i = 0; i < iQuantity; i++)
164    {
165        if ( abValid[i] )
166        {
167            Vec3 kDiff = akPoint[i] - rkCenter;
168            fSumXX += kDiff[Xelt]*kDiff[Xelt];
169            fSumXY += kDiff[Xelt]*kDiff[Yelt];
170            fSumXZ += kDiff[Xelt]*kDiff[Zelt];
171            fSumYY += kDiff[Yelt]*kDiff[Yelt];
172            fSumYZ += kDiff[Yelt]*kDiff[Zelt];
173            fSumZZ += kDiff[Zelt]*kDiff[Zelt];
174        }
175    }
176    fSumXX *= fInvQuantity;
177    fSumXY *= fInvQuantity;
178    fSumXZ *= fInvQuantity;
179    fSumYY *= fInvQuantity;
180    fSumYZ *= fInvQuantity;
181    fSumZZ *= fInvQuantity;
182
183    // compute eigenvectors for covariance matrix
184    Eigen kES(3);
185    kES.Matrix(0,0) = fSumXX;
186    kES.Matrix(0,1) = fSumXY;
187    kES.Matrix(0,2) = fSumXZ;
188    kES.Matrix(1,0) = fSumXY;
189    kES.Matrix(1,1) = fSumYY;
190    kES.Matrix(1,2) = fSumYZ;
191    kES.Matrix(2,0) = fSumXZ;
192    kES.Matrix(2,1) = fSumYZ;
193    kES.Matrix(2,2) = fSumZZ;
194    kES.IncrSortEigenStuff3();
195
196    akAxis[0][Xelt] = kES.GetEigenvector(0,0);
197    akAxis[0][Yelt] = kES.GetEigenvector(1,0);
198    akAxis[0][Zelt] = kES.GetEigenvector(2,0);
199
200    akAxis[1][Xelt] = kES.GetEigenvector(0,1);
201    akAxis[1][Yelt] = kES.GetEigenvector(1,1);
202    akAxis[1][Zelt] = kES.GetEigenvector(2,1);
203
204    akAxis[2][Xelt] = kES.GetEigenvector(0,2);
205    akAxis[2][Yelt] = kES.GetEigenvector(1,2);
206    akAxis[2][Zelt] = kES.GetEigenvector(2,2);
207
208    afExtent[0] = kES.GetEigenvalue(0);
209    afExtent[1] = kES.GetEigenvalue(1);
210    afExtent[2] = kES.GetEigenvalue(2);
211
212    return true;
213}
214
215};
216
217/*
218void MgcGaussPointsFit (int iQuantity, const MgcVector2* akPoint,
219    MgcVector2& rkCenter, MgcVector2 akAxis[2], float afExtent[2])
220{
221    // compute mean of points
222    rkCenter = akPoint[0];
223    int i;
224    for (i = 1; i < iQuantity; i++)
225        rkCenter += akPoint[i];
226    float fInvQuantity = 1.0/iQuantity;
227    rkCenter *= fInvQuantity;
228
229    // compute covariances of points
230    float fSumXX = 0.0, fSumXY = 0.0, fSumYY = 0.0;
231    for (i = 0; i < iQuantity; i++)
232    {
233        MgcVector2 kDiff = akPoint[i] - rkCenter;
234        fSumXX += kDiff.x*kDiff.x;
235        fSumXY += kDiff.x*kDiff.y;
236        fSumYY += kDiff.y*kDiff.y;
237    }
238    fSumXX *= fInvQuantity;
239    fSumXY *= fInvQuantity;
240    fSumYY *= fInvQuantity;
241
242    // solve eigensystem of covariance matrix
243    MgcEigen kES(2);
244    kES.Matrix(0,0) = fSumXX;
245    kES.Matrix(0,1) = fSumXY;
246    kES.Matrix(1,0) = fSumXY;
247    kES.Matrix(1,1) = fSumYY;
248    kES.IncrSortEigenStuff2();
249
250    akAxis[0].x = kES.GetEigenvector(0,0);
251    akAxis[0].y = kES.GetEigenvector(1,0);
252    akAxis[1].x = kES.GetEigenvector(0,1);
253    akAxis[1].y = kES.GetEigenvector(1,1);
254
255    afExtent[0] = kES.GetEigenvalue(0);
256    afExtent[1] = kES.GetEigenvalue(1);
257}
258*/
259
260
261/*
262bool MgcGaussPointsFit (int iQuantity, const MgcVector2* akPoint,
263    const bool* abValid, MgcVector2& rkCenter, MgcVector2 akAxis[2],
264    float afExtent[2])
265{
266    // compute mean of points
267    rkCenter = MgcVector2::ZERO;
268    int i, iValidQuantity = 0;
269    for (i = 0; i < iQuantity; i++)
270    {
271        if ( abValid[i] )
272        {
273            rkCenter += akPoint[i];
274            iValidQuantity++;
275        }
276    }
277    if ( iValidQuantity == 0 )
278        return false;
279
280    float fInvQuantity = 1.0/iValidQuantity;
281    rkCenter *= fInvQuantity;
282
283    // compute covariances of points
284    float fSumXX = 0.0, fSumXY = 0.0, fSumYY = 0.0;
285    for (i = 0; i < iQuantity; i++)
286    {
287        if ( abValid[i] )
288        {
289            MgcVector2 kDiff = akPoint[i] - rkCenter;
290            fSumXX += kDiff.x*kDiff.x;
291            fSumXY += kDiff.x*kDiff.y;
292            fSumYY += kDiff.y*kDiff.y;
293        }
294    }
295    fSumXX *= fInvQuantity;
296    fSumXY *= fInvQuantity;
297    fSumYY *= fInvQuantity;
298
299    // solve eigensystem of covariance matrix
300    MgcEigen kES(2);
301    kES.Matrix(0,0) = fSumXX;
302    kES.Matrix(0,1) = fSumXY;
303    kES.Matrix(1,0) = fSumXY;
304    kES.Matrix(1,1) = fSumYY;
305    kES.IncrSortEigenStuff2();
306
307    akAxis[0].x = kES.GetEigenvector(0,0);
308    akAxis[0].y = kES.GetEigenvector(1,0);
309    akAxis[1].x = kES.GetEigenvector(0,1);
310    akAxis[1].y = kES.GetEigenvector(1,1);
311
312    afExtent[0] = kES.GetEigenvalue(0);
313    afExtent[1] = kES.GetEigenvalue(1);
314
315    return true;
316}
317*/
318
319#endif
Note: See TracBrowser for help on using the repository browser.