source: trunk/packages/vizservers/nanovis/PCASplit.cpp @ 1028

Last change on this file since 1028 was 1028, checked in by gah, 13 years ago

various cleanups

File size: 9.7 KB
Line 
1
2#include <stdio.h>
3#include "PCASplit.h"
4
5#define WANT_STREAM                  // include.h will get stream fns
6#define WANT_MATH                    // include.h will get math fns
7                                     // newmatap.h will get include.h
8#include <newmatap.h>                // need matrix applications
9#include <newmatio.h>                // need matrix output routines
10#include <newmatrc.h>
11
12#ifdef use_namespace
13using namespace NEWMAT;              // access NEWMAT namespace
14#endif
15
16namespace PCA {
17
18PCASplit::PCASplit() :
19    _maxLevel(4),
20    _minDistance(0.5f),
21    _distanceScale(0.2f),
22    _indexCount(0),
23    _finalMaxLevel(0)
24{
25    _indices = new unsigned int[MAX_INDEX];
26    _memClusterChunk1 = new ClusterListNode[MAX_INDEX];
27    _memClusterChunk2 = new ClusterListNode[MAX_INDEX];
28    _curMemClusterChunk = _memClusterChunk1;
29    _memClusterChunkIndex = 0;
30}
31
32PCASplit::~PCASplit()
33{
34    delete [] _indices;
35    delete [] _memClusterChunk1;
36    delete [] _memClusterChunk2;
37}
38
39void
40PCASplit::computeCentroid(Point* data, int count, Vector3& mean)
41{
42    float sumx= 0, sumy = 0, sumz = 0;
43    float size = 0;
44    float sumsize = 0;
45    for (int i = 0; i < count; ++i) {
46        size = data[i].size;
47        sumx += data[i].position.x * size;
48        sumy += data[i].position.y * size;
49        sumz += data[i].position.z * size;
50        sumsize += size;
51    }
52    sumsize = 1.0f / sumsize;
53    mean.set(sumx * sumsize, sumy * sumsize, sumz * sumsize);
54}
55
56void
57PCASplit::computeCovariant(Point* data, int count, const Vector3& mean,
58                           float* m)
59{
60    memset(m, 0, sizeof(float) * 9);
61
62    for (int i = 0; i < count; ++i) {
63        m[0] += (data[i].position.x - mean.x) * (data[i].position.x - mean.x);
64        m[1] += (data[i].position.x - mean.x) * (data[i].position.y - mean.y);
65        m[2] += (data[i].position.x - mean.x) * (data[i].position.z - mean.z);
66        m[4] += (data[i].position.y - mean.y) * (data[i].position.y - mean.y);
67        m[5] += (data[i].position.y - mean.y) * (data[i].position.z - mean.z);
68        m[8] += (data[i].position.z - mean.z) * (data[i].position.z - mean.z);
69    }
70   
71    float invCount = 1.0f / (count - 1);
72    m[0] *= invCount;
73    m[1] *= invCount;
74    m[2] *= invCount;
75    m[4] *= invCount;
76    m[5] *= invCount;
77    m[8] *= invCount;
78    m[3] = m[1];
79    m[6] = m[2];
80    m[7] = m[5];
81}
82
83void
84PCASplit::computeDistortion(Point* data, int count, const Vector3& mean,
85                            float& distortion, float& finalSize)
86{
87    distortion = 0.0f;
88    finalSize = 0.0f;
89       
90    float maxSize = 0.0f;
91    float distance;
92    for (int i = 0; i < count; ++i) {
93        // sum
94        distance = mean.distanceSquare(data[i].position.x, data[i].position.y, data[i].position.z);
95        distortion += distance;
96       
97        if (data[i].size > maxSize) {
98            maxSize = data[i].size;
99        }
100       
101        /*
102          finalSize += data[i].size * sqrt(distance);
103        */
104        if (distance > finalSize) {
105            finalSize = distance;
106        }
107    }
108    // equation 2
109    //finalSize = 0.5f * sqrt (finalSize) / (float) (count - 1) + maxSize;
110    finalSize = sqrt (finalSize) + maxSize;
111}
112       
113void
114PCASplit::init()
115{
116    _curClusterNode = 0;
117    _curClusterCount = 0;
118}
119
120Cluster*
121PCASplit::createClusterBlock(ClusterListNode* clusterList, int count, int level)
122{
123    static int cc = 0;
124    cc += count;
125    Cluster* clusterBlock = new Cluster[count];
126
127    _clusterHeader->numOfClusters[level - 1] = count;
128    _clusterHeader->startPointerCluster[level - 1] = clusterBlock;
129
130    printf("Cluster created %d [in level %d]:total %d\n", count, level, cc);
131    fflush(stdout);
132       
133    int i = 0;
134    ClusterListNode* clusterNode = clusterList;
135    while (clusterNode) {
136        clusterBlock[i].centroid = clusterList->data->centroid_t;
137        clusterBlock[i].level = level;
138        clusterBlock[i].scale = clusterList->data->scale_t;
139       
140        clusterNode = clusterNode->next;
141        ++i;
142    }
143    if (count != i) {
144        printf("ERROR\n");
145    }
146    return clusterBlock;
147}
148
149ClusterAccel* PCASplit::doIt(Point* data, int count)
150{
151    init();
152
153    _clusterHeader = new ClusterAccel(_maxLevel);
154       
155    Cluster* root = new Cluster;
156    Cluster_t* cluster_t = new Cluster_t();
157    cluster_t->points_t = data;
158    cluster_t->numOfPoints_t = count;
159    root->level = 1;
160       
161    _clusterHeader->root = root;
162    _clusterHeader->numOfClusters[0] = 1;
163    _clusterHeader->startPointerCluster[0] = root;
164       
165    Vector3 mean;
166    float distortion, scale;
167       
168    computeCentroid(cluster_t->points_t, cluster_t->numOfPoints_t, mean);
169    computeDistortion(cluster_t->points_t, cluster_t->numOfPoints_t, mean,
170                      distortion, scale);
171       
172    cluster_t->centroid_t = mean;
173    cluster_t->scale_t = scale;
174       
175    float mindistance = _minDistance;
176    int level = 2;
177
178    ClusterListNode* clustNodeList = &(_memClusterChunk2[0]);
179    clustNodeList->data = cluster_t;
180    clustNodeList->next = 0;
181       
182    _curRoot = root;
183    do {
184        analyze(clustNodeList, _curRoot, level, mindistance);
185
186        mindistance *= _distanceScale;
187        ++level;
188
189        // swap memory buffer & init
190        _curMemClusterChunk = (_curMemClusterChunk == _memClusterChunk1)?
191            _memClusterChunk2 : _memClusterChunk1;
192        _memClusterChunkIndex = 0;
193
194        clustNodeList = _curClusterNode;
195
196    } while (level <= _maxLevel);
197
198    return _clusterHeader;
199}
200
201void
202PCASplit::addLeafCluster(Cluster_t* cluster)
203{
204    ClusterListNode* clusterNode =
205        &_curMemClusterChunk[_memClusterChunkIndex++];
206    clusterNode->data = cluster;
207    clusterNode->next = _curClusterNode;
208    _curClusterNode = clusterNode;
209    ++_curClusterCount;
210}
211
212void
213PCASplit::split(Point* data, int count, float limit)
214{
215    Vector3 mean;
216    float m[9];
217
218    computeCentroid(data, count, mean);
219    float scale, distortion;
220    computeDistortion(data, count, mean, distortion, scale);
221       
222    //if (distortion < limit)
223    if (scale < limit) {
224        Cluster_t* cluster_t = new Cluster_t();
225        cluster_t->centroid_t = mean;
226        cluster_t->points_t = data;
227        cluster_t->numOfPoints_t = count;
228        cluster_t->scale_t = scale;
229        addLeafCluster(cluster_t);
230        return;
231    }
232   
233    computeCovariant(data, count, mean, m);
234   
235    SymmetricMatrix A(3);
236    for (int i = 1; i <= 3; ++i) {
237        for (int j = 1; j <= i; ++j) {
238            A(i, j) = m[(i - 1) * 3 + j - 1];
239        }
240    }
241       
242    Matrix U; DiagonalMatrix D;
243    eigenvalues(A,D,U);
244    Vector3 emax(U(1, 3), U(2, 3), U(3, 3));
245   
246    int left = 0, right = count - 1;
247   
248    Point p;
249    for (;left < right;) {
250        while (left < count && emax.dot(data[left].position - mean) >= 0.0f) {
251            ++left;
252        }
253        while (right >= 0 && emax.dot(data[right].position - mean) < 0.0f) {
254            --right;
255        }
256        if (left > right) {
257            break;
258        }
259       
260        p = data[left];
261        data[left] = data[right];
262        data[right] = p;
263        ++left, --right;
264       
265    }
266   
267    if (left == 0 || right == count - 1) {
268        printf("error\n");
269        exit(1);
270    } else {
271        split(data, left, limit);
272        split(data + left, count - left, limit);
273    }
274}
275
276void
277PCASplit::analyze(ClusterListNode* clusterNode, Cluster* parent, int level,
278                  float limit)
279{
280    if (level > _maxLevel) {
281        return;
282    }
283    if (level > _finalMaxLevel) {
284        _finalMaxLevel = level;
285    }
286       
287    init();
288
289    ClusterListNode* clNode = clusterNode;
290       
291    // initialize the indexCount of indices
292    _indexCount = 0;
293    while (clNode) {
294        if (clNode->data) {
295            split(clNode->data->points_t, clNode->data->numOfPoints_t, limit);
296            _indices[_indexCount++] = _curClusterCount;
297        }
298        clNode = clNode->next;
299    }
300       
301    //Vector3 mean;
302    //computeCentroid(cluster->points, cluster->numOfPoints, mean);
303
304    // the process values of split are in _curClusterNode and _curClusterCount
305    ClusterListNode* curClusterNode = _curClusterNode;
306    unsigned int curClusterNodeCount = _curClusterCount;
307
308    if (curClusterNodeCount) {
309       
310        // create and init centroid
311        Cluster* retClusterBlock =
312            createClusterBlock(curClusterNode, curClusterNodeCount, level);
313       
314        _curRoot = retClusterBlock;
315       
316        if (level == _maxLevel) {
317            ClusterListNode* node = curClusterNode;
318            if (_indexCount > 0) {
319                // for parent
320                Point* points = new Point[curClusterNodeCount];
321               
322                parent[0].setChildren(retClusterBlock, _indices[0]);
323                parent[0].setPoints(points, _indices[0]);
324               
325                for (unsigned int i = 0, in = 0; i < curClusterNodeCount; ++i) {
326                    if (i >= _indices[in]) {
327                        in++;
328                        parent[in].setChildren(retClusterBlock + _indices[in - 1], _indices[in] - _indices[in - 1]);
329                        parent[in].setPoints(points + _indices[in - 1], _indices[in] - _indices[in - 1]);
330                    }
331                                       
332                    retClusterBlock[i].scale = node->data->scale_t;
333                    retClusterBlock[i].centroid = node->data->centroid_t;
334                    retClusterBlock[i].points = node->data->points_t;
335                    retClusterBlock[i].numOfPoints = node->data->numOfPoints_t;
336                    retClusterBlock[i].color.set(float(rand()) / RAND_MAX,
337                                                 float(rand()) / RAND_MAX,
338                                                 float(rand()) / RAND_MAX,
339                                                 0.2);
340                   
341                    points[i].position = node->data->centroid_t;
342                    points[i].color.set(1.0f, 1.0f, 1.0f, 0.2f);
343                    points[i].size = node->data->scale_t;
344                   
345                    node = node->next;
346                }
347            }
348        } else {
349            Point* points = new Point[curClusterNodeCount];
350            ClusterListNode* node = curClusterNode;
351           
352            if (_indexCount > 0) {
353                parent[0].setPoints(points, _indices[0]);
354                parent[0].setChildren(retClusterBlock, _indices[0]);
355                for (int k = 1; k < _indexCount; ++k) {
356                    parent[k].setPoints(points + _indices[k - 1],
357                                        _indices[k] - _indices[k - 1]);
358                    parent[k].setChildren(retClusterBlock + _indices[k - 1],
359                                          _indices[k] - _indices[k - 1]);
360                }
361                               
362                // set points of sub-clusters
363                for (unsigned int i = 0; i < curClusterNodeCount; ++i) {
364                    points[i].position = node->data->centroid_t;
365                    points[i].color.set(1.0f, 1.0f, 1.0f, 0.2f);
366                    points[i].size = node->data->scale_t;
367                    node = node->next;
368                }
369               
370            }
371        }
372    }
373}
374
375} /*namespace PCA */
376
Note: See TracBrowser for help on using the repository browser.