ncutCore.h

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 // #include <time.h>
00012 
00013 // #include <arbsmat.h>
00014 // #include <arbgsym.h>
00015 // #include <arerror.h>
00016 
00017 // #include <assert.h>
00018 
00019 namespace ncut
00020 {
00021 
00031 class Ncut
00032 {
00033   public:
00034 
00035     // functions for read-only access to data members ------
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     // class-specific functions ------
00096 
00109     int calculate(unsigned int nEigVecs, double* initVec = NULL);
00110 
00111     // standard functions ------
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     // functions for read-only access to data members ------
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     // class-specific functions ------
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     // standard functions ------
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     // element info
00372     greedy_ptr<Profile<IMG_ELM_TYPE, MSK_ELM_TYPE> > profile_; 
00374     std::map<MSK_ELM_TYPE, unsigned int> indices_; 
00378     // ncut info
00379     greedy_ptr<Ncut> ncut_; 
00381     double eigVal_; 
00383     const double* eigVec_; 
00385     double splitPoint_; 
00386     unsigned int eigVecNum_; 
00389     // child nodes
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     // functions for read & write access to data members ------
00442 
00448     inline Setting* setting()
00449     {return &setting_;};
00450 
00451     // functions for read-only access to data members ------
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     // inactive segments
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 /* template implementations ------------------------------------------------- */
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   // free old data
00798   cleanup();
00799 
00800   // calculate the total association for each element
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   // save the association values in an array
00807   assoc_ = new double [dim_];
00808   for (unsigned int i=0; i<dim_; i++)
00809     assoc_[i] = assocMatrix.val(i,i);
00810 
00811   // calculate the second matrix for the generalized eigenvalue problem
00812   Matrix assocSubSimil(assocMatrix);
00813   assocSubSimil -= *similarity_;
00814 
00815   // convert matrices to band format for use with arpack
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     // construct arpack matrix objects
00826     ARbdSymMatrix<double> assocSubSimilARP(assocSubSimilN, assocSubSimilNsdiag,
00827                                            assocSubSimilBand.get());
00828     ARbdSymMatrix<double> assocARP(assocN, assocNsdiag, assocBand.get());
00829 
00830     // calculate eigenvalues- and vectors of the motion profile using arpack
00831     int ncv = ((2*nEigVecs)+1<dim_-1) ? (2*nEigVecs)+1 : dim_;
00832 //     int ncv = ((2*3)+1<dim_-1) ? (2*3)+1 : dim_;
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     // get eigenvalues and eigenvectors
00846     eigVal_ = new double[nEigVecs];
00847     eigVec_ = new double[nEigVecs*dim_];
00848     nConv_ = problem.EigenValVectors(eigVec_, eigVal_);
00849 //     for (int i=0; i<nEigVecs; i++) {
00850 //       std::cout << "EigVec[" << i << "] - (" << std::flush;
00851 //       for (int j=0; j<dim_; j++) {
00852 //         std::cout  << eigVec_[dim_*i+j] << "," << std::flush;
00853 //       }
00854 //       std::cout << ")" << std::endl << std::flush;
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   // check for self-assignment
00907   if (this == &clone) return *this;
00908 
00909   // free old memory
00910   cleanup();
00911 
00912   // copy ncut info from the clone
00913   this->similarity_ = clone.similarity_;
00914   this->dim_        = clone.dim_;
00915   this->nConv_      = clone.nConv_;
00916 
00917   // initialize the local structures
00918   assoc_  = new double[clone.dim_];
00919   eigVal_ = new double[clone.nConv_];
00920   eigVec_ = new double[clone.nConv_*clone.dim_];
00921 
00922   // copy the total association values
00923   for (unsigned int i=0; i<clone.dim_; i++) {
00924     this->assoc_[i] = clone.assoc_[i];
00925   }
00926 
00927   // copy Eigenvalues and Eigenvectors from the clone
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   // calculate data for this node ---------------
00974 
00975   // save the eigenvector number of this node
00976   eigVecNum_ = eigvn;
00977 
00978   // calculate new ncut if needed
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   // not enough converged eigenvectors, this node will be deleted
00993   if (ncut_->nConv() <= eigvn) return 1;
00994 
00995   // get eigenvalue and eigenvector
00996   eigVal_ = ncut_->eigVal(eigvn);
00997   eigVec_ = &( (ncut_->eigVec())[eigvn*ncut_->dim()] );
00998 
00999   // find a good splitting point
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; // maximum (worst) possible ncut value
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   // find a good splitting point v2 (to be tested)
01039 //   typename std::map<MSK_ELM_TYPE, unsigned int>::iterator elm;
01040 //   typename std::map<MSK_ELM_TYPE, unsigned int>::iterator elm2;
01041 //   double minVal, maxVal;
01042 //   minVal = maxVal = eigVec_[0];
01043 //   for (unsigned int i=0; i<ncut_->dim(); i++) {
01044 //     if (eigVec_[i] < minVal) minVal = eigVec_[i];
01045 //     if (eigVec_[i] > maxVal) maxVal = eigVec_[i];
01046 //   }
01047 //   double step = (maxVal - minVal) / (double) (nStep+1);
01048 //   if (step <= 0.0) return 1;
01049 // 
01050 //   double bestNcut = 2.0; // maximum (worst) possible ncut value
01051 //   double assoc_A, assoc_B, cut, ncutValue, last_split;
01052 //   assoc_A = assoc_B = cut = 0.0;
01053 //   last_split = minVal - (step/2.0);
01054 // 
01055 //   for (double split = minVal+(step/2.0); split<maxVal; split+=step) {
01056 //     if (last_split < minVal) {
01057 //     // first ncut value calculation goes through all elements
01058 //       for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01059 //         if (eigVec_[elm->second] < split)
01060 //           assoc_A += ncut_->assoc(elm->second);
01061 //         else
01062 //           assoc_B += ncut_->assoc(elm->second);
01063 // 
01064 //         for (elm2=elm, elm2++; elm2!=indices_.end(); elm2++) {
01065 //           if ((eigVec_[elm->second]  < split) !=
01066 //               (eigVec_[elm2->second] < split)) {
01067 //             cut += profile_->similarity()->val(elm->second,elm2->second);
01068 //           }
01069 //         }
01070 //       }
01071 //     }
01072 //     else {
01073 //     // following ncut value calculations only consider elements that changed
01074 //     // sides
01075 //       for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01076 //         if (eigVec_[elm->second]<split && eigVec_[elm->second]>last_split) {
01077 //         // elm changed sides from B to A
01078 //           assoc_A += ncut_->assoc(elm->second);
01079 //           assoc_B -= ncut_->assoc(elm->second);
01080 // 
01081 //           for (elm2=elm, elm2++; elm2!=indices_.end(); elm2++) {
01082 //             if (eigVec_[elm2->second] > split)
01083 //             // elm2 is still on the B side
01084 //               cut += profile_->similarity()->val(elm->second,elm2->second);
01085 //             else
01086 //             // both elements are on the A side
01087 //               cut -= profile_->similarity()->val(elm->second,elm2->second);
01088 //           }
01089 //         }
01090 //         else {
01091 //           for (elm2=elm, elm2++; elm2!=indices_.end(); elm2++) {
01092 //             if (eigVec_[elm2->second] < split &&
01093 //                 eigVec_[elm2->second] > last_split) {
01094 //             // elm2 changed sides from B to A
01095 //               if (eigVec_[elm->second] > split)
01096 //               // elm is still on the B side
01097 //                 cut += profile_->similarity()->val(elm->second,elm2->second);
01098 //               else
01099 //               // both elements are on the A side
01100 //                 cut -= profile_->similarity()->val(elm->second,elm2->second);
01101 //             }
01102 //           }
01103 //         }
01104 //       }
01105 //     }
01106 //     ncutValue = (cut/assoc_A) + (cut/assoc_B);
01107 //     if (ncutValue < bestNcut) {
01108 //       splitPoint_ = split;
01109 //       bestNcut = ncutValue;
01110 //     }
01111 //     last_split = split;
01112 //   }
01113 #ifdef debugout
01114   std::cout << "---splitPoint " << splitPoint_ << std::endl << std::flush;
01115 #endif
01116 
01117   // all done if this node has no child nodes
01118   if (depth >= maxDepth-1) return 0;
01119 
01120 
01121   // prepare data for child nodes ---------------------
01122 
01123   // calculate eigenvector numbers for both child nodes
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   // the element indices map for the left and right child nodes
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   // child nodes will only be more eigenvectors of the same ncut
01135   // (profiles stay the same)
01136 
01137     // construct the element indices lists for both child nodes
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     // this node doesn't do anything, it will be deleted
01146     if (left_indices.empty() || right_indices.empty()) return 1;
01147 
01148     // profiles stay the same
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     // construct the left child node (ncut_ is ignored if left_eigvn == 1)
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     // construct the right child node (ncut_ is ignored if right_eigvn == 1)
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   // both child node will be a new ncut
01169   // (profiles change)
01170 
01171     // gather element indices for left and right child nodes in a list
01172     std::list<unsigned int> left_elements;
01173     std::list<unsigned int> right_elements;
01174     for (elm=indices_.begin(); elm!=indices_.end(); elm++) {
01175 //       std::cout << "eigVec_[" << elm->second << "] - " << eigVec_[elm->second] << " - splitPoint " << splitPoint_ << std::endl << std::flush;
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     // this node doesn't do anything, it will be deleted
01183     if (left_elements.empty() || right_elements.empty()) return 1;
01184 
01185     // split the similarity matrix of this node twice (for both child nodes)
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     // construct the element indices map for both child nodes
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     // construct the element identifiers map for both child nodes
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     // construct the profiles for both child nodes
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     // construct the left child node (ncut_ is ignored if left_eigvn == 1)
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     // construct the right child node (ncut_ is ignored if right_eigvn == 1)
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; // everything done with greedy_ptr
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   // copy simple members
01323   this->indices_    = clone.indices_;
01324   this->eigVal_     = clone.eigVal_;
01325   this->splitPoint_ = clone.splitPoint_;
01326   this->eigVecNum_  = clone.eigVecNum_;
01327 
01328   // copy the profile:
01329   // if the clone owns it, copy the profile
01330   // if the clone does not own it and a profile pointer was handed down, use it
01331   // if no profile pointer was handed down, use the clone's profile pointer
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   // copy the ncut:
01340   // if the clone owns it, copy the ncut
01341   // if the clone does not own it and an ncut pointer was handed down, use it
01342   // if no ncut pointer was handed down, use the clone's ncut pointer
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   // copy the child nodes and make sure they use pointers
01359   // to this node's profile and ncut, not to the clone's.
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   // copy the setting
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   // get pointers to both child nodes
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   // this is a (half-)leaf node
01393 
01394     // get current segment numbers of both segments
01395     // and create vectors to hold the element id's
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     // save element id's for the two segments in separate vectors
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   // recurse to child nodes
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 // TODO: use already calculated profile similarity matrix instead of
01456 //       calculating a new one (to be tested)
01457 //   Matrix difference(*profile->similarity());
01458 
01459   // calculate the differences of each pair of segments
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   // merge segments until one of the stopping criteria is met
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     // search for the pair of segments with the smallest difference
01485     double minDif = 600.0; // dif shouldn't be bigger than 500
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     // if the smallest difference exceeds the merge threshold stop merging
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     // else merge the two segments with smallest difference
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     // update the difference matrix
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   // delete inactive segments
01532   // (segments that have been merged into another segment)
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     // is neighbor and not reached
01564       if (segProfile.idxMask(elm) == startSegIdx ) {
01565       // same segment
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   // split parts of segments that aren't touching
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     // mark elements as reached if they can be reached from this element
01590     // without crossing other segments
01591     connected(profile, unassigned.begin()->second, reached);
01592 
01593     // save those elements marked as reached as one segment and remove them
01594     // from the unassigned list
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   // save information needed to calculate the ncut starting vector
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   // push frame and mask to the element profile
01639   if (elmProfile_.push(frame, mask) == 1)
01640     return 1; // frame window is not filled yet
01641 
01642   // calculate the ncut starting vector (to be tested)
01643   double* initVec = NULL;
01644 //   if (dendogram_ != NULL) {
01645 // 
01646 //     // calculate the matrix of the associated non-generalized
01647 //     // eigenvalue problem
01648 //     Matrix a(elmProfile_.nElm());
01649 //     double inv_sqrt_assoc [elmProfile_.nElm()];
01650 //     for (unsigned int i=0; i<elmProfile_.nElm(); i++) {
01651 //       if (dendogram_->ncut()->assoc(i) == 0.0)
01652 //         throw ncutCritical("invalid total association value (0), "
01653 //                            "Segmentation left in inconsistent state");
01654 //       inv_sqrt_assoc[i] = 1.0/sqrt(dendogram_->ncut()->assoc(i));
01655 //     }
01656 //     for (unsigned int i=0; i<elmProfile_.nElm(); i++)
01657 //       for (unsigned int j=0; j<elmProfile_.nElm(); j++)
01658 //         a.val(i,j) = ( ((i==j)?dendogram_->ncut()->assoc(i):0.0) -
01659 //                       dendogram_->ncut()->similarity()->val(i,j)) *
01660 //                      inv_sqrt_assoc[j] * inv_sqrt_assoc[i];
01661 // 
01662 //     // calculate the eigenvector of the associated non-generalized
01663 //     // eigenvalue problem of the last iteration,
01664 //     // shifted to match the indices in this iteration
01665 //     initVec = new double [elmProfile_.nElm()];
01666 //     unsigned int elm, elm_last;
01667 //     for (elm_last=nElm_front, elm=0; elm_last<nElm_last; elm_last++, elm++)
01668 //       initVec[elm] = dendogram_->eigVec()[elm_last] *
01669 //                      sqrt(assoc_last[elm_last]);
01670 //     delete [] assoc_last;
01671 //     for (elm=elm;elm<elmProfile_.nElm(); elm++)
01672 //       initVec[elm] = 0.0;
01673 // 
01674 //     // transform the non-generalized eigenvector of the last iteration
01675 //     // with the matrix of the non-generalized eigenvalue problem of this
01676 //     // iteration
01677 //     a.transform(initVec);
01678 //   }
01679 
01680   // construct the ncut dendogram
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   // fetch segment information from nodes
01696   std::vector<std::vector<MSK_ELM_TYPE> > elements;
01697   getInfo(dendogram_, elements);
01698 
01699   // merge segments
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   // split parts of segments that aren't touching
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   // TODO: this does not create a completely independent copy, all setting
01764   //       pointers will still point to clone.setting_
01765 
01766   // check for self-assignment
01767   if (this == &clone) return *this;
01768 
01769   // free old memory
01770   cleanup();
01771 
01772   // copy the element profile and the segment profiles
01773   this->elmProfile_ = clone.elmProfile_;
01774   this->segProfile_ = clone.segProfile_;
01775 
01776   // copy the dendogram with pointers to this profile (not the clone's)
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 } // namespace ncut
01792 
01793 #endif // NCUT_CORE_H

Generated on Thu Jun 22 14:47:20 2006 for ncut.kdevelop by  doxygen 1.4.6