← Back to branch summary

~ubuntu-branches/ubuntu/vivid/meshlab/vivid

« back to all changes in this revision

Viewing changes to meshlab/src/meshlabplugins/filter_feature_alignment/feature_msc.h

  • Committer: Bazaar Package Importer
  • Author(s): Teemu Ikonen
  • Date: 2009-10-08 16:40:41 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20091008164041-0c2ealqv8b8uc20c
Tags: 1.2.2-1
http://bugs.debian.org/539544
* New upstream version
* Do not build filter_isoparametrization because liblevmar dependency
  is not (yet) in Debian
* Fix compilation with gcc-4.4, thanks to Jonathan Liu for the patch
  (closes: #539544)
* rules: Add compiler variables to the qmake call (for testing with new
  GCC versions)
* io_3ds.pro: Make LIBS and INCLUDEPATH point to Debian version of lib3ds
* io_epoch.pro: Make LIBS point to Debian version of libbz2
* control:
  - Move Homepage URL to the source package section
  - Update to standards-version 3.8.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
meshlab/src/meshlabplugins/filter_feature_alignment/feature_msc.h
 
1
#include <stdio.h>
 
2
#include <vcg/complex/trimesh/base.h>
 
3
#include <vcg/complex/trimesh/update/flag.h>
 
4
#include <vcg/complex/trimesh/update/curvature.h>
 
5
#include <vcg/complex/trimesh/update/quality.h>
 
6
#include <vcg/complex/trimesh/update/topology.h>
 
7
#include <vcg/complex/trimesh/clean.h>
 
8
#include <vcg/complex/trimesh/stat.h>
 
9
#include <vcg/math/histogram.h>
 
10
#include <vcg/complex/trimesh/smooth.h>
 
11
#include <vcg/complex/trimesh/clustering.h>
 
12
 
 
13
#include <meshlabplugins/filter_mls/apss.h>
 
14
#include <meshlabplugins/filter_mls/implicits.h>
 
15
 
 
16
using namespace vcg;
 
17
using namespace std;
 
18
using namespace GaelMls;
 
19
 
 
20
//class for a multiscale feature based on mean curvature: curvature is computed at differents levels of smoothness
 
21
template<class MESH_TYPE, int dim>
 
22
class SmoothCurvatureFeature
 
23
{
 
24
    public:        
 
25
 
 
26
    class Parameters
 
27
    {
 
28
        public:       
 
29
 
 
30
        enum CurvatureType {GAUSSIAN, MEAN, ABSOLUTE};
 
31
 
 
32
        struct Item
 
33
        {
 
34
            CurvatureType type;
 
35
            int scale;
 
36
            float lower_bound, upper_bound;
 
37
 
 
38
            Item(CurvatureType _type, int _scale, float _lower_bound = 0.0f, float _upper_bound = 1.0f):
 
39
                 type(_type),scale(_scale),lower_bound(_lower_bound),upper_bound(_upper_bound){}
 
40
        };
 
41
 
 
42
        vector<Item> featureDesc;
 
43
 
 
44
        bool add(CurvatureType cType, int smoothStep, float lower_bound = 0.0f, float upper_bound = 1.0f){
 
45
 
 
46
            assert(smoothStep>=0 & featureDesc.size()<dim & lower_bound>=0.0f & upper_bound<=1.0f & lower_bound<upper_bound);
 
47
            featureDesc.push_back(Item(cType, smoothStep, lower_bound, upper_bound));
 
48
            return true;
 
49
        }
 
50
    };
 
51
 
 
52
    typedef MESH_TYPE MeshType;
 
53
    typedef SmoothCurvatureFeature<MeshType,dim> FeatureType;
 
54
    typedef typename FeatureType::Parameters ParamType;
 
55
    typedef typename MeshType::ScalarType ScalarType;
 
56
    typedef typename MeshType::VertexType VertexType;
 
57
    typedef typename MeshType::VertexIterator VertexIterator;
 
58
    typedef typename MeshType::template PerVertexAttributeHandle<FeatureType> PVAttributeHandle;
 
59
 
 
60
    Point3<ScalarType> pos;
 
61
    Point3<ScalarType> normal;
 
62
    float description[dim];           
 
63
 
 
64
    SmoothCurvatureFeature();
 
65
    SmoothCurvatureFeature(VertexType& v);
 
66
    static char* getName();
 
67
    static int getRequirements();
 
68
    static float getNullValue();
 
69
    static bool isNullValue(float);
 
70
    static int getFeatureDimension();
 
71
    static void SetupParameters(ParamType& param );
 
72
    static bool HasBeenComputed(MeshType& m);
 
73
    static bool ComputeFeature(MeshType&, ParamType& param, CallBackPos *cb=NULL);
 
74
    static bool Subset(int, MeshType&, vector<FeatureType*>&, int, CallBackPos *cb=NULL);
 
75
    static bool CheckPersistency(FeatureType f);
 
76
 
 
77
    private:
 
78
    static MESH_TYPE* CreateSamplingMesh();
 
79
    static int SetupSamplingStructures(MeshType& m, typename MeshType::template PerVertexAttributeHandle<FeatureType>& fh, MeshType* samplingMesh, vector<FeatureType*>* vecFeatures);
 
80
    static pair<float,float> FindMinMax(vector<FeatureType*>& vec, int scale);
 
81
    static void SortByBinCount(vector<FeatureType*>& vecFeatures, float min, float max, int histSize, vector<FeatureType*>& sorted);
 
82
    static int SpatialSelection(vector<FeatureType*>& container, vector<FeatureType*>& vec, int k, float radius);
 
83
    static void PreCleaning(MeshType& m);
 
84
};
 
85
 
 
86
template<class MESH_TYPE, int dim>
 
87
inline SmoothCurvatureFeature<MESH_TYPE,dim>::SmoothCurvatureFeature(){}
 
88
 
 
89
template<class MESH_TYPE, int dim>
 
90
inline SmoothCurvatureFeature<MESH_TYPE,dim>::SmoothCurvatureFeature(VertexType& v):pos(v.P()),normal(v.N())
 
91
{
 
92
    normal.Normalize();
 
93
    for(int i=0; i<dim; i++) SmoothCurvatureFeature::getNullValue();
 
94
}
 
95
 
 
96
template<class MESH_TYPE, int dim> inline char* SmoothCurvatureFeature<MESH_TYPE,dim>::getName()
 
97
{
 
98
    return "SmoothCurvatureFeature";
 
99
}
 
100
 
 
101
/* Provides needed attribute to compute feature. A detailed list follows:
 
102
    MM_FACEFACETOPO required by RemoveNonManifoldVertex and RemoveNonManifoldFace
 
103
    MM_FACEFLAGBORDER required by RemoveNonManifoldFace
 
104
    MM_VERTCURV required by curvature computation
 
105
    MM_VERTQUALITY required by curvature and histogram computation
 
106
*/
 
107
template<class MESH_TYPE, int dim> inline int SmoothCurvatureFeature<MESH_TYPE,dim>::getRequirements()
 
108
{       
 
109
    return (MeshModel::MM_VERTCURV |
 
110
            MeshModel::MM_VERTQUALITY |
 
111
            MeshModel::MM_FACEFACETOPO |
 
112
            MeshModel::MM_FACEFLAGBORDER);
 
113
}
 
114
 
 
115
template<class MESH_TYPE, int dim> inline bool SmoothCurvatureFeature<MESH_TYPE,dim>::isNullValue(float val)
 
116
{
 
117
    return ( ((val==FeatureType::getNullValue()) | (math::IsNAN(val))) );
 
118
}
 
119
 
 
120
template<class MESH_TYPE, int dim> inline float SmoothCurvatureFeature<MESH_TYPE,dim>::getNullValue()
 
121
{
 
122
    return -std::numeric_limits<float>::max();
 
123
}
 
124
 
 
125
template<class MESH_TYPE, int dim> inline int SmoothCurvatureFeature<MESH_TYPE,dim>::getFeatureDimension()
 
126
{
 
127
    return dim;
 
128
}
 
129
 
 
130
//check persistence beetween scales: return true if description is valid for all scales, false otherwise
 
131
template<class MESH_TYPE, int dim> bool SmoothCurvatureFeature<MESH_TYPE,dim>::CheckPersistency(FeatureType f)
 
132
{
 
133
    for(int i = 0; i<FeatureType::getFeatureDimension(); i++)
 
134
    {
 
135
        if( FeatureType::isNullValue(f.description[i]) ) return false;
 
136
    }
 
137
    return true;
 
138
}
 
139
 
 
140
//parameters must be ordered according to smooth iterations
 
141
template<class MESH_TYPE, int dim> void SmoothCurvatureFeature<MESH_TYPE,dim>::SetupParameters(ParamType& param)
 
142
{
 
143
    param.add(Parameters::GAUSSIAN, 5, 0.4f, 0.9f);
 
144
    param.add(Parameters::MEAN, 5, 0.4f, 0.9f);
 
145
    param.add(Parameters::GAUSSIAN, 10, 0.4f, 0.9f);
 
146
    param.add(Parameters::MEAN, 10, 0.4f, 0.9f);
 
147
    param.add(Parameters::GAUSSIAN, 15, 0.4f, 0.9f);
 
148
    param.add(Parameters::MEAN, 15, 0.4f, 0.9f);
 
149
    assert(int(param.featureDesc.size())==getFeatureDimension());
 
150
}
 
151
 
 
152
template<class MESH_TYPE, int dim> void SmoothCurvatureFeature<MESH_TYPE,dim>::PreCleaning(MeshType& m)
 
153
{
 
154
    tri::Clean<MeshType>::RemoveZeroAreaFace(m);
 
155
    tri::Clean<MeshType>::RemoveDuplicateFace(m);
 
156
    tri::Clean<MeshType>::RemoveDuplicateVertex(m);
 
157
    tri::Clean<MeshType>::RemoveNonManifoldFace(m);
 
158
    //tri::Clean<MeshType>::RemoveNonManifoldVertex(m);
 
159
    tri::Clean<MeshType>::RemoveUnreferencedVertex(m);
 
160
    tri::Allocator<MeshType>::CompactVertexVector(m);
 
161
}
 
162
 
 
163
template<class MESH_TYPE, int dim> bool SmoothCurvatureFeature<MESH_TYPE,dim>::HasBeenComputed(MeshType &m)
 
164
{
 
165
    //checks if the attribute exist
 
166
    return tri::HasPerVertexAttribute(m,std::string(FeatureType::getName()));
 
167
}
 
168
 
 
169
template<class MESH_TYPE, int dim> bool SmoothCurvatureFeature<MESH_TYPE,dim>::ComputeFeature(MeshType &m, ParamType& param, CallBackPos *cb)
 
170
{
 
171
    //variables needed for progress bar callback
 
172
    float progBar = 0.0f;
 
173
    float offset = 100.0f/((FeatureType::getFeatureDimension()+2)*m.VertexNumber());
 
174
    if(cb) cb(0,"Computing features...");
 
175
 
 
176
    //allocates a custom per vertex attribute in which we can store pointers to features in the heap.
 
177
    PVAttributeHandle fh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(m, true);
 
178
    if(!tri::Allocator<MeshType>::IsValidHandle(m,fh)) return false;
 
179
 
 
180
    //clear the mesh to avoid wrong values during curvature computations
 
181
    PreCleaning(m);
 
182
 
 
183
    //copy vertex coords of the mesh before smoothing
 
184
    vector<Point3<ScalarType> > oldVertCoords(m.VertexNumber());
 
185
    for(unsigned int i = 0; i<m.vert.size(); ++i){
 
186
        oldVertCoords[i] = m.vert[i].P();
 
187
 
 
188
        if(cb){ progBar+=offset; cb((int)progBar,"Computing features..."); }
 
189
    }
 
190
 
 
191
    assert(int(oldVertCoords.size())==m.VertexNumber());
 
192
 
 
193
    //creates a feature for each vertex
 
194
    for(VertexIterator vi = m.vert.begin(); vi!=m.vert.end(); ++vi) fh[vi] = FeatureType(*vi);
 
195
 
 
196
    //loop trough scale levels
 
197
    int smooth_step = 0, smooth_accum = 0;
 
198
    for (unsigned int i = 0; i<FeatureType::getFeatureDimension(); i++, smooth_accum+=smooth_step)
 
199
    {
 
200
        smooth_step = param.featureDesc[i].scale - smooth_accum;
 
201
 
 
202
        tri::Smooth<MeshType>::VertexCoordLaplacian(m, smooth_step);//smooth mesh
 
203
        tri::UpdateCurvature<MeshType>::MeanAndGaussian(m);  //compute curvature
 
204
        //copy curvature values in the quality attributes; then build the histogram and take lower and upper bounds.
 
205
        switch(param.featureDesc[i].type){
 
206
            case Parameters::GAUSSIAN:{
 
207
                tri::UpdateQuality<MeshType>::VertexFromGaussianCurvature(m);
 
208
                break;
 
209
            }
 
210
            case Parameters::MEAN:{
 
211
                tri::UpdateQuality<MeshType>::VertexFromMeanCurvature(m);
 
212
                break;
 
213
            }
 
214
            case Parameters::ABSOLUTE:{
 
215
                tri::UpdateQuality<MeshType>::VertexFromAbsoluteCurvature(m);
 
216
                break;
 
217
            }
 
218
            default: assert(0);
 
219
        }
 
220
        Histogram<ScalarType> hist = Histogram<ScalarType>();
 
221
        tri::Stat<MeshType>::ComputePerVertexQualityHistogram(m, hist);
 
222
        float vmin = hist.Percentile(param.featureDesc[i].lower_bound); float vmax = hist.Percentile(param.featureDesc[i].upper_bound);
 
223
 
 
224
        //If curvature is beetween bounds and vertex is not a boundary, curvature is stored in the feature, otherwise the feature is set to an empty value.
 
225
        for(VertexIterator vi = m.vert.begin(); vi!=m.vert.end(); ++vi)
 
226
        {
 
227
            if( (!(*vi).IsB()) & ((*vi).Q()<=vmin) | ((*vi).Q()>=vmax) & (!FeatureType::isNullValue((*vi).Q())) )
 
228
                fh[vi].description[i] = (*vi).Q();
 
229
            else
 
230
                fh[vi].description[i] = FeatureType::getNullValue();
 
231
 
 
232
            if(cb){ progBar+=offset; cb((int)progBar,"Computing features..."); }
 
233
        }
 
234
    }
 
235
    //restore old coords
 
236
    for(unsigned int i = 0; i<m.vert.size(); ++i){
 
237
        m.vert[i].P() = oldVertCoords[i];
 
238
 
 
239
        if(cb){ progBar+=offset; cb((int)progBar,"Computing features..."); }
 
240
    }
 
241
 
 
242
    return true;
 
243
}
 
244
 
 
245
template<class MESH_TYPE, int dim>
 
246
MESH_TYPE* SmoothCurvatureFeature<MESH_TYPE,dim>::CreateSamplingMesh()
 
247
{
 
248
    MeshType* m = new MeshType();
 
249
    PVAttributeHandle fh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(*m, true);
 
250
    if(!tri::Allocator<MeshType>::IsValidHandle(*m,fh)){
 
251
        if(m) delete m;
 
252
        return NULL;
 
253
    }
 
254
    return m;
 
255
}
 
256
 
 
257
template<class MESH_TYPE, int dim>
 
258
int SmoothCurvatureFeature<MESH_TYPE,dim>::SetupSamplingStructures(MeshType& m, typename MeshType::template PerVertexAttributeHandle<FeatureType>& fh, MeshType* samplingMesh, vector<FeatureType*>* vecFeatures)
 
259
{
 
260
    int countFeatures = 0;
 
261
    PVAttributeHandle pmfh;
 
262
    if(samplingMesh){
 
263
        pmfh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(*samplingMesh);
 
264
        if(!tri::Allocator<MeshType>::IsValidHandle(*samplingMesh,pmfh)) return 0;
 
265
    }
 
266
 
 
267
    //fill the vector with all persistent features.
 
268
    for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi){
 
269
        //check persistence beetween scales: if feature is persistent, add a pointer in vecFeatures
 
270
        if( FeatureType::CheckPersistency(fh[vi])){
 
271
            countFeatures++;  //increment counter of valid features
 
272
            if(vecFeatures) vecFeatures->push_back(&(fh[vi]));
 
273
            if(samplingMesh){
 
274
                tri::Allocator<MeshType>::AddVertices(*samplingMesh, 1);
 
275
                samplingMesh->vert.back().ImportLocal(*vi);
 
276
                pmfh[samplingMesh->vert.back()] = fh[vi];
 
277
            }
 
278
        }
 
279
    }
 
280
 
 
281
    return countFeatures;
 
282
}
 
283
 
 
284
template<class MESH_TYPE, int dim>
 
285
bool SmoothCurvatureFeature<MESH_TYPE,dim>::Subset(int k, MeshType &m, vector<FeatureType*> &vecSubset, int sampType, CallBackPos *cb)
 
286
{
 
287
    //variables needed for progress bar callback
 
288
    float progBar = 0.0f;
 
289
    float offset = 100.0f/(m.VertexNumber() + k);
 
290
 
 
291
    //if attribute doesn't exist, return; else we can get a handle to the attribute
 
292
    PVAttributeHandle fh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(m);
 
293
    if(!tri::Allocator<MeshType>::IsValidHandle(m,fh)) return false;
 
294
 
 
295
    //create a vector to hold valid features that later will be sampled
 
296
    vector<FeatureType*>* vecFeatures = NULL;
 
297
    MeshType* poissonMesh = NULL;
 
298
    if(sampType==0) vecFeatures = new vector<FeatureType*>();
 
299
    else poissonMesh = CreateSamplingMesh();
 
300
 
 
301
    //fill the vector with all persistent features.
 
302
    int countFeatures = SetupSamplingStructures(m, fh, poissonMesh, vecFeatures);
 
303
    if(countFeatures<k) k = countFeatures;  //we can't extract more of what we got!
 
304
 
 
305
    //perform different kinds of sampling
 
306
    FeatureType** sampler = NULL;
 
307
    switch(sampType){
 
308
        case 0:{ //uniform sampling: uses vecFeatures
 
309
            sampler = FeatureAlignment<MeshType, FeatureType>::FeatureUniform(*vecFeatures, &k);
 
310
            break;
 
311
        }
 
312
        case 1:{ //poisson disk sampling: uses poissonMesh
 
313
            sampler = FeatureAlignment<MeshType, FeatureType>::FeaturePoisson(*poissonMesh, &k);
 
314
            break;
 
315
        }
 
316
        default: assert(0);
 
317
    }
 
318
 
 
319
    //store features into the returned vector
 
320
    for(int i=0; i<k; ++i){
 
321
        vecSubset.push_back(sampler[i]);
 
322
        if(cb){ progBar+=offset; cb(int(progBar),"Extracting features..."); }
 
323
    }
 
324
 
 
325
    if(vecFeatures) delete vecFeatures;   //clear useless data
 
326
    if(poissonMesh) delete poissonMesh;
 
327
    if(sampler) delete[] sampler;
 
328
 
 
329
    return true;
 
330
 
 
331
    ////UNCOMMENT FOLLOW TO GET A RARE FEATURE SELECTION////
 
332
/*
 
333
    int histSize = 10000;  //number of bins of the histogram
 
334
    vector<FeatureType*>* sorted = new vector<FeatureType*>();  //vector that holds features sorted by bin cont
 
335
 
 
336
    //compute min val e max val between all features; min and max are needed to bound the histogram
 
337
    pair<float,float> minmax = FindMinMax(*vecFeatures, 0);
 
338
 
 
339
    //fill multimap with features sorted by bin count
 
340
    SortByBinCount(*vecFeatures, minmax.first, minmax.second, histSize, *sorted);
 
341
 
 
342
    //select the first k entries from mulptimap, and put them into vecSubset
 
343
    typename vector<FeatureType*>::iterator it = sorted->begin();
 
344
    ///UNCOMMENT THIS LOOP FOR A SORTED SELECTION
 
345
    for(int i=0; i<k & it!=sorted->end(); i++, it++)
 
346
    {
 
347
        (*it)->selected = true;        //mark features as selected
 
348
        vecSubset.push_back(*it);
 
349
    }
 
350
    ///UNCOMMENT THIS LOOP FOR A EQUALIZED SELECTION
 
351
    int off = int(sorted->size()/k);
 
352
    for(it = sorted->begin(); it<sorted->end(); it=it+off)
 
353
    {
 
354
        (*it)->selected = true;        //mark features as selected
 
355
        vecSubset.push_back(*it);
 
356
    }
 
357
 
 
358
    //delete vector and set pointer to NULL for safety
 
359
    if(vecFeatures) delete vecFeatures;   //clear useless data
 
360
    sorted->clear();
 
361
    delete sorted;
 
362
    sorted = NULL;
 
363
    return true;
 
364
*/
 
365
}
 
366
 
 
367
//scan the vector of features and return a pair containig the min and max description values
 
368
template<class MESH_TYPE, int dim> pair<float,float> SmoothCurvatureFeature<MESH_TYPE,dim>::FindMinMax(vector<FeatureType*>& vec, int scale)
 
369
{
 
370
    assert(scale>=0 && scale<FeatureType::GetFeatureDimension());
 
371
 
 
372
    typename vector<FeatureType*>::iterator vi;
 
373
    pair<float,float> minmax = make_pair(numeric_limits<float>::max(),-numeric_limits<float>::max());
 
374
    for(vi = vec.begin(); vi!=vec.end(); ++vi)
 
375
    {
 
376
        if( !FeatureType::isNullValue((*vi)->description[scale]))  //invalid values are discarded
 
377
        {
 
378
            if( (*vi)->description[scale] < minmax.first) minmax.first=(*vi)->description[scale];
 
379
            if( (*vi)->description[scale] > minmax.second ) minmax.second=(*vi)->description[scale];
 
380
        }
 
381
    }
 
382
    return minmax;
 
383
}
 
384
 
 
385
//fill a vector that holds features pointers sorted by bin count; i.e first all the very infrequent features and last all the very frequent ones.
 
386
template<class MESH_TYPE, int dim> void SmoothCurvatureFeature<MESH_TYPE,dim>::SortByBinCount(vector<FeatureType*>& vecFeatures, float min, float max, int histSize, vector<FeatureType*>& sorted)
 
387
{
 
388
    //vector to hold bins ranges
 
389
    vector<float>* bins = new vector<float>(histSize, 0);
 
390
    //vector to count content of each bin
 
391
    vector<int> *hToC = new vector<int>(histSize, 0);
 
392
    //multimap to keep track of features for each bin
 
393
    multimap<int, FeatureType* >* hToF = new multimap<int, FeatureType* >();
 
394
    //multimap to store pairs of (count, FeatureType*). multimap is naturally sorted by count.
 
395
    multimap<int, FeatureType* >* cToF = new multimap<int, FeatureType* >();
 
396
 
 
397
    //offset beetween min and max is divided into histSize uniform bins
 
398
    for(unsigned int i = 0; i<bins->size(); i++) bins->at(i) = min + i*(max-min)/float(histSize);
 
399
 
 
400
    //build histogram, hToF and hToC; all with just one features scan
 
401
    typename multimap<int, FeatureType* >::iterator hfIt = hToF->begin();
 
402
    typename vector<FeatureType*>::iterator vi;
 
403
    for(vi = vecFeatures.begin(); vi!=vecFeatures.end(); ++vi)
 
404
    {
 
405
        float val = (*vi)->description[0];
 
406
        // lower_bound returns an iterator pointing to the first element "not less than" val, or end() if every element is less than val.
 
407
        typename vector<float>::iterator it = lower_bound(bins->begin(),bins->end(),val);
 
408
        //if val is out of range, skip iteration and take in account next val
 
409
        if(it==bins->begin() || it==bins->end()) continue;
 
410
        //determine in which bin got to insert val
 
411
        int binId = (it - bins->begin())-1;
 
412
        assert(binId>=0);
 
413
        assert (bins->at(binId) < val);
 
414
        assert (val <= bins->at(binId+1));
 
415
 
 
416
        //increment bin count
 
417
        hToC->at(binId)++;
 
418
        //remember in which bin has been inserted this feature
 
419
        hfIt = hToF->insert(hfIt,make_pair(binId, (*vi)));
 
420
    }
 
421
 
 
422
    //delete bins and set pointer to NULL for safety
 
423
    bins->clear();
 
424
    delete bins;
 
425
    bins = NULL;
 
426
 
 
427
    //for all bin index insert in the multimap an entry with key bin count and value a feature. Entries are naturally
 
428
    //sorted by key in the multimap
 
429
    typename multimap<int, FeatureType* >::iterator cfIt = cToF->begin();
 
430
    pair< typename multimap<int, FeatureType* >::iterator, typename multimap<int, FeatureType* >::iterator> range;
 
431
    for(unsigned int i=0; i<hToC->size(); i++)
 
432
    {
 
433
        //if bin count is zero, then skip; nothing to put in the multimap
 
434
        if(hToC->at(i)!=0)
 
435
        {
 
436
            range = hToF->equal_range(i);
 
437
            assert(range.first!=range.second);
 
438
            for (; range.first!=range.second; ++range.first)
 
439
            {
 
440
                cfIt = cToF->insert( cfIt, make_pair(hToC->at(i), (*range.first).second) );
 
441
            }
 
442
        }
 
443
    }
 
444
 
 
445
    typename multimap<int, FeatureType* >::iterator it;
 
446
    for(it = cToF->begin(); it != cToF->end(); it++)
 
447
        sorted.push_back((*it).second);
 
448
 
 
449
    assert(sorted.size()==cToF->size());
 
450
 
 
451
    //delete all ausiliary structures and set pointers to NULL for safety
 
452
    hToF->clear();
 
453
    delete hToF;
 
454
    hToF = NULL;
 
455
    hToC->clear();
 
456
    delete hToC;
 
457
    hToC = NULL;
 
458
    cToF->clear();
 
459
    delete cToF;
 
460
    cToF = NULL;
 
461
}
 
462
 
 
463
template<class MESH_TYPE, int dim> int SmoothCurvatureFeature<MESH_TYPE,dim>::SpatialSelection(vector<FeatureType*>& container, vector<FeatureType*>& vec, int k, float radius)
 
464
{
 
465
    //variables to manage the kDTree which works on features position
 
466
    ANNpointArray   dataPts = NULL;     // data points
 
467
    ANNpoint        queryPnt = NULL;    // query points
 
468
    ANNkd_tree*     kdTree = NULL;      // search structure
 
469
 
 
470
    queryPnt = annAllocPt(3);
 
471
 
 
472
    typename vector<FeatureType*>::iterator ci = container.begin();
 
473
    while(ci != container.end() && vec.size()<k)
 
474
    {
 
475
        if(vec.size()==0)
 
476
        {
 
477
            vec.push_back(*ci);
 
478
            (*ci)->selected = true;
 
479
            ci++;
 
480
            continue;
 
481
        }
 
482
 
 
483
        if(dataPts){ annDeallocPts(dataPts); dataPts = NULL; }
 
484
        if(kdTree){ delete kdTree; kdTree = NULL; }
 
485
        dataPts = annAllocPts(vec.size(), 3);
 
486
 
 
487
        for(int i = 0; i < vec.size(); i++)
 
488
        {
 
489
            for (int j = 0; j < 3; j++)
 
490
            {
 
491
                (dataPts[i])[j] = (ANNcoord)(vec.at(i)->pos[j]);
 
492
            }
 
493
        }
 
494
        //build search structure
 
495
        kdTree = new ANNkd_tree(dataPts,vec.size(),3);
 
496
 
 
497
        for (int j = 0; j < 3; j++)
 
498
        {
 
499
            queryPnt[j] = (ANNcoord)((*ci)->pos[j]);
 
500
        }
 
501
 
 
502
        //if there aren't features yet selected in the range distance
 
503
        if(!kdTree->annkFRSearch(queryPnt, math::Sqr(radius), 0, NULL, NULL, 0.0))
 
504
        {
 
505
             vec.push_back(*ci);
 
506
             (*ci)->selected = true;
 
507
        }
 
508
        ci++;
 
509
    }
 
510
 
 
511
    if(dataPts){ annDeallocPts(dataPts); dataPts = NULL; }
 
512
    if(queryPnt){ annDeallocPt(queryPnt); queryPnt = NULL; }
 
513
    if(kdTree){ delete kdTree; kdTree = NULL; }
 
514
 
 
515
    return vec.size();
 
516
}
 
517
 
 
518
//class for a multiscale feature based on APSS curvature. Works with point clouds.
 
519
template<class MESH_TYPE, int dim>
 
520
class APSSCurvatureFeature
 
521
{
 
522
    public:
 
523
 
 
524
    class Parameters
 
525
    {
 
526
        public:
 
527
 
 
528
        enum CurvatureType { MEAN, GAUSSIAN, K1, K2, APSS};
 
529
 
 
530
        struct Item
 
531
        {
 
532
            public:
 
533
            CurvatureType type;
 
534
            float lower_bound, upper_bound, scale;  //indicates how much sub sample the mesh. 1 = whole mesh, 0.5 = half mesh, etc
 
535
 
 
536
            Item(CurvatureType _type, float _scale, float _lower_bound = 0.0f, float _upper_bound = 1.0f):
 
537
                 type(_type),scale(_scale),lower_bound(_lower_bound),upper_bound(_upper_bound){}
 
538
        };
 
539
 
 
540
        vector<Item> featureDesc;
 
541
 
 
542
        bool add(CurvatureType cType, float subSampAmount, float lower_bound = 0.0f, float upper_bound = 1.0f)
 
543
        {
 
544
            assert(subSampAmount>=0 & subSampAmount<=1 & featureDesc.size()<dim & lower_bound>=0.0f & upper_bound<=1.0f & lower_bound<upper_bound);
 
545
            featureDesc.push_back(Item(cType, subSampAmount, lower_bound, upper_bound));
 
546
            return true;
 
547
        }
 
548
    };
 
549
 
 
550
    typedef MESH_TYPE MeshType;
 
551
    typedef APSSCurvatureFeature<MeshType,dim> FeatureType;
 
552
    typedef typename FeatureType::Parameters ParamType;
 
553
    typedef typename MeshType::ScalarType ScalarType;
 
554
    typedef typename MeshType::VertexType VertexType;
 
555
    typedef typename MeshType::VertexIterator VertexIterator;
 
556
    typedef typename MeshType::template PerVertexAttributeHandle<FeatureType> PVAttributeHandle;
 
557
 
 
558
    Point3<ScalarType> pos;
 
559
    Point3<ScalarType> normal;
 
560
    float description[dim];
 
561
 
 
562
    APSSCurvatureFeature();
 
563
    APSSCurvatureFeature(VertexType& v);
 
564
    static char* getName();
 
565
    static int getRequirements();
 
566
    static float getNullValue();
 
567
    static bool isNullValue(float);
 
568
    static int getFeatureDimension();
 
569
    static void SetupParameters(ParamType& param );
 
570
    static bool HasBeenComputed(MeshType& m);
 
571
    static bool ComputeFeature(MeshType&, ParamType& param, CallBackPos *cb=NULL);
 
572
    static bool Subset(int, MeshType&, vector<FeatureType*>&, int, CallBackPos *cb=NULL);
 
573
    static bool CheckPersistency(FeatureType f);
 
574
 
 
575
    private:
 
576
    static bool ComputeAPSS(MeshType& m, int type, float filterScale, int maxProjIter, float projAcc, float sphPar, bool accNorm, bool selectionOnly);
 
577
    static MESH_TYPE* CreateSamplingMesh();
 
578
    static int SetupSamplingStructures(MeshType& m, typename MeshType::template PerVertexAttributeHandle<FeatureType>& fh, MeshType* samplingMesh, vector<FeatureType*>* vecFeatures);
 
579
    static pair<float,float> FindMinMax(vector<FeatureType*>& vec, int scale);
 
580
    static void SortByBinCount(vector<FeatureType*>& vecFeatures, float min, float max, int histSize, vector<FeatureType*>& sorted);
 
581
    static int SpatialSelection(vector<FeatureType*>& container, vector<FeatureType*>& vec, int k, float radius);
 
582
    static void PreCleaning(MeshType& m);
 
583
};
 
584
 
 
585
template<class MESH_TYPE, int dim>
 
586
inline APSSCurvatureFeature<MESH_TYPE,dim>::APSSCurvatureFeature(){}
 
587
 
 
588
template<class MESH_TYPE, int dim>
 
589
inline APSSCurvatureFeature<MESH_TYPE,dim>::APSSCurvatureFeature(VertexType& v):pos(v.P()),normal(v.N())
 
590
{
 
591
    normal.Normalize();
 
592
    for(int i=0; i<dim; i++) APSSCurvatureFeature::getNullValue();
 
593
}
 
594
 
 
595
template<class MESH_TYPE, int dim> inline char* APSSCurvatureFeature<MESH_TYPE,dim>::getName()
 
596
{
 
597
    return "APSSCurvatureFeature";
 
598
}
 
599
 
 
600
/* Provides needed attribute to compute feature. A detailed list follows:
 
601
    MM_VERTCURV required by curvature computation
 
602
    MM_VERTQUALITY required by curvature and histogram computation
 
603
    MM_VERTRADIUS required by APSS curvature computation
 
604
    MM_VERTCURVDIR required by APSS curvature computation
 
605
*/
 
606
template<class MESH_TYPE, int dim> inline int APSSCurvatureFeature<MESH_TYPE,dim>::getRequirements()
 
607
{
 
608
    return (MeshModel::MM_VERTCURVDIR | MeshModel::MM_VERTQUALITY | MeshModel::MM_VERTRADIUS);
 
609
}
 
610
 
 
611
template<class MESH_TYPE, int dim> inline bool APSSCurvatureFeature<MESH_TYPE,dim>::isNullValue(float val)
 
612
{
 
613
    return ( ((val==FeatureType::getNullValue()) | (math::IsNAN(val))) );
 
614
}
 
615
 
 
616
template<class MESH_TYPE, int dim> inline float APSSCurvatureFeature<MESH_TYPE,dim>::getNullValue()
 
617
{
 
618
    return -std::numeric_limits<float>::max();
 
619
}
 
620
 
 
621
template<class MESH_TYPE, int dim> inline int APSSCurvatureFeature<MESH_TYPE,dim>::getFeatureDimension()
 
622
{
 
623
    return dim;
 
624
}
 
625
 
 
626
//check persistence beetween scales: return true if description is valid for all scales, false otherwise
 
627
template<class MESH_TYPE, int dim> bool APSSCurvatureFeature<MESH_TYPE,dim>::CheckPersistency(FeatureType f)
 
628
{
 
629
    for(int i = 0; i<FeatureType::getFeatureDimension(); i++)
 
630
    {
 
631
        if( FeatureType::isNullValue(f.description[i]) ) return false;
 
632
    }
 
633
    return true;
 
634
}
 
635
 
 
636
//parameters must be ordered according to smooth iterations
 
637
template<class MESH_TYPE, int dim> void APSSCurvatureFeature<MESH_TYPE,dim>::SetupParameters(ParamType& param)
 
638
{
 
639
    param.add(Parameters::MEAN, 1.0f, 0.3f, 0.9f);
 
640
    param.add(Parameters::MEAN, 0.75f, 0.3f, 0.9f);
 
641
    param.add(Parameters::MEAN, 0.5f, 0.3f, 0.9f);
 
642
    assert(int(param.featureDesc.size())==getFeatureDimension());
 
643
}
 
644
 
 
645
template<class MESH_TYPE, int dim> void APSSCurvatureFeature<MESH_TYPE,dim>::PreCleaning(MeshType& m)
 
646
{
 
647
    //if we are not working on point clouds, clean up faces
 
648
    if(m.fn>0)
 
649
    {
 
650
        tri::Clean<MeshType>::RemoveZeroAreaFace(m);
 
651
        tri::Clean<MeshType>::RemoveDuplicateFace(m);
 
652
        tri::Clean<MeshType>::RemoveDuplicateVertex(m);
 
653
        tri::Clean<MeshType>::RemoveUnreferencedVertex(m);
 
654
    }
 
655
    tri::Allocator<MeshType>::CompactVertexVector(m);
 
656
}
 
657
 
 
658
template<class MESH_TYPE, int dim> bool APSSCurvatureFeature<MESH_TYPE,dim>::ComputeAPSS(MeshType& m, int type = 0, float filterScale = 2.0f, int maxProjIter = 15, float projAcc = 1e-4f, float sphPar = 1.0f, bool accNorm = true, bool selectionOnly=true)
 
659
{
 
660
    // create the MLS surface
 
661
    APSS<MeshType>* mls = new APSS<MeshType>(m);
 
662
 
 
663
    // We require a per vertex radius so as a first thing compute it
 
664
    mls->computeVertexRaddi();
 
665
 
 
666
    mls->setFilterScale(filterScale);
 
667
    mls->setMaxProjectionIters(maxProjIter);
 
668
    mls->setProjectionAccuracy(projAcc);
 
669
    mls->setSphericalParameter(sphPar);
 
670
    mls->setGradientHint(accNorm ? GaelMls::MLS_DERIVATIVE_ACCURATE : GaelMls::MLS_DERIVATIVE_APPROX);
 
671
 
 
672
    uint size = m.vert.size();
 
673
    vcg::Point3f grad;
 
674
    vcg::Matrix33f hess;
 
675
 
 
676
    //computes curvatures and statistics
 
677
    for (unsigned int i = 0; i< size; i++)
 
678
    {
 
679
        if ( (!selectionOnly) || (m.vert[i].IsS()) )
 
680
        {
 
681
            Point3f p = mls->project(m.vert[i].P());
 
682
            float c = 0;
 
683
 
 
684
            if (type==Parameters::APSS) c = mls->approxMeanCurvature(p);
 
685
            else
 
686
            {
 
687
                int errorMask;
 
688
                grad = mls->gradient(p, &errorMask);
 
689
                if (errorMask == MLS_OK && grad.Norm() > 1e-8)
 
690
                {
 
691
                    hess = mls->hessian(p, &errorMask);
 
692
                    implicits::WeingartenMap<float> W(grad,hess);
 
693
 
 
694
                    m.vert[i].PD1() = W.K1Dir();
 
695
                    m.vert[i].PD2() = W.K2Dir();
 
696
                    m.vert[i].K1() =  W.K1();
 
697
                    m.vert[i].K2() =  W.K2();
 
698
 
 
699
                    switch(type){
 
700
                        case Parameters::MEAN: c = W.MeanCurvature(); break;
 
701
                        case Parameters::GAUSSIAN: c = W.GaussCurvature(); break;
 
702
                        case Parameters::K1: c = W.K1(); break;
 
703
                        case Parameters::K2: c = W.K2(); break;
 
704
                        default: assert(0 && "invalid curvature type");
 
705
                    }
 
706
                }
 
707
                assert(!math::IsNAN(c) && "You should never try to compute Histogram with Invalid Floating points numbers (NaN)");
 
708
            }
 
709
            m.vert[i].Q() = c;
 
710
        }
 
711
    }
 
712
 
 
713
    delete mls;
 
714
 
 
715
    return true;
 
716
}
 
717
 
 
718
template<class MESH_TYPE, int dim> bool APSSCurvatureFeature<MESH_TYPE,dim>::HasBeenComputed(MeshType &m)
 
719
{
 
720
    //checks if the attribute exist
 
721
    return tri::HasPerVertexAttribute(m,std::string(FeatureType::getName()));
 
722
}
 
723
 
 
724
template<class MESH_TYPE, int dim> bool APSSCurvatureFeature<MESH_TYPE,dim>::ComputeFeature(MeshType &m, ParamType& param, CallBackPos *cb)
 
725
{
 
726
    //variables needed for progress bar callback
 
727
    float progBar = 0.0f;
 
728
    float offset = 100.0f/((FeatureType::getFeatureDimension())*m.VertexNumber());
 
729
    if(cb) cb(0,"Computing features...");
 
730
 
 
731
    //allocates a custom per vertex attribute in which we can store pointers to features in the heap.
 
732
    PVAttributeHandle fh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(m, true);
 
733
    if(!tri::Allocator<MeshType>::IsValidHandle(m,fh)) return false;
 
734
 
 
735
    //clear the mesh to avoid wrong values during curvature computations
 
736
    PreCleaning(m);
 
737
 
 
738
    //creates a feature for each vertex
 
739
    for(VertexIterator vi = m.vert.begin(); vi!=m.vert.end(); ++vi) fh[vi] = FeatureType(*vi);
 
740
 
 
741
    //loop trough scale levels
 
742
    for (unsigned int i = 0; i<FeatureType::getFeatureDimension(); i++)
 
743
    {
 
744
        //compute the amount of verteces desidered for this scale
 
745
        int targetSize = int(m.VertexNumber()*param.featureDesc[i].scale);
 
746
        //set up the clustering structure and perform a logaritmic search in order to get a number
 
747
        //of clustered vertex very close to the target number
 
748
        tri::Clustering<MeshType, tri::NearestToCenter<MeshType> > ClusteringGrid;
 
749
        bool done = false;
 
750
        bool onlySelected = (i==0) ? false : true; //we subsample the whole mesh just the first time, then we subsample from the previous scale!
 
751
        int size = 10*targetSize, errSize = int(0.02f*targetSize), inf = 0, sup = 0;
 
752
        do{
 
753
            ClusteringGrid.Init(m.bbox,size);
 
754
            ClusteringGrid.AddPointSet(m,onlySelected);
 
755
            int sel = ClusteringGrid.CountPointSet();
 
756
            if(sel<targetSize-errSize){
 
757
                inf = size;
 
758
                if(sup) size+=(sup-inf)/2;
 
759
                else size*=2;
 
760
            }
 
761
            else if(sel>targetSize+errSize){
 
762
                sup = size;
 
763
                if(inf) size-=(sup-inf)/2;
 
764
                else size/=2;
 
765
            }
 
766
            else done=true;
 
767
        }while(!done);
 
768
 
 
769
        //perform clustering. this set some vertesec of m as selected
 
770
        ClusteringGrid.SelectPointSet(m);
 
771
 
 
772
        ComputeAPSS(m);  //compute curvature
 
773
 
 
774
        //compute curvature histogram just on selected verteces
 
775
        Histogram<ScalarType> hist = Histogram<ScalarType>();
 
776
        tri::Stat<MeshType>::ComputePerVertexQualityHistogram(m, hist, true);
 
777
        float vmin = hist.Percentile(param.featureDesc[i].lower_bound); float vmax = hist.Percentile(param.featureDesc[i].upper_bound);
 
778
 
 
779
        //If curvature is beetween bounds and vertex is selected and it is not a boundary, curvature is stored in the feature, otherwise the feature is set to an empty value.
 
780
        for(VertexIterator vi = m.vert.begin(); vi!=m.vert.end(); ++vi)
 
781
        {
 
782
            if( ((*vi).IsS()) & ((!(*vi).IsB()) & ((*vi).Q()<=vmin) | ((*vi).Q()>=vmax) & (!FeatureType::isNullValue((*vi).Q()))) )
 
783
                fh[vi].description[i] = (*vi).Q();
 
784
            else
 
785
                fh[vi].description[i] = FeatureType::getNullValue();
 
786
 
 
787
            if(cb){ progBar+=offset; cb((int)progBar,"Computing features..."); }
 
788
        }
 
789
    }
 
790
    return true;
 
791
}
 
792
 
 
793
template<class MESH_TYPE, int dim>
 
794
MESH_TYPE* APSSCurvatureFeature<MESH_TYPE,dim>::CreateSamplingMesh()
 
795
{
 
796
    MeshType* m = new MeshType();
 
797
    PVAttributeHandle fh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(*m, true);
 
798
    if(!tri::Allocator<MeshType>::IsValidHandle(*m,fh)){
 
799
        if(m) delete m;
 
800
        return NULL;
 
801
    }
 
802
    return m;
 
803
}
 
804
 
 
805
template<class MESH_TYPE, int dim>
 
806
int APSSCurvatureFeature<MESH_TYPE,dim>::SetupSamplingStructures(MeshType& m, typename MeshType::template PerVertexAttributeHandle<FeatureType>& fh, MeshType* samplingMesh, vector<FeatureType*>* vecFeatures)
 
807
{
 
808
    int countFeatures = 0;
 
809
    PVAttributeHandle pmfh;
 
810
    if(samplingMesh){
 
811
        pmfh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(*samplingMesh);
 
812
        if(!tri::Allocator<MeshType>::IsValidHandle(*samplingMesh,pmfh)) return 0;
 
813
    }
 
814
 
 
815
    //fill the vector with all persistent features.
 
816
    for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi){
 
817
        //check persistence beetween scales: if feature is persistent, add a pointer in vecFeatures
 
818
        if( FeatureType::CheckPersistency(fh[vi])){
 
819
            countFeatures++;  //increment counter of valid features
 
820
            if(vecFeatures) vecFeatures->push_back(&(fh[vi]));
 
821
            if(samplingMesh){
 
822
                tri::Allocator<MeshType>::AddVertices(*samplingMesh, 1);
 
823
                samplingMesh->vert.back().ImportLocal(*vi);
 
824
                pmfh[samplingMesh->vert.back()] = fh[vi];
 
825
            }
 
826
        }
 
827
    }
 
828
 
 
829
    return countFeatures;
 
830
}
 
831
 
 
832
template<class MESH_TYPE, int dim>
 
833
bool APSSCurvatureFeature<MESH_TYPE,dim>::Subset(int k, MeshType &m, vector<FeatureType*> &vecSubset, int sampType, CallBackPos *cb)
 
834
{
 
835
    //variables needed for progress bar callback
 
836
    float progBar = 0.0f;
 
837
    float offset = 100.0f/(m.VertexNumber() + k);
 
838
 
 
839
    //if attribute doesn't exist, return; else we can get a handle to the attribute
 
840
    PVAttributeHandle fh = FeatureAlignment<MeshType, FeatureType>::GetFeatureAttribute(m);
 
841
    if(!tri::Allocator<MeshType>::IsValidHandle(m,fh)) return false;
 
842
 
 
843
    //create a vector to hold valid features that later will be sampled
 
844
    vector<FeatureType*>* vecFeatures = NULL;
 
845
    MeshType* poissonMesh = NULL;
 
846
    if(sampType==0) vecFeatures = new vector<FeatureType*>();
 
847
    else poissonMesh = CreateSamplingMesh();
 
848
 
 
849
    //fill the vector with all persistent features.
 
850
    int countFeatures = SetupSamplingStructures(m, fh, poissonMesh, vecFeatures);
 
851
    if(countFeatures<k) k = countFeatures;  //we can't extract more of what we got!
 
852
 
 
853
    //perform different kinds of sampling
 
854
    FeatureType** sampler = NULL;
 
855
    switch(sampType){
 
856
        case 0:{ //uniform sampling: uses vecFeatures
 
857
            sampler = FeatureAlignment<MeshType, FeatureType>::FeatureUniform(*vecFeatures, &k);
 
858
            break;
 
859
        }
 
860
        case 1:{ //poisson disk sampling: uses poissonMesh
 
861
            sampler = FeatureAlignment<MeshType, FeatureType>::FeaturePoisson(*poissonMesh, &k);
 
862
            break;
 
863
        }
 
864
        default: assert(0);
 
865
    }
 
866
 
 
867
    //store features into the returned vector
 
868
    for(int i=0; i<k; ++i){
 
869
        vecSubset.push_back(sampler[i]);
 
870
        if(cb){ progBar+=offset; cb(int(progBar),"Extracting features..."); }
 
871
    }
 
872
 
 
873
    if(vecFeatures) delete vecFeatures;   //clear useless data
 
874
    if(poissonMesh) delete poissonMesh;
 
875
    if(sampler) delete[] sampler;
 
876
 
 
877
    return true;
 
878
 
 
879
    ////UNCOMMENT FOLLOW TO GET A RARE FEATURE SELECTION////
 
880
/*
 
881
    int histSize = 10000;  //number of bins of the histogram
 
882
    vector<FeatureType*>* sorted = new vector<FeatureType*>();  //vector that holds features sorted by bin cont
 
883
 
 
884
    //compute min val e max val between all features; min and max are needed to bound the histogram
 
885
    pair<float,float> minmax = FindMinMax(*vecFeatures, 0);
 
886
 
 
887
    //fill multimap with features sorted by bin count
 
888
    SortByBinCount(*vecFeatures, minmax.first, minmax.second, histSize, *sorted);
 
889
 
 
890
    //select the first k entries from mulptimap, and put them into vecSubset
 
891
    typename vector<FeatureType*>::iterator it = sorted->begin();
 
892
    ///UNCOMMENT THIS LOOP FOR A SORTED SELECTION
 
893
    for(int i=0; i<k & it!=sorted->end(); i++, it++)
 
894
    {
 
895
        (*it)->selected = true;        //mark features as selected
 
896
        vecSubset.push_back(*it);
 
897
    }
 
898
    ///UNCOMMENT THIS LOOP FOR A EQUALIZED SELECTION
 
899
    int off = int(sorted->size()/k);
 
900
    for(it = sorted->begin(); it<sorted->end(); it=it+off)
 
901
    {
 
902
        (*it)->selected = true;        //mark features as selected
 
903
        vecSubset.push_back(*it);
 
904
    }
 
905
 
 
906
    //delete vector and set pointer to NULL for safety
 
907
    if(vecFeatures) delete vecFeatures;   //clear useless data
 
908
    sorted->clear();
 
909
    delete sorted;
 
910
    sorted = NULL;
 
911
    return true;
 
912
*/
 
913
}
 
914
 
 
915
//scan the vector of features and return a pair containig the min and max description values
 
916
template<class MESH_TYPE, int dim> pair<float,float> APSSCurvatureFeature<MESH_TYPE,dim>::FindMinMax(vector<FeatureType*>& vec, int scale)
 
917
{
 
918
    assert(scale>=0 && scale<FeatureType::GetFeatureDimension());
 
919
 
 
920
    typename vector<FeatureType*>::iterator vi;
 
921
    pair<float,float> minmax = make_pair(numeric_limits<float>::max(),-numeric_limits<float>::max());
 
922
    for(vi = vec.begin(); vi!=vec.end(); ++vi)
 
923
    {
 
924
        if( !FeatureType::isNullValue((*vi)->description[scale]))  //invalid values are discarded
 
925
        {
 
926
            if( (*vi)->description[scale] < minmax.first) minmax.first=(*vi)->description[scale];
 
927
            if( (*vi)->description[scale] > minmax.second ) minmax.second=(*vi)->description[scale];
 
928
        }
 
929
    }
 
930
    return minmax;
 
931
}
 
932
 
 
933
//fill a vector that holds features pointers sorted by bin count; i.e first all the very infrequent features and last all the very frequent ones.
 
934
template<class MESH_TYPE, int dim> void APSSCurvatureFeature<MESH_TYPE,dim>::SortByBinCount(vector<FeatureType*>& vecFeatures, float min, float max, int histSize, vector<FeatureType*>& sorted)
 
935
{
 
936
    //vector to hold bins ranges
 
937
    vector<float>* bins = new vector<float>(histSize, 0);
 
938
    //vector to count content of each bin
 
939
    vector<int> *hToC = new vector<int>(histSize, 0);
 
940
    //multimap to keep track of features for each bin
 
941
    multimap<int, FeatureType* >* hToF = new multimap<int, FeatureType* >();
 
942
    //multimap to store pairs of (count, FeatureType*). multimap is naturally sorted by count.
 
943
    multimap<int, FeatureType* >* cToF = new multimap<int, FeatureType* >();
 
944
 
 
945
    //offset beetween min and max is divided into histSize uniform bins
 
946
    for(unsigned int i = 0; i<bins->size(); i++) bins->at(i) = min + i*(max-min)/float(histSize);
 
947
 
 
948
    //build histogram, hToF and hToC; all with just one features scan
 
949
    typename multimap<int, FeatureType* >::iterator hfIt = hToF->begin();
 
950
    typename vector<FeatureType*>::iterator vi;
 
951
    for(vi = vecFeatures.begin(); vi!=vecFeatures.end(); ++vi)
 
952
    {
 
953
        float val = (*vi)->description[0];
 
954
        // lower_bound returns an iterator pointing to the first element "not less than" val, or end() if every element is less than val.
 
955
        typename vector<float>::iterator it = lower_bound(bins->begin(),bins->end(),val);
 
956
        //if val is out of range, skip iteration and take in account next val
 
957
        if(it==bins->begin() || it==bins->end()) continue;
 
958
        //determine in which bin got to insert val
 
959
        int binId = (it - bins->begin())-1;
 
960
        assert(binId>=0);
 
961
        assert (bins->at(binId) < val);
 
962
        assert (val <= bins->at(binId+1));
 
963
 
 
964
        //increment bin count
 
965
        hToC->at(binId)++;
 
966
        //remember in which bin has been inserted this feature
 
967
        hfIt = hToF->insert(hfIt,make_pair(binId, (*vi)));
 
968
    }
 
969
 
 
970
    //delete bins and set pointer to NULL for safety
 
971
    bins->clear();
 
972
    delete bins;
 
973
    bins = NULL;
 
974
 
 
975
    //for all bin index insert in the multimap an entry with key bin count and value a feature. Entries are naturally
 
976
    //sorted by key in the multimap
 
977
    typename multimap<int, FeatureType* >::iterator cfIt = cToF->begin();
 
978
    pair< typename multimap<int, FeatureType* >::iterator, typename multimap<int, FeatureType* >::iterator> range;
 
979
    for(unsigned int i=0; i<hToC->size(); i++)
 
980
    {
 
981
        //if bin count is zero, then skip; nothing to put in the multimap
 
982
        if(hToC->at(i)!=0)
 
983
        {
 
984
            range = hToF->equal_range(i);
 
985
            assert(range.first!=range.second);
 
986
            for (; range.first!=range.second; ++range.first)
 
987
            {
 
988
                cfIt = cToF->insert( cfIt, make_pair(hToC->at(i), (*range.first).second) );
 
989
            }
 
990
        }
 
991
    }
 
992
 
 
993
    typename multimap<int, FeatureType* >::iterator it;
 
994
    for(it = cToF->begin(); it != cToF->end(); it++)
 
995
        sorted.push_back((*it).second);
 
996
 
 
997
    assert(sorted.size()==cToF->size());
 
998
 
 
999
    //delete all ausiliary structures and set pointers to NULL for safety
 
1000
    hToF->clear();
 
1001
    delete hToF;
 
1002
    hToF = NULL;
 
1003
    hToC->clear();
 
1004
    delete hToC;
 
1005
    hToC = NULL;
 
1006
    cToF->clear();
 
1007
    delete cToF;
 
1008
    cToF = NULL;
 
1009
}
 
1010
 
 
1011
template<class MESH_TYPE, int dim> int APSSCurvatureFeature<MESH_TYPE,dim>::SpatialSelection(vector<FeatureType*>& container, vector<FeatureType*>& vec, int k, float radius)
 
1012
{
 
1013
    //variables to manage the kDTree which works on features position
 
1014
    ANNpointArray   dataPts = NULL;     // data points
 
1015
    ANNpoint        queryPnt = NULL;    // query points
 
1016
    ANNkd_tree*     kdTree = NULL;      // search structure
 
1017
 
 
1018
    queryPnt = annAllocPt(3);
 
1019
 
 
1020
    typename vector<FeatureType*>::iterator ci = container.begin();
 
1021
    while(ci != container.end() && vec.size()<k)
 
1022
    {
 
1023
        if(vec.size()==0)
 
1024
        {
 
1025
            vec.push_back(*ci);
 
1026
            (*ci)->selected = true;
 
1027
            ci++;
 
1028
            continue;
 
1029
        }
 
1030
 
 
1031
        if(dataPts){ annDeallocPts(dataPts); dataPts = NULL; }
 
1032
        if(kdTree){ delete kdTree; kdTree = NULL; }
 
1033
        dataPts = annAllocPts(vec.size(), 3);
 
1034
 
 
1035
        for(int i = 0; i < vec.size(); i++)
 
1036
        {
 
1037
            for (int j = 0; j < 3; j++)
 
1038
            {
 
1039
                (dataPts[i])[j] = (ANNcoord)(vec.at(i)->pos[j]);
 
1040
            }
 
1041
        }
 
1042
        //build search structure
 
1043
        kdTree = new ANNkd_tree(dataPts,vec.size(),3);
 
1044
 
 
1045
        for (int j = 0; j < 3; j++)
 
1046
        {
 
1047
            queryPnt[j] = (ANNcoord)((*ci)->pos[j]);
 
1048
        }
 
1049
 
 
1050
        //if there aren't features yet selected in the range distance
 
1051
        if(!kdTree->annkFRSearch(queryPnt, math::Sqr(radius), 0, NULL, NULL, 0.0))
 
1052
        {
 
1053
             vec.push_back(*ci);
 
1054
             (*ci)->selected = true;
 
1055
        }
 
1056
        ci++;
 
1057
    }
 
1058
 
 
1059
    if(dataPts){ annDeallocPts(dataPts); dataPts = NULL; }
 
1060
    if(queryPnt){ annDeallocPt(queryPnt); queryPnt = NULL; }
 
1061
    if(kdTree){ delete kdTree; kdTree = NULL; }
 
1062
 
 
1063
    return vec.size();
 
1064
}