source: nanovis/branches/1.1/PCASplit.cpp @ 5722

Last change on this file since 5722 was 4889, checked in by ldelgass, 5 years ago

Merge r3611:3618 from trunk

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