00001 #ifndef NCUT_CORE_H
00002 #define NCUT_CORE_H 1
00003
00004 #include "motionProfile.h"
00005
00006 #define ARPACK_SILENT_MODE 1
00007
00008 #include <arbsmat.h>
00009 #include <arbgsym.h>
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 namespace ncut
00020 {
00021
00031 class Ncut
00032 {
00033 public:
00034
00035
00036
00041 inline const Matrix* similarity() const
00042 {return similarity_;};
00048 inline double assoc(unsigned int i) const
00049 {return assoc_[i];};
00054 inline const double* assoc() const
00055 {return assoc_;};
00060 inline unsigned int dim() const
00061 {return dim_;};
00067 inline double eigVal(unsigned int i) const
00068 {return eigVal_[i];};
00073 inline const double* eigVal() const
00074 {return eigVal_;};
00080 inline double eigVec(unsigned int i) const
00081 {return eigVec_[i];};
00086 inline const double* eigVec() const
00087 {return eigVec_;};
00092 inline unsigned int nConv() const
00093 {return nConv_;};
00094
00095
00096
00109 int calculate(unsigned int nEigVecs, double* initVec = NULL);
00110
00111
00112
00119 Ncut(const Matrix* similarity);
00120
00128 Ncut(const Ncut& clone);
00129
00133 virtual ~Ncut();
00134
00141 virtual Ncut& operator=(const Ncut& clone);
00142
00143 protected:
00144
00145 const Matrix* similarity_;
00147 unsigned int dim_;
00150 double* assoc_;
00152 double* eigVal_;
00154 double* eigVec_;
00156 unsigned int nConv_;
00159 private:
00160
00164 void cleanup();
00165 };
00166
00188 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
00189 class NcutNode
00190 {
00191 public:
00192
00193
00194
00199 inline const greedy_ptr<const Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >
00200 profile() const
00201 {return profile_;};
00206 inline const std::map<MSK_ELM_TYPE, unsigned int>* indices() const
00207 {return &indices_;};
00212 inline const greedy_ptr<const Ncut> ncut() const
00213 {return ncut_;};
00218 inline double eigVal() const
00219 {return eigVal_;};
00224 inline const double* eigVec() const
00225 {return eigVec_;};
00230 inline double splitPoint() const
00231 {return splitPoint_;};
00236 inline double eigVecNum() const
00237 {return eigVecNum_;};
00242 inline const greedy_ptr<const NcutNode> leftNode() const
00243 {return leftNode_;};
00248 inline const greedy_ptr<const NcutNode> rightNode() const
00249 {return rightNode_;};
00250
00251
00252
00274 int calculate(unsigned int maxDepth, unsigned int depth,
00275 unsigned int nEigVecs, unsigned int eigvn,
00276 double* initVec = NULL);
00277
00287 bool isLeft(unsigned int idx) const;
00288
00289
00290
00304 NcutNode(std::auto_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& profile,
00305 const std::map<MSK_ELM_TYPE, unsigned int>& indices,
00306 greedy_ptr<Ncut>& ncut, const Setting* setting);
00307
00321 NcutNode(greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& profile,
00322 const std::map<MSK_ELM_TYPE, unsigned int>& indices,
00323 greedy_ptr<Ncut>& ncut, const Setting* setting);
00324
00333 NcutNode(const NcutNode& clone);
00334
00347 NcutNode(const NcutNode& clone,
00348 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& pclone,
00349 greedy_ptr<Ncut>& nclone);
00350
00356 virtual ~NcutNode();
00357
00367 virtual NcutNode& operator=(const NcutNode& clone);
00368
00369 protected:
00370
00371
00372 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > profile_;
00374 std::map<MSK_ELM_TYPE, unsigned int> indices_;
00378
00379 greedy_ptr<Ncut> ncut_;
00381 double eigVal_;
00383 const double* eigVec_;
00385 double splitPoint_;
00386 unsigned int eigVecNum_;
00389
00390 greedy_ptr<NcutNode> leftNode_;
00394 greedy_ptr<NcutNode> rightNode_;
00399 const Setting* setting_;
00401 private:
00402
00418 NcutNode&
00419 copy(const NcutNode& clone,
00420 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& pclone,
00421 greedy_ptr<Ncut>& nclone);
00422 };
00423
00436 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
00437 class Segmentation
00438 {
00439 public:
00440
00441
00442
00448 inline Setting* setting()
00449 {return &setting_;};
00450
00451
00452
00458 inline const Image<IMG_ELM_TYPE>* frame(unsigned int num) const
00459 {return elmProfile_.frame(num);};
00465 inline const Image<MSK_ELM_TYPE>* elmMask(unsigned int num) const
00466 {return elmProfile_.mask(num);};
00473 inline const MSK_ELM_TYPE& elmMask(unsigned int num, unsigned int pxl) const
00474 {return elmProfile_.mask(num, pxl);};
00481 inline const IMG_ELM_TYPE& frameMask(unsigned int num,
00482 unsigned int pxl) const
00483 {return elmProfile_.frame(num, pxl);};
00488 inline unsigned int width() const
00489 {return elmProfile_.width();};
00494 inline unsigned int height() const
00495 {return elmProfile_.height();};
00500 inline unsigned int nChl() const
00501 {return elmProfile_.nChl();};
00506 inline unsigned int nElm() const
00507 {return elmProfile_.nElm();};
00508
00514 inline double locationX(unsigned int seg) const
00515 {return segProfile_.locationX(seg);};
00521 inline double locationY(unsigned int seg) const
00522 {return segProfile_.locationY(seg);};
00523
00528 inline unsigned int nSeg() const
00529 {return segProfile_.nSeg();};
00534 inline unsigned int nFrm() const
00535 {return elmProfile_.nFrm();};
00540 inline unsigned int wSize() const
00541 {return elmProfile_.wSize();};
00542
00547 inline const std::map<MSK_ELM_TYPE, unsigned int>* indices() const
00548 {return segProfile_.indices();}
00549
00555 inline unsigned int idx(const MSK_ELM_TYPE& segId) const
00556 {return segProfile_.idx(segId);};
00562 inline const MSK_ELM_TYPE& id(unsigned int segIdx) const
00563 {return segProfile_.id(segIdx);};
00569 inline const MSK_ELM_TYPE& segMask(const MSK_ELM_TYPE& elmId) const
00570 {return segProfile_.mask(elmId);};
00576 inline const std::vector<unsigned int>* pixels(unsigned int elmIdx) const
00577 {return elmProfile_.pixels(elmIdx);};
00583 inline double nPixels(int seg) const
00584 {return segProfile_.nPixels(seg);};
00590 inline const std::vector<unsigned int>* elements(unsigned int segIdx) const
00591 {return segProfile_.elements(segIdx);};
00597 inline unsigned int frameOf(unsigned int elmIdx) const
00598 {return elmProfile_.frameOf(elmIdx);};
00604 inline unsigned int age(const MSK_ELM_TYPE& segId) const
00605 {return segProfile_.age(segId);};
00610 inline const std::map<MSK_ELM_TYPE, unsigned int>* age() const
00611 {return segProfile_.age();};
00612
00618 inline const bool inactive(const MSK_ELM_TYPE& segId) const
00619 {return segProfile_.inactive(segId);};
00624 inline const std::map<MSK_ELM_TYPE, bool>* inactive() const
00625 {return segProfile_.inactive();};
00631 inline const MSK_ELM_TYPE& simSeg(const MSK_ELM_TYPE& segId) const
00632 {return segProfile_.simSeg(segId);};
00637 inline const NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>* dendogram() const
00638 {return dendogram_;};
00643 inline const ElementProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>* elmProfile() const
00644 {return &elmProfile_;};
00649 inline const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>* segProfile() const
00650 {return segProfile_;};
00656 inline const Setting* setting() const
00657 {return &setting_;};
00658
00672 int push(const Image<IMG_ELM_TYPE>* frame,
00673 const Image<MSK_ELM_TYPE>* mask);
00674
00685 Segmentation(unsigned int wSize = 3, Setting setting = Setting());
00686
00692 Segmentation(const Segmentation& clone);
00693
00697 virtual ~Segmentation();
00698
00705 virtual Segmentation& operator=(const Segmentation& clone);
00706
00707 protected:
00708
00709 Setting setting_;
00711 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>* dendogram_;
00714 ElementProfile<IMG_ELM_TYPE, MSK_ELM_TYPE> elmProfile_;
00716 SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE> segProfile_;
00719 private:
00720
00730 void getInfo(const NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>* node,
00731 std::vector<std::vector<MSK_ELM_TYPE> >& elements);
00732
00733
00745 std::vector<std::vector<MSK_ELM_TYPE> >
00746 merge(const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>& profile);
00747
00761 void connected(const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>& segProfile,
00762 unsigned int startElmIdx, std::vector<bool>& reached);
00763
00764
00776 std::vector<std::vector<MSK_ELM_TYPE> >
00777 split(const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>& profile);
00778
00782 void cleanup();
00783 };
00784
00785
00786
00787
00788 int Ncut::calculate(unsigned int nEigVecs, double* initVec)
00789 {
00790 nEigVecs = 4;
00791 if (similarity_ == NULL)
00792 throw ncutException("invalid similarity matrix (NULL)");
00793 if (dim_ == 0 || similarity_->dim() != dim_)
00794 throw ncutException("invalid dimension, must be > 0 and equal to "
00795 "similarity matrix dimension");
00796
00797
00798 cleanup();
00799
00800
00801 Matrix assocMatrix(dim_);
00802 for (unsigned int i=0; i<dim_; i++)
00803 for (unsigned int j=0; j<dim_; j++)
00804 assocMatrix.val(i,i) += similarity_->val(i,j);
00805
00806
00807 assoc_ = new double [dim_];
00808 for (unsigned int i=0; i<dim_; i++)
00809 assoc_[i] = assocMatrix.val(i,i);
00810
00811
00812 Matrix assocSubSimil(assocMatrix);
00813 assocSubSimil -= *similarity_;
00814
00815
00816 int assocSubSimilN;
00817 int assocSubSimilNsdiag;
00818 auto_array<double> assocSubSimilBand(assocSubSimil.band(assocSubSimilN,
00819 assocSubSimilNsdiag));
00820 int assocN;
00821 int assocNsdiag;
00822 auto_array<double> assocBand(assocMatrix.band(assocN, assocNsdiag));
00823
00824 try {
00825
00826 ARbdSymMatrix<double> assocSubSimilARP(assocSubSimilN, assocSubSimilNsdiag,
00827 assocSubSimilBand.get());
00828 ARbdSymMatrix<double> assocARP(assocN, assocNsdiag, assocBand.get());
00829
00830
00831 int ncv = ((2*nEigVecs)+1<dim_-1) ? (2*nEigVecs)+1 : dim_;
00832
00833 ARluSymGenEig<double> problem(nEigVecs, assocSubSimilARP, assocARP, "SA",
00834 ncv, 0.0, nEigVecs*100000, initVec);
00835
00836 for (unsigned int i=0; i<dim_; i++)
00837 {
00838 for (unsigned int j=0; j<dim_; j++)
00839 {
00840 std::cout << "," << assocSubSimil.val(i,j) << " ";
00841 }
00842 std::cout << std::endl << std::flush;
00843 }
00844
00845
00846 eigVal_ = new double[nEigVecs];
00847 eigVec_ = new double[nEigVecs*dim_];
00848 nConv_ = problem.EigenValVectors(eigVec_, eigVal_);
00849
00850
00851
00852
00853
00854
00855
00856 }
00857 catch (ArpackError& e) {
00858 throw ncutException(
00859 (std::string("ArpackError: error code ") +
00860 std::string(1,e.Status()) ).c_str() );
00861 }
00862
00863 return 0;
00864 }
00865
00866 Ncut::Ncut(const Matrix* similarity)
00867 : similarity_(similarity),
00868 dim_( (similarity==NULL) ? 0 : similarity->dim() ),
00869 assoc_(NULL), eigVal_(NULL), eigVec_(NULL), nConv_(0)
00870 {
00871 if (similarity_ == NULL)
00872 throw ncutException("invalid similarity matrix (NULL)");
00873 return;
00874 }
00875
00876 Ncut::Ncut(const Ncut& clone)
00877 : similarity_(NULL), dim_(0),
00878 assoc_(NULL), eigVal_(NULL), eigVec_(NULL), nConv_(0)
00879 {
00880 *this = clone;
00881 }
00882
00883 void Ncut::cleanup()
00884 {
00885 if (assoc_ != NULL) {
00886 delete [] assoc_;
00887 assoc_ = NULL;
00888 }
00889 if (eigVal_ != NULL) {
00890 delete [] eigVal_;
00891 eigVal_ = NULL;
00892 }
00893 if (eigVec_ != NULL) {
00894 delete [] eigVec_;
00895 eigVec_ = NULL;
00896 }
00897 }
00898
00899 Ncut::~Ncut()
00900 {
00901 cleanup();
00902 }
00903
00904 Ncut& Ncut::operator=(const Ncut& clone)
00905 {
00906
00907 if (this == &clone) return *this;
00908
00909
00910 cleanup();
00911
00912
00913 this->similarity_ = clone.similarity_;
00914 this->dim_ = clone.dim_;
00915 this->nConv_ = clone.nConv_;
00916
00917
00918 assoc_ = new double[clone.dim_];
00919 eigVal_ = new double[clone.nConv_];
00920 eigVec_ = new double[clone.nConv_*clone.dim_];
00921
00922
00923 for (unsigned int i=0; i<clone.dim_; i++) {
00924 this->assoc_[i] = clone.assoc_[i];
00925 }
00926
00927
00928 for (unsigned i=0; i<clone.nConv_; i++) {
00929 this->eigVal_[i] = clone.eigVal_[i];
00930
00931 for (unsigned j=0; j<clone.dim_; j++) {
00932 this->eigVec_[i*j] = clone.eigVec_[i*j];
00933 }
00934 }
00935
00936 return *this;
00937 }
00938
00939 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
00940 bool NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::isLeft(unsigned int idx) const
00941 {
00942 if (eigVec_ == NULL || ncut_.get() == NULL || idx >= ncut_->dim() || idx < 0)
00943 throw ncutException("invalid element index");
00944 if (eigVec_[idx] > splitPoint_)
00945 return true;
00946 else
00947 return false;
00948 }
00949
00950
00951 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
00952 int NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
00953 calculate(unsigned int maxDepth, unsigned int depth,
00954 unsigned int nEigVecs, unsigned int eigvn,
00955 double* initVec)
00956 {
00957 #ifdef debugout
00958 std::cout << "::calculate depth " << depth << " eigvn " << eigvn
00959 << " indices_.size() " << indices_.size() << std::endl << std::flush;
00960 #endif
00961
00962 if (profile_.get() == NULL)
00963 throw ncutException("invalid profile (NULL)");
00964 if (depth >= maxDepth ||
00965 eigvn >= nEigVecs ||
00966 eigvn < 1)
00967 throw ncutException("invalid arguments");
00968
00969 unsigned int nStep = setting_->number_of_eigenvector_splittest_steps;
00970 if (nStep < 1)
00971 throw ncutException("number of split-test steps invalid, must be >= 1");
00972
00973
00974
00975
00976 eigVecNum_ = eigvn;
00977
00978
00979 if (eigvn == 1 || ncut_.get() == NULL) {
00980 ncut_ = new Ncut(profile_->similarity());
00981 try {
00982 ncut_->calculate(nEigVecs, initVec);
00983 }
00984 catch (ncutException& e) {
00985 #ifdef debugout
00986 std::cerr << e.what() << std::endl << std::flush;
00987 #endif
00988 return 1;
00989 }
00990 }
00991
00992
00993 if (ncut_->nConv() <= eigvn) return 1;
00994
00995
00996 eigVal_ = ncut_->eigVal(eigvn);
00997 eigVec_ = &( (ncut_->eigVec())[eigvn*ncut_->dim()] );
00998
00999
01000 typename std::map<MSK_ELM_TYPE, unsigned int>::iterator elm;
01001 typename std::map<MSK_ELM_TYPE, unsigned int>::iterator elm2;
01002 double minVal, maxVal;
01003 minVal = maxVal = eigVec_[0];
01004 for (unsigned int i=0; i<ncut_->dim(); i++) {
01005 if (eigVec_[i] < minVal) minVal = eigVec_[i];
01006 if (eigVec_[i] > maxVal) maxVal = eigVec_[i];
01007 }
01008 double step = (maxVal - minVal) / (double) nStep;
01009 if (step <= 0.0) return 1;
01010
01011 double bestNcut = 2.0;
01012 double assoc_A, assoc_B, cut, ncutValue;
01013 for (double split = minVal+(step/2.0); split<maxVal; split+=step) {
01014 assoc_A = assoc_B = cut = 0.0;
01015 for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01016 if (eigVec_[elm->second] < split)
01017 assoc_A += ncut_->assoc(elm->second);
01018 else
01019 assoc_B += ncut_->assoc(elm->second);
01020
01021 for (elm2=elm, elm2++; elm2!=indices_.end(); elm2++) {
01022 if ((eigVec_[elm->second] < split) !=
01023 (eigVec_[elm2->second] < split)) {
01024 cut += profile_->similarity()->val(elm->second,elm2->second);
01025 }
01026 }
01027 }
01028 ncutValue = (cut/assoc_A) + (cut/assoc_B);
01029 if (ncutValue < bestNcut) {
01030 splitPoint_ = split;
01031 bestNcut = ncutValue;
01032 }
01033 }
01034 #ifdef debugout
01035 std::cout << "splitPoint - " << splitPoint_ << std::endl << std::flush;
01036 #endif
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113 #ifdef debugout
01114 std::cout << "---splitPoint " << splitPoint_ << std::endl << std::flush;
01115 #endif
01116
01117
01118 if (depth >= maxDepth-1) return 0;
01119
01120
01121
01122
01123
01124 unsigned int left_eigvn = (eigvn*2);
01125 unsigned int right_eigvn = left_eigvn+1;
01126 if (left_eigvn >= nEigVecs || right_eigvn >= nEigVecs);
01127 left_eigvn = right_eigvn = 1;
01128
01129
01130 std::map<MSK_ELM_TYPE, unsigned int> left_indices;
01131 std::map<MSK_ELM_TYPE, unsigned int> right_indices;
01132
01133 if (left_eigvn != 1 && right_eigvn != 1) {
01134
01135
01136
01137
01138 for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01139 if (eigVec_[elm->second] > splitPoint_)
01140 left_indices[elm->first] = elm->second;
01141 else if (eigVec_[elm->second] < splitPoint_)
01142 right_indices[elm->first] = elm->second;
01143 }
01144
01145
01146 if (left_indices.empty() || right_indices.empty()) return 1;
01147
01148
01149 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > left_profile;
01150 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > right_profile;
01151 left_profile = right_profile = profile_.get();
01152
01153
01154 leftNode_ = new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(left_profile,
01155 left_indices,
01156 ncut_, setting_);
01157 int r = leftNode_->calculate(maxDepth,depth+1,nEigVecs,left_eigvn);
01158 if (r == 1) leftNode_.reset(NULL);
01159
01160
01161 rightNode_ = new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(right_profile,
01162 right_indices,
01163 ncut_, setting_);
01164 r = rightNode_->calculate(maxDepth,depth+1,nEigVecs,right_eigvn);
01165 if (r == 1) rightNode_.reset(NULL);
01166 }
01167 else {
01168
01169
01170
01171
01172 std::list<unsigned int> left_elements;
01173 std::list<unsigned int> right_elements;
01174 for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01175
01176 if (eigVec_[elm->second] > splitPoint_)
01177 left_elements.push_back(elm->second);
01178 else if (eigVec_[elm->second] < splitPoint_)
01179 right_elements.push_back(elm->second);
01180 }
01181
01182
01183 if (left_elements.empty() || right_elements.empty()) return 1;
01184
01185
01186 std::vector<unsigned int> left_indexMap;
01187 std::vector<unsigned int> right_indexMap;
01188 Matrix left_sim(0);
01189 Matrix right_sim(0);
01190 Matrix dummy(0);
01191 profile_->similarity()->split(left_sim, dummy, left_indexMap,
01192 left_elements);
01193 profile_->similarity()->split(right_sim, dummy, right_indexMap,
01194 right_elements);
01195
01196
01197 for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01198 if (eigVec_[elm->second] > splitPoint_) {
01199 left_indices[elm->first] = left_indexMap[elm->second];
01200 }
01201 else if (eigVec_[elm->second] < splitPoint_) {
01202 right_indices[elm->first] = right_indexMap[elm->second];
01203 }
01204 }
01205
01206
01207 std::vector<MSK_ELM_TYPE> left_identifiers(left_indices.size());
01208 std::vector<MSK_ELM_TYPE> right_identifiers(right_indices.size());
01209 typename std::map<MSK_ELM_TYPE, unsigned int>::iterator indicesItr;
01210 for (indicesItr = left_indices.begin();
01211 indicesItr != left_indices.end(); indicesItr++) {
01212 left_identifiers[indicesItr->second] = indicesItr->first;
01213 }
01214 for (indicesItr = right_indices.begin();
01215 indicesItr != right_indices.end(); indicesItr++) {
01216 right_identifiers[indicesItr->second] = indicesItr->first;
01217 }
01218
01219
01220 std::auto_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >
01221 left_profile(new Profile<IMG_ELM_TYPE,MSK_ELM_TYPE>(left_sim,
01222 left_indices,
01223 left_identifiers,
01224 setting_));
01225 std::auto_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >
01226 right_profile(new Profile<IMG_ELM_TYPE,MSK_ELM_TYPE>(right_sim,
01227 right_indices,
01228 right_identifiers,
01229 setting_));
01230
01231
01232 leftNode_ = new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(left_profile,
01233 left_indices,
01234 ncut_, setting_);
01235 int r = leftNode_->calculate(maxDepth,depth+1,nEigVecs,left_eigvn);
01236 if (r == 1) leftNode_.reset(NULL);
01237
01238
01239 rightNode_ = new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(right_profile,
01240 right_indices,
01241 ncut_, setting_);
01242 r = rightNode_->calculate(maxDepth,depth+1,nEigVecs,right_eigvn);
01243 if (r == 1) rightNode_.reset(NULL);
01244 }
01245
01246 return 0;
01247 }
01248
01249 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01250 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01251 NcutNode(std::auto_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& profile,
01252 const std::map<MSK_ELM_TYPE, unsigned int>& indices,
01253 greedy_ptr<Ncut>& ncut, const Setting* setting)
01254 : profile_(profile), indices_(indices), ncut_(ncut), eigVal_(0),
01255 eigVec_(NULL), splitPoint_(0.0), eigVecNum_(0),
01256 leftNode_(NULL), rightNode_(NULL), setting_(setting)
01257 {
01258 if (profile_.get() == NULL)
01259 throw ncutException("invalid profile (NULL)");
01260
01261 return;
01262 }
01263
01264 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01265 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01266 NcutNode(greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& profile,
01267 const std::map<MSK_ELM_TYPE, unsigned int>& indices,
01268 greedy_ptr<Ncut>& ncut, const Setting* setting)
01269 : profile_(profile), indices_(indices), ncut_(ncut), eigVal_(0),
01270 eigVec_(NULL), splitPoint_(0.0), eigVecNum_(0),
01271 leftNode_(NULL), rightNode_(NULL), setting_(setting)
01272 {
01273 if (profile_.get() == NULL)
01274 throw ncutException("invalid profile (NULL)");
01275
01276 return;
01277 }
01278
01279 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01280 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::~NcutNode()
01281 {
01282 return;
01283 }
01284
01285 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01286 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01287 NcutNode(const NcutNode& clone)
01288 : profile_(NULL), indices_(), ncut_(NULL), eigVal_(0),
01289 eigVec_(NULL), splitPoint_(0.0), eigVecNum_(0),
01290 leftNode_(NULL), rightNode_(NULL), setting_(clone.setting_)
01291 {
01292 *this = clone;
01293 }
01294
01295 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01296 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01297 NcutNode(const NcutNode& clone,
01298 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& pclone,
01299 greedy_ptr<Ncut>& nclone)
01300 : profile_(NULL), indices_(), ncut_(NULL), eigVal_(0),
01301 eigVec_(NULL), splitPoint_(0.0), eigVecNum_(0),
01302 leftNode_(NULL), rightNode_(NULL), setting_(clone.setting_)
01303 {
01304 copy(clone,pclone,nclone);
01305 }
01306
01307 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01308 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>& NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01309 operator=(const NcutNode& clone)
01310 {
01311 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > pclone(NULL);
01312 greedy_ptr<Ncut> nclone(NULL);
01313 return copy(clone,pclone,nclone);
01314 }
01315
01316 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01317 NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>& NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01318 copy(const NcutNode& clone,
01319 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> >& pclone,
01320 greedy_ptr<Ncut>& nclone)
01321 {
01322
01323 this->indices_ = clone.indices_;
01324 this->eigVal_ = clone.eigVal_;
01325 this->splitPoint_ = clone.splitPoint_;
01326 this->eigVecNum_ = clone.eigVecNum_;
01327
01328
01329
01330
01331
01332 if (clone.profile_.owns())
01333 this->profile_ = new Profile<IMG_ELM_TYPE, MSK_ELM_TYPE>(*clone.profile_);
01334 else if (pclone.get() != NULL)
01335 this->profile_ = pclone;
01336 else
01337 this->profile_ = clone.profile_;
01338
01339
01340
01341
01342
01343 if (clone.ncut_.owns()) {
01344 this->ncut_ = new Ncut(*clone.ncut_);
01345 this->eigVec_ =
01346 &(this->ncut_->eigVec()[this->eigVecNum_*this->ncut_->dim()]);
01347 }
01348 else if (nclone.get() != NULL) {
01349 this->ncut_ = nclone;
01350 this->eigVec_ =
01351 &(this->ncut_->eigVec()[this->eigVecNum_*this->ncut_->dim()]);
01352 }
01353 else {
01354 this->ncut_ = clone.ncut_;
01355 this->eigVec_ = clone.eigVec_;
01356 }
01357
01358
01359
01360 if (clone.leftNode_.get() != NULL)
01361 this->leftNode_ =
01362 new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(*clone.leftNode_,
01363 this->profile_,
01364 this->ncut_);
01365 if (clone.rightNode_.get() != NULL)
01366 this->rightNode_ =
01367 new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(*clone.rightNode_,
01368 this->profile_,
01369 this->ncut_);
01370
01371
01372 this->setting_ = clone.setting_;
01373
01374 return *this;
01375 }
01376
01377 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01378 void Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01379 getInfo(const NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>* node,
01380 std::vector<std::vector<MSK_ELM_TYPE> >& elements)
01381 {
01382 if (node == NULL || node->eigVec() == NULL)
01383 throw ncutException("invalid node, node corrupt");
01384
01385
01386 const NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>* leftChild =
01387 node->leftNode().get();
01388 const NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>* rightChild =
01389 node->rightNode().get();
01390
01391 if (leftChild == NULL || rightChild == NULL) {
01392
01393
01394
01395
01396 unsigned int leftSegNum;
01397 unsigned int rightSegNum;
01398 if (leftChild == NULL) {
01399 leftSegNum = elements.size();
01400 elements.push_back(std::vector<MSK_ELM_TYPE>());
01401 }
01402 if (rightChild == NULL) {
01403 rightSegNum = elements.size();
01404 elements.push_back(std::vector<MSK_ELM_TYPE>());
01405 }
01406
01407
01408 unsigned int segNum;
01409 typename std::map<MSK_ELM_TYPE, unsigned int>::const_iterator elmItr;
01410 for (elmItr = node->indices()->begin();
01411 elmItr != node->indices()->end(); elmItr++) {
01412
01413 if (node->eigVec()[elmItr->second] > node->splitPoint()) {
01414 if (leftChild != NULL) continue;
01415 segNum = leftSegNum;
01416 }
01417 else {
01418 if (rightChild != NULL) continue;
01419 segNum = rightSegNum;
01420 }
01421
01422 elements[segNum].push_back(elmItr->first);
01423 }
01424 }
01425
01426
01427 if (leftChild != NULL)
01428 getInfo(leftChild, elements);
01429 if (rightChild != NULL)
01430 getInfo(rightChild, elements);
01431 }
01432
01433 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01434 std::vector<std::vector<MSK_ELM_TYPE> >
01435 Segmentation<IMG_ELM_TYPE,MSK_ELM_TYPE>::
01436 merge(const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>& profile)
01437 {
01438 unsigned int minSegCount = setting_.minimum_segment_count;
01439 double mergeThreshold = setting_.merge_threshold;
01440
01441 std::vector<std::vector<MSK_ELM_TYPE> > elements = *profile.elements();
01442
01443 if (profile.nSeg() <= 1)
01444 return elements;
01445 if (profile.similarity()->dim() != profile.nSeg())
01446 throw ncutException(
01447 "invalid segment profile (similarity matrix not calculated)." );
01448 if (profile.elmProfile() == NULL)
01449 throw ncutException("invalid element profile (NULL)");
01450 if (minSegCount <= 1)
01451 throw ncutException("invalid minimum segment count, must be >= 2");
01452 if (mergeThreshold < 0.0 || mergeThreshold >= 500.0)
01453 throw ncutException("invalid merge threshold, must be in [0,500)");
01454
01455
01456
01457
01458
01459
01460 Matrix difference(profile.nSeg());
01461 try {
01462 for (unsigned int seg1=0; seg1<profile.nSeg(); seg1++) {
01463 for (unsigned int seg2=seg1+1; seg2<profile.nSeg(); seg2++) {
01464 difference.val(seg2,seg1) = difference.val(seg1,seg2) =
01465 profile.dif(seg1, seg2);
01466 }
01467 }
01468 }
01469 catch (ncutException& e) {
01470 throw ncutCritical(
01471 (std::string("Error: Segmentation left in inconsistent state: ")+
01472 e.what() ).c_str() );
01473 }
01474
01475
01476 bool active [profile.nSeg()];
01477 for (unsigned int seg=0; seg<profile.nSeg(); seg++)
01478 active[seg] = true;
01479 unsigned int mergeSeg1;
01480 unsigned int mergeSeg2;
01481 unsigned int segCount = profile.nSeg();
01482 while(segCount >= minSegCount) {
01483
01484
01485 double minDif = 600.0;
01486 for (unsigned int seg1=0; seg1<profile.nSeg(); seg1++) {
01487 if (!active[seg1]) continue;
01488 for (unsigned int seg2=seg1+1; seg2<profile.nSeg(); seg2++) {
01489 if (!active[seg2]) continue;
01490 if (difference.val(seg1,seg2) < minDif) {
01491 mergeSeg1 = seg1;
01492 mergeSeg2 = seg2;
01493 minDif = difference.val(seg1,seg2);
01494 }
01495 }
01496 }
01497
01498
01499 #ifdef debugout
01500 std::cout << "would merge " << minDif << std::endl << std::flush;
01501 #endif
01502
01503 if (minDif > mergeThreshold) break;
01504 #ifdef debugout
01505 std::cout << " ... and will" << minDif << std::endl << std::flush;
01506 #endif
01507
01508
01509 for (unsigned int elm=0; elm<elements[mergeSeg2].size(); elm++)
01510 elements[mergeSeg1].push_back(elements[mergeSeg2][elm]);
01511 active[mergeSeg2] = false;
01512
01513
01514 try {
01515 for (unsigned int seg=0; seg<elements.size(); seg++) {
01516 if ((!active[seg]) || seg == mergeSeg1) continue;
01517
01518 difference.val(seg,mergeSeg1) = difference.val(mergeSeg1,seg) =
01519 profile.dif(mergeSeg1, seg, &elements);
01520 }
01521 }
01522 catch (ncutException& e) {
01523 throw ncutCritical(
01524 (std::string("Error: Segmentation left in inconsistent state: ")+
01525 e.what() ).c_str() );
01526 }
01527
01528 segCount--;
01529 }
01530
01531
01532
01533 unsigned int segIdx;
01534 typename std::vector<std::vector<MSK_ELM_TYPE> >::iterator seg;
01535 for (seg=elements.begin(), segIdx=0; seg!=elements.end(); segIdx++) {
01536 if (!active[segIdx])
01537 seg = elements.erase(seg);
01538 else
01539 seg++;
01540 }
01541
01542 return elements;
01543 }
01544
01545 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01546 void Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01547 connected(const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>& segProfile,
01548 unsigned int startElmIdx, std::vector<bool>& reached)
01549 {
01550 if (startElmIdx >= segProfile.nElm())
01551 throw ncutException(
01552 "invalid start element index or invalid segment profile" );
01553
01554 const Matrix* neighbors = segProfile.elmProfile()->neighbors();
01555 unsigned int nElm = segProfile.nElm();
01556
01557 unsigned int startSegIdx = segProfile.idxMask(startElmIdx);
01558
01559 reached[startElmIdx] = true;
01560
01561 for (unsigned int elm=0; elm<nElm; elm++) {
01562 if (neighbors->val(startElmIdx,elm) > 0.5 && !reached[elm]) {
01563
01564 if (segProfile.idxMask(elm) == startSegIdx ) {
01565
01566 connected(segProfile, elm, reached);
01567 }
01568 }
01569 }
01570 }
01571
01572 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01573 std::vector<std::vector<MSK_ELM_TYPE> >
01574 Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01575 split(const SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE>& profile)
01576 {
01577 if (profile.elmProfile() == NULL)
01578 throw ncutException("invalid element profile (NULL)");
01579
01580
01581 std::vector<std::vector<MSK_ELM_TYPE> > elements;
01582 std::map<MSK_ELM_TYPE, unsigned int>
01583 unassigned(*profile.elmProfile()->indices());
01584 std::vector<bool> reached(profile.nElm(), false);
01585 typename std::map<MSK_ELM_TYPE, unsigned int>::iterator elmItr;
01586 typename std::map<MSK_ELM_TYPE, unsigned int>::iterator eraseItr;
01587 while (!unassigned.empty()) {
01588
01589
01590
01591 connected(profile, unassigned.begin()->second, reached);
01592
01593
01594
01595 elements.push_back(std::vector<MSK_ELM_TYPE>());
01596 for (elmItr = unassigned.begin(); elmItr!=unassigned.end();) {
01597 if (reached[elmItr->second]) {
01598 elements.back().push_back(elmItr->first);
01599 eraseItr = elmItr++;
01600 unassigned.erase(eraseItr);
01601 }
01602 else {
01603 elmItr++;
01604 }
01605 }
01606 }
01607
01608 return elements;
01609 }
01610
01611 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01612 int Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01613 push(const Image<IMG_ELM_TYPE>* frame, const Image<MSK_ELM_TYPE>* mask)
01614 {
01615 unsigned int depth = setting_.ncut_tree_depth;
01616 unsigned int nEigVecs = setting_.number_of_eigenvectors;
01617 bool doMergeSegments = (setting_.merge_segments == 1);
01618 bool doSplitSegments = (setting_.connected_segments == 1);
01619
01620 if (depth == 0)
01621 throw ncutException("invalid depth, must be >= 1");
01622 if (nEigVecs < 2)
01623 throw ncutException(
01624 "invalid number of eigenvectors per ncut, must be >= 2" );
01625
01626
01627 unsigned int nElm_last = 0;
01628 unsigned int nElm_front = 0;
01629 double* assoc_last = NULL;
01630 if (dendogram_ != NULL) {
01631 nElm_last = elmProfile_.nElm();
01632 nElm_front = elmProfile_.maxIdx(0)+1;
01633 assoc_last = new double [nElm_last];
01634 for (unsigned int elm_last=0; elm_last<nElm_last; elm_last++)
01635 assoc_last[elm_last] = dendogram_->ncut()->assoc(elm_last);
01636 }
01637
01638
01639 if (elmProfile_.push(frame, mask) == 1)
01640 return 1;
01641
01642
01643 double* initVec = NULL;
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681 delete dendogram_;
01682 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > elmProfile_p(
01683 *dynamic_cast<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE>*>(&elmProfile_) );
01684 greedy_ptr<Ncut> ncut_p (NULL);
01685 dendogram_ = new NcutNode<IMG_ELM_TYPE, MSK_ELM_TYPE>(elmProfile_p,
01686 *elmProfile_.indices(),
01687 ncut_p, &setting_);
01688 if (dendogram_->calculate(depth,0,nEigVecs,1,initVec) == 1) {
01689 delete dendogram_;
01690 delete [] initVec;
01691 throw ncutException("Couldn't find any ncut");
01692 }
01693 delete [] initVec;
01694
01695
01696 std::vector<std::vector<MSK_ELM_TYPE> > elements;
01697 getInfo(dendogram_, elements);
01698
01699
01700 if (doMergeSegments) {
01701
01702 SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE> preMergeProfile(&setting_);
01703 preMergeProfile.push(&elmProfile_, elements);
01704 preMergeProfile.calcSim();
01705
01706 elements = merge(preMergeProfile);
01707 }
01708
01709
01710 if (doSplitSegments) {
01711
01712 SegmentProfile<IMG_ELM_TYPE, MSK_ELM_TYPE> preSplitProfile(&setting_);
01713 preSplitProfile.push(&elmProfile_, elements);
01714
01715 elements = split(preSplitProfile);
01716 }
01717
01718 segProfile_.push(&elmProfile_, elements);
01719
01720 #ifdef debugout
01721 for (unsigned int i=0; i<nSeg(); i++) {
01722 std::cout << "id " << id(i) << " area " << segProfile_.nPixels(i) << std::endl << std::flush;
01723 }
01724 #endif
01725
01726 return 0;
01727 }
01728
01729 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01730 Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01731 Segmentation(unsigned int wSize, Setting setting)
01732 : setting_(setting), dendogram_(NULL), elmProfile_(wSize, &setting_),
01733 segProfile_(&setting_)
01734 {
01735 return;
01736 }
01737
01738 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01739 Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01740 Segmentation(const Segmentation& clone)
01741 : setting_(clone.setting_), dendogram_(NULL), elmProfile_(3), segProfile_()
01742 {
01743 *this=clone;
01744 }
01745
01746 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01747 void Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::cleanup()
01748 {
01749 delete dendogram_;
01750 }
01751
01752 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01753 Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::~Segmentation()
01754 {
01755 cleanup();
01756 }
01757
01758 template <class IMG_ELM_TYPE, class MSK_ELM_TYPE>
01759 Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>&
01760 Segmentation<IMG_ELM_TYPE, MSK_ELM_TYPE>::
01761 operator=(const Segmentation& clone)
01762 {
01763
01764
01765
01766
01767 if (this == &clone) return *this;
01768
01769
01770 cleanup();
01771
01772
01773 this->elmProfile_ = clone.elmProfile_;
01774 this->segProfile_ = clone.segProfile_;
01775
01776
01777 greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > pclone(&this->elmProfile_);
01778 greedy_ptr<Ncut> nclone(NULL);
01779 if (clone.dendogram_ != NULL)
01780 this->dendogram_ =
01781 new NcutNode<IMG_ELM_TYPE,MSK_ELM_TYPE>(*clone.dendogram_,
01782 pclone, nclone);
01783 else
01784 this->dendogram_ = NULL;
01785
01786 this->setting_ = clone.setting_;
01787
01788 return *this;
01789 }
01790
01791 }
01792
01793 #endif // NCUT_CORE_H