Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/McAulayQuatieri.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca>
00003 **
00004 ** This program is free software; you can redistribute it and/or modify
00005 ** it under the terms of the GNU General Public License as published by
00006 ** the Free Software Foundation; either version 2 of the License, or
00007 ** (at your option) any later version.
00008 **
00009 ** This program is distributed in the hope that it will be useful,
00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 ** GNU General Public License for more details.
00013 **
00014 ** You should have received a copy of the GNU General Public License
00015 ** along with this program; if not, write to the Free Software
00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 */
00018 
00019 #include "McAulayQuatieri.h"
00020 #include "../common_source.h"
00021 
00022 #include <marsyas/peakView.h>
00023 #include <marsyas/NumericLib.h>
00024 
00025 using std::ostringstream;
00026 using std::abs;
00027 
00028 using namespace Marsyas;
00029 
00030 McAulayQuatieri::McAulayQuatieri(mrs_string name):MarSystem("McAulayQuatieri", name)
00031 {
00032   addControls();
00033   nextGroup_ = 0;
00034 }
00035 
00036 McAulayQuatieri::McAulayQuatieri(const McAulayQuatieri& a) : MarSystem(a)
00037 {
00038   ctrl_reset_ = getctrl("mrs_bool/reset");
00039   ctrl_useGroups_ = getctrl("mrs_bool/useGroups");
00040   ctrl_useMemory_ = getctrl("mrs_bool/useMemory");
00041   ctrl_delta_ = getctrl("mrs_real/delta");
00042   ctrl_matchThres_ = getctrl("mrs_real/matchThres");
00043 
00044   nextGroup_ = a.nextGroup_;
00045 }
00046 
00047 McAulayQuatieri::~McAulayQuatieri()
00048 {
00049 }
00050 
00051 MarSystem*
00052 McAulayQuatieri::clone() const
00053 {
00054   return new McAulayQuatieri(*this);
00055 }
00056 
00057 void
00058 McAulayQuatieri::addControls()
00059 {
00060   addctrl("mrs_bool/reset", false);
00061   setctrlState("mrs_bool/reset", true);
00062 
00063   addctrl("mrs_bool/useMemory", false);
00064   //setctrlState("mrs_bool/useMemory", true);
00065 
00066   addctrl("mrs_bool/useGroups", false);
00067   //setctrlState("mrs_bool/useGroups", true);
00068 
00069   addctrl("mrs_real/delta", 0.5); //[TODO][!]
00070   addctrl("mrs_real/matchThres", 0.5);//[TODO][!]
00071 }
00072 
00073 void
00074 McAulayQuatieri::myUpdate(MarControlPtr sender)
00075 {
00076   MRSDIAG("McAulayQuatieri.cpp - McAulayQuatieri:myUpdate");
00077 
00078   MarSystem::myUpdate(sender);
00079 
00080   if(ctrl_reset_->to<mrs_bool>())
00081   {
00082     ctrl_reset_->setValue(false, NOUPDATE);
00083     memory_.stretch(0,0);
00084     nextGroup_ = 0;
00085   }
00086 }
00087 
00088 mrs_real
00089 McAulayQuatieri::peakTrack(realvec& vec, mrs_natural frame, mrs_natural grpOne, mrs_natural grpTwo)
00090 {
00091   mrs_real dist;
00092   mrs_natural candidate;
00093   mrs_natural lastMatched = -1;
00094   mrs_natural matchedTracks = 0;
00095 
00096   mrs_real delta = ctrl_delta_->to<mrs_real>();
00097 
00098   if(frame+1 >= vec.getCols())
00099   {
00100     MRSERR("McAulayQuatieri::peakTrack - frame index is bigger than the input vector!");
00101     return -1.0;
00102   }
00103 
00104   peakView tmpPeakView(vec);
00105 
00106   //get the trackID for any future track to be born (in STEP 3 - see below)
00107   mrs_natural nextTrack = tmpPeakView.getFrameNumPeaks(0, grpOne);
00108 
00109   //iterate over peaks in current frame
00110   for(mrs_natural n = 0; n < tmpPeakView.getFrameNumPeaks(frame, grpOne); ++n)
00111   {
00112     mrs_real lastdist = MAXREAL;
00113     candidate = -1;
00114 
00115     // STEP 1
00116     // find a candidate match on the next frame for each peak (i.e. track) in current frame
00117     for(mrs_natural m = lastMatched + 1; m < tmpPeakView.getFrameNumPeaks(frame+1, grpTwo); ++m)
00118     {
00119       //set track parameter of all peaks of next frame to -1 so we know later
00120       //which ones were not matched (=> BIRTH of new tracks)
00121       tmpPeakView(m, peakView::pkTrack, frame+1, grpTwo) = -1.0;
00122 
00123       dist = abs(tmpPeakView(n, peakView::pkFrequency, frame, grpOne) - tmpPeakView(m, peakView::pkFrequency, frame+1, grpTwo));
00124       if (dist < delta && dist < lastdist)
00125       {
00126         //found a candidate!
00127         lastdist  = dist;
00128         candidate = m;
00129       }
00130     }
00131 
00132     // STEP 2
00133     // must confirm candidate (if any)
00134     if(candidate >= 0) //check if a candidate was found
00135     {
00136       //confirm if this is not the last peak in current frame
00137       if(n < tmpPeakView.getFrameNumPeaks(frame, grpOne)-1)
00138       {
00139         //check the next remaining peak in current frame and see if it is a better match for the found candidate
00140         dist = abs(tmpPeakView(n+1, peakView::pkFrequency, frame, grpOne) - tmpPeakView(candidate, peakView::pkFrequency, frame+1, grpTwo));
00141         if(dist < lastdist)
00142         {
00143           // it is a better match! Check two additional conditions:
00144           // 1. an unmatched lower freq candidate should exist
00145           // 2. it is inside the frequency interval specified by delta
00146           if(candidate - 1 > lastMatched)
00147           {
00148             if(abs(tmpPeakView(n, peakView::pkFrequency, frame, grpOne) - tmpPeakView(candidate-1, peakView::pkFrequency, frame+1, grpTwo)) < delta)
00149             {
00150               //found a peak to continue the track -> confirm candidate!
00151               tmpPeakView(candidate-1, peakView::pkTrack, frame+1, grpTwo) = tmpPeakView(n, peakView::pkTrack, frame, grpOne);
00152               lastMatched = candidate-1;
00153               matchedTracks++;
00154             }
00155           }
00156         }
00157         else
00158         {
00159           //no better match than this one, so confirm candidate!
00160           tmpPeakView(candidate, peakView::pkTrack, frame+1, grpTwo) = tmpPeakView(n, peakView::pkTrack, frame, grpOne);
00161           lastMatched = candidate;
00162           matchedTracks++;
00163         }
00164       }
00165       else
00166       {
00167         //if this was the last peak in current frame, so inherently it was the best match.
00168         //Candidate is therefore automatically confirmed and can be propagated.
00169         tmpPeakView(candidate, peakView::pkTrack, frame+1, grpTwo) = tmpPeakView(n, peakView::pkTrack, frame, grpOne);
00170         lastMatched = candidate;
00171         matchedTracks++;
00172       }
00173     }
00174   } //end of loop on peaks of current frame
00175 
00176   // STEP 3
00177   // check for any unmatched peaks in the next frame and give BIRTH to new tracks!
00178   for(mrs_natural m = 0; m < tmpPeakView.getFrameNumPeaks(frame+1, grpTwo); ++m)
00179   {
00180     if(tmpPeakView(m, peakView::pkTrack, frame+1, grpTwo) == -1.0)
00181       tmpPeakView(m, peakView::pkTrack, frame+1, grpTwo) = nextTrack++; //BIRTH of new track
00182   }
00183 
00184   return matchedTracks;
00185 }
00186 
00187 void
00188 McAulayQuatieri::myProcess(realvec& in, realvec& out)
00189 {
00190   mrs_natural t,o,c;
00191   t=0;
00192   o=0;
00193   c=0;
00194 
00195   realvec* outPtr;
00196 
00197   out(o,t) = in(o,t);          //    ??????
00198 
00199   //if we want to use memory and we already have data from
00200   //past inputs (i.e. memory is not empty)...
00201   if(ctrl_useMemory_->to<mrs_bool>() && memory_.getSize() != 0)
00202   {
00203     //concatenate memory column vector with current input
00204     //so we can continue peak tracking from previous input
00205     tmp_.stretch(onObservations_, onSamples_+1);
00206     for(o = 0; o < onObservations_; ++o)
00207       tmp_(o, 0) = memory_(o);
00208     for(o = 0; o < onObservations_; ++o)
00209       for(c = 0; c < onSamples_; ++c)
00210         tmp_(o,c+1) = in(o,c);
00211     outPtr = &tmp_;
00212 
00213     //attempt matching of groups between the frame in memory
00214     //and the first frame from current input
00215     if(ctrl_useGroups_->to<mrs_bool>())
00216     {
00217       peakView inPV(in);
00218       mrs_realvec inFirstFrame;
00219       in.getCol(0, inFirstFrame);
00220       peakView inFirstFramePV(inFirstFrame);
00221       peakView memPV(memory_);
00222       peakView tmpPV(tmp_);
00223 
00224       //mrs_natural numInGroups = inPV.getNumGroups();
00225       mrs_natural numInFirstFrameGroups = inFirstFramePV.getNumGroups();
00226       mrs_natural numMemGroups = memPV.getNumGroups();
00227 
00228       //we must update the group numbers of the groups
00229       //in tmp realvec (i.e. [mem|current input])
00230       //so they do not clash with the previous ones
00231       if(nextGroup_ > 0)
00232         for(mrs_natural f=1; f < tmpPV.getNumFrames(); ++f)
00233           for(mrs_natural p = 0; p < tmpPV.getFrameNumPeaks(f); ++p)
00234             tmpPV(p, peakView::pkGroup, f) = tmpPV(p, peakView::pkGroup, f) + nextGroup_;
00235 
00236       // Try matching previous groups (in memory from last input)
00237       // with groups in current input
00238 
00239       // create a tmp copy of the frame in memory and the first frame
00240       // of current input, so we can do the group matching without
00241       // destroying the input values
00242       realvec frames2Match(inObservations_, 2);
00243 
00244       // calculate the matching score for all pairs of groups under matching
00245       realvec matchScores(numInFirstFrameGroups, numMemGroups);
00246       for(mrs_natural mg=0; mg < numMemGroups; ++mg)
00247       {
00248         for(mrs_natural ig = nextGroup_; ig < nextGroup_ + numInFirstFrameGroups; ++ig)
00249         {
00250           //since peakTrack(...) is destructible, we must reset frames2Match everytime... [!]
00251           for(o=0; o<inObservations_; ++o)
00252             for(c=0; c < 2; ++c)
00253               frames2Match(o, c) = tmp_(o, c);
00254 
00255           //use McAulay-Quatieri num of successful peak continuations as a score
00256           //for the group matching (may be replaced by some other metric in future)
00257           matchScores(ig-nextGroup_, mg) = peakTrack(frames2Match, 0, ig, mg);
00258         }
00259       }
00260 
00261       //Given the matchScores, try to find the optimal assignment
00262       //of the groups coming from previous input (stored in memory)
00263       //and the new groups in the current input
00264       //(using, for e.g. the hungarian method)
00265       realvec assignedGrp(numInFirstFrameGroups);
00266 
00267       //convert matchScores to costs
00268       mrs_real maxScore = matchScores.maxval();
00269       for(o=0; o < matchScores.getRows(); ++o)
00270         for(c=0; c < matchScores.getCols(); ++ c)
00271           matchScores(o,c) = maxScore - matchScores(o,c);
00272 
00273       NumericLib::hungarianAssignment(matchScores, assignedGrp); 
00274 
00275       // given the assignments, try to propagate the group IDs
00276       // to the groups in the current input
00277       mrs_natural ig;
00278       for(mrs_natural f=1; f < tmpPV.getNumFrames(); ++f)
00279       {
00280         for(mrs_natural p = 0; p < tmpPV.getFrameNumPeaks(f); ++p)
00281         {
00282           //get input group ID (converted to the range [0:...])
00283           ig = (mrs_natural)(tmpPV(p, peakView::pkGroup, f)) - nextGroup_;
00284 
00285           if(assignedGrp(ig) > -1) //a match was found for this group (ig)
00286           {
00287             //check if match is higher than the specified threshold
00288             if((maxScore - matchScores(ig, (mrs_natural)assignedGrp(ig))) / memPV.getFrameNumPeaks(0,(mrs_natural)assignedGrp(ig)) > ctrl_matchThres_->to<mrs_real>())
00289             {
00290               //match confirmed --> propagate group ID
00291               tmpPV(p, peakView::pkGroup, f) = assignedGrp(ig);
00292             }
00293             else //match below threshold --> set as new group
00294             {
00295               tmpPV(p, peakView::pkGroup, f) = nextGroup_;
00296               assignedGrp(ig) = nextGroup_;
00297               nextGroup_++;
00298             }
00299           }
00300           else //no match found for this group! --> set as new group
00301           {
00302             tmpPV(p, peakView::pkGroup, f) = nextGroup_;
00303             assignedGrp(ig) = nextGroup_;
00304             nextGroup_++;
00305           }
00306         }
00307       }
00308     }
00309   }
00310   else
00311   {
00312     //no need to concatenate memory information with
00313     //current input. Just do it inplace in the output realvec (avoid extra copy)!
00314     outPtr = &out;
00315   }
00316 
00317   peakView tmpPeakView(*outPtr);
00318 
00320   mrs_natural numGroups;
00321   mrs_natural g;
00322   if(ctrl_useGroups_->to<mrs_bool>())
00323   {
00324     numGroups = tmpPeakView.getNumGroups();
00325     g = 0;
00326   }
00327   else
00328   {
00329     numGroups = 0;
00330     g = -1;
00331   }
00332   //iterate over groups (if any or enabled)
00333   for(; g < numGroups; ++g)
00334   {
00335     //if no memory being used (or no memory stored yet), we must use peaks
00336     //in first frame to give birth to new tracks
00337     if(!ctrl_useMemory_->to<mrs_bool>() || memory_.getSize() == 0)
00338     {
00339       for(mrs_natural n = 0; n < tmpPeakView.getFrameNumPeaks(0, g); ++n)
00340         tmpPeakView(n, peakView::pkTrack, 0) = (mrs_real) n;
00341     }
00342 
00343     //iterate over input frames
00344     for(mrs_natural f=0; f < tmpPeakView.getNumFrames()-1; ++f)
00345       peakTrack(*outPtr, f, g, g);
00346   }
00347 
00348   //if using memory...
00349   if(ctrl_useMemory_->to<mrs_bool>())
00350   {
00351     if(memory_.getSize() != 0)
00352     {
00353       //if using a non-empty memory, we should now fill the trackID and GroupID parameters
00354       //computed above (and stored in the tmp realvec) into the actual output
00355       peakView outPeakView(out);
00356       for(mrs_natural f=0; f < outPeakView.getNumFrames(); ++f)
00357         for(mrs_natural p = 0; p < outPeakView.getFrameNumPeaks(f); ++p)
00358         {
00359           outPeakView(p, peakView::pkTrack, f) = tmpPeakView(p, peakView::pkTrack, f+1);
00360           outPeakView(p, peakView::pkGroup, f) = tmpPeakView(p, peakView::pkGroup, f+1);
00361         }
00362     }
00363 
00364     //store the last frame of current output for next time
00365     memory_.stretch(onObservations_, 1);
00366     for(o = 0; o < onObservations_; ++o)
00367       memory_(o, 0) = out(o, onSamples_-1);
00368   }
00369 }