cannam@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@0: 
cannam@0: /*
cannam@0:     Vamp feature extraction plugin using the MATCH audio alignment
cannam@0:     algorithm.
cannam@0: 
cannam@0:     Centre for Digital Music, Queen Mary, University of London.
Chris@236:     Copyright (c) 2007-2020 Simon Dixon, Chris Cannam, and Queen Mary
Chris@230:     University of London, Copyright (c) 2014-2015 Tido GmbH.
cannam@0:     
cannam@0:     This program is free software; you can redistribute it and/or
cannam@0:     modify it under the terms of the GNU General Public License as
cannam@0:     published by the Free Software Foundation; either version 2 of the
cannam@0:     License, or (at your option) any later version.  See the file
cannam@0:     COPYING included with this distribution for more information.
cannam@0: */
cannam@0: 
Chris@220: //#define DEBUG_FINDER 1
Chris@220: //#define PERFORM_ERROR_CHECKS 1
Chris@220: 
cannam@0: #include "Finder.h"
cannam@0: 
Chris@30: #include "Path.h"
Chris@30: 
Chris@30: #include <algorithm>
Chris@92: #include <iomanip>
Chris@30: 
Chris@72: using namespace std;
cannam@0: 
Chris@72: Finder::Finder(Matcher *pm)
cannam@0: {
Chris@72:     m_m = pm;
Chris@72:     m_duration1 = -1;
Chris@72:     m_duration2 = -1;
cannam@0: } // constructor
cannam@0: 
cannam@0: Finder::~Finder()
cannam@0: {
cannam@0: }
cannam@0: 
Chris@60: void
Chris@154: Finder::setMatcher(Matcher *pm)
Chris@154: {
Chris@155:     cerr << "Finder::setMatcher: finder " << this << ", matcher " << pm << endl;
Chris@154:     m_m = pm;
Chris@154: }
Chris@154: 
Chris@154: void
Chris@60: Finder::setDurations(int d1, int d2)
Chris@60: {
Chris@140: #ifdef DEBUG_FINDER
Chris@140:     cerr << "*** setDurations: " << d1 << ", " << d2 << endl;
Chris@140: #endif
Chris@72:     m_duration1 = d1;
Chris@72:     m_duration2 = d2;
Chris@60: }
Chris@60: 
Chris@154: bool
Chris@191: Finder::getBestRowCost(int row, int &bestCol, normpathcost_t &min)
Chris@154: {
Chris@172:     if (!m_m->isRowAvailable(row)) return false;
Chris@205:     pair<int, int> colRange = m_m->getColRangeForRow(row);
Chris@172:     if (colRange.first >= colRange.second) return false;
Chris@154:     for (int index = colRange.first; index < colRange.second; index++) {
Chris@191:         normpathcost_t tmp = m_m->getNormalisedPathCost(row, index);
Chris@154:         if (index == colRange.first || tmp < min) {
Chris@154:             min = tmp;
Chris@154:             bestCol = index;
Chris@154:         }
Chris@154:     }
Chris@154:     return true;
Chris@154: }    
Chris@154: 
Chris@154: bool
Chris@191: Finder::getBestColCost(int col, int &bestRow, normpathcost_t &min)
Chris@154: {
Chris@154:     if (!m_m->isColAvailable(col)) return false;
Chris@205:     pair<int, int> rowRange = m_m->getRowRangeForCol(col);
Chris@154:     if (rowRange.first >= rowRange.second) return false;
Chris@191:     bestRow = rowRange.first;
Chris@154:     for (int index = rowRange.first; index < rowRange.second; index++) {
Chris@191:         normpathcost_t tmp = m_m->getNormalisedPathCost(index, col);
Chris@191:         if (index == rowRange.first) {
Chris@191:             min = tmp;
Chris@191:         } else if (tmp < min) {
Chris@154:             min = tmp;
Chris@154:             bestRow = index;
Chris@154:         }
Chris@154:     }
Chris@154:     return true;
Chris@154: }    
Chris@154: 
Chris@147: void
Chris@147: Finder::getBestEdgeCost(int row, int col,
Chris@147:                         int &bestRow, int &bestCol,
Chris@191:                         normpathcost_t &min)
cannam@0: {
Chris@191:     min = m_m->getNormalisedPathCost(row, col);
Chris@72:     
Chris@147:     bestRow = row;
Chris@147:     bestCol = col;
Chris@72: 
Chris@205:     pair<int, int> rowRange = m_m->getRowRangeForCol(col);
Chris@72:     if (rowRange.second > row+1) {
Chris@72:         rowRange.second = row+1;	// don't cheat by looking at future :)
Chris@72:     }
Chris@72:     for (int index = rowRange.first; index < rowRange.second; index++) {
Chris@191:         normpathcost_t tmp = m_m->getNormalisedPathCost(index, col);
cannam@0:         if (tmp < min) {
cannam@0:             min = tmp;
cannam@0:             bestRow = index;
cannam@0:         }
cannam@0:     }
Chris@72: 
Chris@205:     pair<int, int> colRange = m_m->getColRangeForRow(row);
Chris@72:     if (colRange.second > col+1) {
Chris@72:         colRange.second = col+1;	// don't cheat by looking at future :)
Chris@72:     }
Chris@72:     for (int index = colRange.first; index < colRange.second; index++) {
Chris@191:         normpathcost_t tmp = m_m->getNormalisedPathCost(row, index);
cannam@0:         if (tmp < min) {
cannam@0:             min = tmp;
cannam@0:             bestCol = index;
cannam@0:             bestRow = row;
cannam@0:         }
cannam@0:     }
Chris@147: }
Chris@72: 
Chris@181: advance_t
Chris@171: Finder::getExpandDirection()
Chris@171: {
Chris@171:     return getExpandDirection(m_m->getFrameCount() - 1,
Chris@171:                               m_m->getOtherFrameCount() - 1);
Chris@171: }
Chris@171: 
Chris@181: advance_t
Chris@147: Finder::getExpandDirection(int row, int col)
Chris@147: {
Chris@147:     // To determine which direction to expand the search area in, we
Chris@147:     // look at the path costs along the leading edges of the search
Chris@147:     // area (the final row and column within the area). We find the
Chris@147:     // lowest path cost within the final row, and the lowest within
Chris@147:     // the final column, and we compare them. If the row is cheaper
Chris@147:     // then we expand by adding another row next to it; if the column
Chris@147:     // is cheaper then we expand by adding another column next to
Chris@147:     // it. (The overall lowest path cost across the row and column
Chris@147:     // represents the best alignment we have within the entire search
Chris@147:     // area given the data available and the assumption that the piece
Chris@147:     // is not ending yet.)
Chris@147: 
Chris@147:     int bestRow = row;
Chris@147:     int bestCol = col;
Chris@215:     normpathcost_t bestCost = INVALID_PATHCOST;
Chris@147: 
Chris@155: //    cerr << "Finder " << this << "::getExpandDirection: ";
Chris@155:     
Chris@147:     getBestEdgeCost(row, col, bestRow, bestCol, bestCost);
Chris@147: 
Chris@147: //    cerr << "at [" << row << "," << col << "] (cost " << m_m->getPathCost(row, col) << ") blocksize = " << m_m->getBlockSize() << " best is [" << bestRow << "," << bestCol << "] (cost " << bestCost << ")" << endl;
Chris@135:     
Chris@45:     if (bestRow == row) {
Chris@45:         if (bestCol == col) {
Chris@181:             return AdvanceBoth;
Chris@45:         } else {
Chris@181:             return AdvanceThis;
Chris@45:         }
Chris@45:     } else if (bestCol == col) {
Chris@181:         return AdvanceOther;
Chris@45:     } else {
Chris@181:         return AdvanceNone;
Chris@45:     }
Chris@73: }
cannam@0: 
cannam@0: void
cannam@0: Finder::recalculatePathCostMatrix(int r1, int c1, int r2, int c2) 
cannam@0: {
Chris@72:     int prevRowStart = 0, prevRowStop = 0;
Chris@72: 
Chris@72:     for (int r = r1; r <= r2; r++) {
Chris@72: 
Chris@205:         pair<int, int> colRange = m_m->getColRangeForRow(r);
Chris@72: 
Chris@72:         int rowStart = max(c1, colRange.first);
Chris@72:         int rowStop = min(c2 + 1, colRange.second);
Chris@72:         
Chris@72:         for (int c = rowStart; c < rowStop; c++) {
Chris@72: 
Chris@181:             advance_t dir = AdvanceNone;
Chris@188:             pathcost_t straightIncrement = m_m->getDistance(r, c);
Chris@188:             pathcost_t diagIncrement = pathcost_t(straightIncrement *
Chris@188:                                                   m_m->getDiagonalWeight());
Chris@72: 
Chris@72:             if (r > r1) {	// not first row
Chris@215:                 pathcost_t min = INVALID_PATHCOST;
Chris@72:                 if ((c > prevRowStart) && (c <= prevRowStop)) {
Chris@72:                     // diagonal from (r-1,c-1)
Chris@188:                     min = m_m->getPathCost(r-1, c-1) + diagIncrement;
Chris@181:                     dir = AdvanceBoth;
Chris@72:                 }
Chris@72:                 if ((c >= prevRowStart) && (c < prevRowStop)) {
Chris@72:                     // vertical from (r-1,c)
Chris@188:                     pathcost_t cost = m_m->getPathCost(r-1, c) + straightIncrement;
Chris@215:                     if ((min == INVALID_PATHCOST) || (cost < min)) {
Chris@72:                         min = cost;
Chris@181:                         dir = AdvanceThis;
Chris@72:                     }
Chris@72:                 }
Chris@72:                 if (c > rowStart) {
Chris@72:                     // horizontal from (r,c-1)
Chris@188:                     pathcost_t cost = m_m->getPathCost(r, c-1) + straightIncrement;
Chris@215:                     if ((min == INVALID_PATHCOST) || (cost < min)) {
Chris@72:                         min = cost;
Chris@181:                         dir = AdvanceOther;
Chris@72:                     }
Chris@72:                 }
Chris@72:                 
Chris@72:                 m_m->setPathCost(r, c, dir, min);
Chris@72: 
Chris@72:             } else if (c > rowStart) {	// first row
Chris@72:                 // horizontal from (r,c-1)
Chris@181:                 m_m->setPathCost(r, c, AdvanceOther,
Chris@188:                                  m_m->getPathCost(r, c-1) + straightIncrement);
Chris@72:             }
Chris@72:         }
Chris@72: 
Chris@72:         prevRowStart = rowStart;
Chris@72:         prevRowStop = rowStop;
cannam@0:     }
Chris@72: } 
Chris@30: 
Chris@82: #ifdef PERFORM_ERROR_CHECKS
Chris@81: Finder::ErrorPosition
Chris@81: Finder::checkPathCostMatrix() 
Chris@81: {
Chris@81:     ErrorPosition err;
Chris@81: 
Chris@81:     int r1 = 0;
Chris@81:     int c1 = 0;
Chris@81:     int r2 = m_m->getFrameCount() - 1;
Chris@81:     int c2 = m_m->getOtherFrameCount() - 1;
Chris@81: 
Chris@81:     if (r2 < r1 || c2 < c1) {
Chris@81:         return err;
Chris@81:     }
Chris@81: 
Chris@81:     int prevRowStart = 0, prevRowStop = 0;
Chris@81: 
Chris@81:     for (int r = r1; r <= r2; r++) {
Chris@81: 
Chris@205:         pair<int, int> colRange = m_m->getColRangeForRow(r);
Chris@81: 
Chris@81:         int rowStart = max(c1, colRange.first);
Chris@81:         int rowStop = min(c2 + 1, colRange.second);
Chris@81:         
Chris@81:         for (int c = rowStart; c < rowStop; c++) {
Chris@81: 
Chris@188:             advance_t dir = AdvanceNone;
Chris@212:             pathcost_t updateTo = INVALID_PATHCOST;
Chris@188:             distance_t distance = m_m->getDistance(r, c);
Chris@188:             pathcost_t straightIncrement = distance;
Chris@188:             pathcost_t diagIncrement = pathcost_t(distance * m_m->getDiagonalWeight());
Chris@188:             err.distance = distance;
Chris@81: 
Chris@95:             if (r > r1) { // not first row
Chris@215:                 pathcost_t min = INVALID_PATHCOST;
Chris@81:                 if ((c > prevRowStart) && (c <= prevRowStop)) {
Chris@81:                     // diagonal from (r-1,c-1)
Chris@188:                     min = m_m->getPathCost(r-1, c-1) + diagIncrement;
Chris@81:                     err.prevCost = m_m->getPathCost(r-1, c-1);
Chris@181:                     dir = AdvanceBoth;
Chris@81:                 }
Chris@81:                 if ((c >= prevRowStart) && (c < prevRowStop)) {
Chris@81:                     // vertical from (r-1,c)
Chris@188:                     pathcost_t cost = m_m->getPathCost(r-1, c) + straightIncrement;
Chris@215:                     if ((min == INVALID_PATHCOST) || (cost < min)) {
Chris@81:                         min = cost;
Chris@81:                         err.prevCost = m_m->getPathCost(r-1, c);
Chris@181:                         dir = AdvanceThis;
Chris@81:                     }
Chris@81:                 }
Chris@81:                 if (c > rowStart) {
Chris@81:                     // horizontal from (r,c-1)
Chris@188:                     pathcost_t cost = m_m->getPathCost(r, c-1) + straightIncrement;
Chris@215:                     if ((min == INVALID_PATHCOST) || (cost < min)) {
Chris@81:                         min = cost;
Chris@81:                         err.prevCost = m_m->getPathCost(r, c-1);
Chris@181:                         dir = AdvanceOther;
Chris@81:                     }
Chris@81:                 }
Chris@81: 
Chris@81:                 updateTo = min;
Chris@81: 
Chris@82:             } else { // first row
Chris@82: 
Chris@82:                 if (c > rowStart) {
Chris@82:                     // horizontal from (r,c-1)
Chris@188:                     updateTo = m_m->getPathCost(r, c-1) + straightIncrement;
Chris@83:                     err.prevCost = m_m->getPathCost(r, c-1);
Chris@181:                     dir = AdvanceOther;
Chris@82:                 }
Chris@81:             }
Chris@81: 
Chris@181:             if (dir != AdvanceNone) {
Chris@86:                 if (m_m->getAdvance(r, c) != dir) {
Chris@86:                     err.type = ErrorPosition::WrongAdvance;
Chris@86:                     err.r = r;
Chris@86:                     err.c = c;
Chris@86:                     err.costWas = m_m->getPathCost(r, c);
Chris@86:                     err.costShouldBe = updateTo;
Chris@86:                     err.advanceWas = m_m->getAdvance(r, c);
Chris@86:                     err.advanceShouldBe = dir;
Chris@86:                     return err;
Chris@86:                 }
Chris@84:                 if (m_m->getPathCost(r, c) != updateTo) {
Chris@84:                     err.type = ErrorPosition::WrongCost;
Chris@84:                     err.r = r;
Chris@84:                     err.c = c;
Chris@84:                     err.costWas = m_m->getPathCost(r, c);
Chris@84:                     err.costShouldBe = updateTo;
Chris@84:                     err.advanceWas = m_m->getAdvance(r, c);
Chris@84:                     err.advanceShouldBe = dir;
Chris@82:                     return err;
Chris@82:                 }
Chris@82:             } else {
Chris@82:                 // AdvanceNone should occur only at r = r1, c = c1
Chris@82:                 if (r != r1 || c != c1) {
Chris@82:                     err.type = ErrorPosition::NoAdvance;
Chris@82:                     err.r = r;
Chris@82:                     err.c = c;
Chris@82:                     err.costWas = m_m->getPathCost(r, c);
Chris@82:                     err.costShouldBe = updateTo;
Chris@84:                     err.advanceWas = m_m->getAdvance(r, c);
Chris@84:                     err.advanceShouldBe = dir;
Chris@82:                     return err;
Chris@82:                 }
Chris@81:             }
Chris@81:         }
Chris@81: 
Chris@81:         prevRowStart = rowStart;
Chris@81:         prevRowStop = rowStop;
Chris@81:     }
Chris@81: 
Chris@81:     return err;
Chris@82: }
Chris@81: 
Chris@92: void
Chris@92: Finder::checkAndReport()
Chris@30: {
Chris@92:     cerr << "Finder: Checking path-cost matrix..." << endl;
Chris@82:     ErrorPosition err = checkPathCostMatrix();
Chris@92:     if (err.type == ErrorPosition::NoError) {
Chris@92:         cerr << "No errors found" << endl;
Chris@92:     } else {
Chris@82:         cerr << "\nWARNING: Checking path-cost matrix returned mismatch:" << endl;
Chris@92:         cerr << "Type: " << err.type << ": ";
Chris@92:         switch (err.type) {
Chris@92:         case ErrorPosition::NoError: break;
Chris@92:         case ErrorPosition::WrongCost: cerr << "WrongCost"; break;
Chris@92:         case ErrorPosition::WrongAdvance: cerr << "WrongAdvance"; break;
Chris@92:         case ErrorPosition::NoAdvance: cerr << "NoAdvance"; break;
Chris@92:         }
Chris@92:         cerr << endl;
Chris@84:         cerr << "At row " << err.r << ", column " << err.c
Chris@84:              << "\nShould be advancing "
Chris@84:              << Matcher::advanceToString(err.advanceShouldBe)
Chris@84:              << ", advance in matrix is "
Chris@84:              << Matcher::advanceToString(err.advanceWas)
Chris@83:              << "\nPrev cost " << err.prevCost
Chris@191:              << " plus distance " << distance_print_t(err.distance)
Chris@191:              << " [perhaps diagonalised] gives "
Chris@84:              << err.costShouldBe << ", matrix contains " << err.costWas
Chris@83:              << endl;
Chris@83:         cerr << "Note: diagonal weight = " << m_m->getDiagonalWeight() << endl;
Chris@83:         cerr << endl;
Chris@92: 
Chris@95:         int w(4);
Chris@95:         int ww(15);
Chris@92: 
Chris@92:         cerr << "Distance matrix leading up to this point:" << endl;
Chris@95:         cerr << setprecision(12) << setw(w) << "";
Chris@92:         for (int i = -4; i <= 0; ++i) {
Chris@95:             cerr << setw(ww) << i;
Chris@92:         }
Chris@92:         cerr << endl;
Chris@92:         for (int j = -4; j <= 0; ++j) {
Chris@92:             cerr << setw(w) << j;
Chris@92:             for (int i = -4; i <= 0; ++i) {
Chris@191:                 cerr << setw(ww)
Chris@191:                      << distance_print_t(m_m->getDistance(err.r + j, err.c + i));
Chris@92:             }
Chris@92:             cerr << endl;
Chris@92:         }
Chris@92:         cerr << endl;
Chris@92: 
Chris@92:         cerr << "Cost matrix leading up to this point:" << endl;
Chris@92:         cerr << setw(w) << "";
Chris@92:         for (int i = -4; i <= 0; ++i) {
Chris@95:             cerr << setw(ww) << i;
Chris@92:         }
Chris@92:         cerr << endl;
Chris@92:         for (int j = -4; j <= 0; ++j) {
Chris@92:             cerr << setw(w) << j;
Chris@92:             for (int i = -4; i <= 0; ++i) {
Chris@95:                 cerr << setw(ww) << m_m->getPathCost(err.r + j, err.c + i);
Chris@92:             }
Chris@92:             cerr << endl;
Chris@92:         }
Chris@92:         cerr << endl;
Chris@82:     }
Chris@92: }
Chris@92: #endif
Chris@92: 
Chris@182: pathcost_t
Chris@163: Finder::getOverallCost()
Chris@163: {
Chris@163:     int ex = m_m->getOtherFrameCount() - 1;
Chris@163:     int ey = m_m->getFrameCount() - 1;
Chris@163: 
Chris@163:     if (ex < 0 || ey < 0) {
Chris@163:         return 0;
Chris@163:     }
Chris@163: 
Chris@165:     return m_m->getPathCost(ey, ex);
Chris@163: }
Chris@163: 
Chris@92: int
Chris@92: Finder::retrievePath(bool smooth, vector<int> &pathx, vector<int> &pathy)
Chris@92: {
Chris@92:     pathx.clear();
Chris@92:     pathy.clear();
Chris@92: 
Chris@92: #ifdef PERFORM_ERROR_CHECKS
Chris@92:     checkAndReport();
Chris@82: #endif
Chris@82: 
Chris@72:     int ex = m_m->getOtherFrameCount() - 1;
Chris@72:     int ey = m_m->getFrameCount() - 1;
Chris@69: 
Chris@69:     if (ex < 0 || ey < 0) {
Chris@69:         return 0;
Chris@69:     }
Chris@66:     
Chris@66:     int x = ex;
Chris@66:     int y = ey;
Chris@30: 
Chris@140: #ifdef DEBUG_FINDER
Chris@140:     cerr << "*** retrievePath: smooth = " << smooth << endl;
Chris@140:     cerr << "*** retrievePath: before: x = " << x << ", y = " << y << endl;
Chris@140: #endif
Chris@140: 
Chris@72:     if (m_duration2 > 0 && m_duration2 < m_m->getOtherFrameCount()) {
Chris@72:         x = m_duration2 - 1;
Chris@60:     }
Chris@72:     if (m_duration1 > 0 && m_duration1 < m_m->getFrameCount()) {
Chris@72:         y = m_duration1 - 1;
Chris@60:     }
Chris@60: 
Chris@72:     if (!m_m->isAvailable(y, x)) {
Chris@66:         // Path did not pass through the expected end point --
Chris@66:         // probably means the pieces are substantially different in
Chris@66:         // the later bits. Reset the expected end point to the end of
Chris@66:         // both files including any trailing silence.
Chris@66:         cerr << "NOTE: Path did not pass through expected end point, inputs are probably significantly different" << endl;
Chris@66:         x = ex;
Chris@66:         y = ey;
Chris@66:     }
Chris@66: 
Chris@55:     recalculatePathCostMatrix(0, 0, y, x);
Chris@55: 
Chris@140: #ifdef DEBUG_FINDER
Chris@140:     cerr << "*** retrievePath: start: x = " << x << ", y = " << y << endl;
Chris@140: #endif
Chris@66:     
Chris@72:     while (m_m->isAvailable(y, x) && (x > 0 || y > 0)) {
Chris@30: 
Chris@33: //        cerr << "x = " << x << ", y = " << y;
Chris@33:         
Chris@30:         pathx.push_back(x);
Chris@30:         pathy.push_back(y);
Chris@30: 
Chris@72:         switch (m_m->getAdvance(y, x)) {
Chris@181:         case AdvanceThis:
Chris@70: //            cerr << ", going down (dist = " << getDistance() << ")" << endl;
Chris@33:             y--;
Chris@33:             break;
Chris@181:         case AdvanceOther:
Chris@70: //            cerr << ", going left (dist = " << getDistance() << ")" << endl;
Chris@33:             x--;
Chris@33:             break;
Chris@181:         case AdvanceBoth:
Chris@70: //            cerr << ", going diag (dist = " << getDistance() << ")" << endl;
Chris@33:             x--;
Chris@33:             y--;
Chris@33:             break;
Chris@181:         case AdvanceNone: // this would indicate a bug, but we wouldn't want to hang
Chris@69:             cerr << "WARNING: Neither matcher advanced in path backtrack at (" << x << "," << y << ")" << endl;
Chris@33:             if (x > y) {
Chris@33:                 x--;
Chris@33:             } else {
Chris@33:                 y--;
Chris@33:             }
Chris@33:             break;
Chris@30:         }
Chris@30:     }
Chris@30: 
Chris@72:     if (x > 0 || y > 0) {
Chris@72:         cerr << "WARNING: Ran out of available path at (" << y << "," << x
Chris@72:              << ")!" << endl;
Chris@72:     }
Chris@72:     
Chris@72:     reverse(pathx.begin(), pathx.end());
Chris@72:     reverse(pathy.begin(), pathy.end());
Chris@30: 
Chris@31:     if (smooth) {
Chris@180:         int smoothedLen = Path().smooth
Chris@180:             (pathx, pathy, static_cast<int>(pathx.size()));
Chris@31:         return smoothedLen;
Chris@31:     } else {
Chris@180:         return static_cast<int>(pathx.size());
Chris@31:     }
Chris@30: }
Chris@30: