    //=============================================================================
    //  MusE Reader
    //  Music Score Reader
    //  $Id$
    //
    //  Copyright (C) 2010-2011 Werner Schweer
    //
    //  This program is free software; you can redistribute it and/or modify
    //  it under the terms of the GNU General Public License version 2.
    //
    //  This program is distributed in the hope that it will be useful,
    //  but WITHOUT ANY WARRANTY; without even the implied warranty of
    //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    //  GNU General Public License for more details.
    //
    //  You should have received a copy of the GNU General Public License
    //  along with this program; if not, write to the Free Software
    //  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
    //=============================================================================

    #include "omrpage.h"
    #include "image.h"
    #include "utils.h"
    #include "omr.h"
    #ifdef OCR
    #include "ocr.h"
    #endif
    #include "libmscore/score.h"
    #include "libmscore/text.h"
    #include "libmscore/measurebase.h"
    #include "libmscore/box.h"
    #include "libmscore/sym.h"
    #include "libmscore/note.h"
    #include "pattern.h"

    namespace Ms {
        //static const double noteTH = 1.0;
        static const double timesigTH = 0.7;
        //static const double clefTH = 0.7;
        static const double keysigTH = 0.8;

        struct Hv {
            int x;
            int val;
            Hv(int a, int b) : x(a), val(b) {}
            bool operator< (const Hv& a) const { return a.val < val; }
        };

    struct Peak {
          int x;
          double val;
          int sym;

          bool operator<(const Peak& p) const { return p.val < val; }
          Peak(int _x, double _val) : x(_x), val(_val) {}
          Peak(int _x, double _val, int s) : x(_x), val(_val), sym(s) {}
          };

    //---------------------------------------------------------
    //   Lv
    //    line + value
    //---------------------------------------------------------

    struct Lv {
          int line;
          double val;
          Lv(int a, double b) : line(a), val(b) {}
          bool operator< (const Lv& a) const { return a.val < val; }
          };

    //---------------------------------------------------------
    //   OmrPage
    //---------------------------------------------------------

    OmrPage::OmrPage(Omr* parent)
          {
          _omr = parent;
          cropL = cropR = cropT = cropB = 0;
          }

    //---------------------------------------------------------
    //   dot
    //---------------------------------------------------------

    bool OmrPage::dot(int x, int y) const
          {
          const uint* p = scanLine(y) + (x / 32);
          return (*p) & (0x1 << (x % 32));
          }

    //---------------------------------------------------------
    //   isBlack
    //---------------------------------------------------------

    bool OmrPage::isBlack(int x, int y) const
          {
          QRgb c = _image.pixel(x,y);
          return (qGray(c) < 100);
          }

    //---------------------------------------------------------
    //   read
    //---------------------------------------------------------

    void OmrPage::read()
          {
          //removeBorder();
          crop();
          slice();
          deSkew();
          crop();
          slice();
          getStaffLines();
          getRatio();
          }

    struct SysState {
          int status;// 0 for start_staff, 1 for end_staff
          int index;//staff index
          };

    struct BAR_STATE {
          int x;
          int status;//0 represents white space, 1 represents bar
          };

    //---------------------------------------------------------
    //   searchBarLines
    //---------------------------------------------------------
    float OmrPage::searchBarLines(int start_staff, int end_staff)
          {
          OmrStaff& r1 = staves[start_staff];
          OmrStaff& r2 = staves[end_staff];

          int x1 = r1.x();
          int x2 = x1 + r1.width();
          int y1 = r1.y();
          int y2 = r2.y() + r2.height();

          int th = 0;
          for (int i = start_staff; i <= end_staff; ++i) {
                th += staves[i].height() / 2;
                }

          int vpw = x2 - x1;
          float vp[vpw];
          memset(vp, 0, sizeof(float) * vpw);

          //using note constraints
          //searchNotes();

          int note_constraints[x2 - x1];
          for (OmrNote* n : r1.notes()) {
                for(int x = n->x(); x <= n->x() + n->width(); ++x)
                      note_constraints[x - x1] = 1;
                }
          for (OmrNote* n : r2.notes()) {
                for(int x = n->x(); x <= n->x() + n->width(); ++x)
                      note_constraints[x - x1] = 1;
                }

          //
          // compute vertical projections
          //
          for (int x = x1; x < x2; ++x) {
                int dots = 0;
                for (int y = y1; y < y2; ++y) {
                      if (this->dot(x, y))
                            ++dots;
                      }
                if (!note_constraints[x - x1])
                      vp[x - x1] = dots - th;
                else
                      vp[x - x1] = -HUGE_VAL;
                }

          float scores[x2 - x1 + 1][2];
          BAR_STATE pred[x2 - x1 + 1][2];
          BAR_STATE bs;

          //initialization
          bs.x = -1; bs.status = -1;
          for (int x = x1; x <= x2; ++x) {
                int i = x - x1;
                for (int status = 0; status <= 1; ++status) {
                      scores[i][status] = -HUGE_VAL;
                      pred[i][status] = bs;
                      }
                }
          scores[0][0] = 0;

          //forward pass
          for (int x = x1; x <= x2; ++x) {
                int i = x - x1;
                for (int cur = 0; cur <= 1; ++cur) {
                      //current state
                      if (scores[i][cur] < -1000)
                            continue;

                      if (cur) {
                            scores[i][cur] += vp[i];
                            int next = 0;
                            int step = 3 * _spatium;//constraints between adjacent barlines
                            if (step + i <= x2 - x1) {
                                  if (scores[step + i][next] < scores[i][cur]) {
                                        scores[step + i][next] = scores[i][cur];
                                        bs.x = x; bs.status = cur;
                                        pred[step + i][next] = bs;
                                        }
                                  }
                            }
                      else {
                            for (int next = 0; next <= 1; ++next) {
                                  int step = 1;
                                  if (step + i <= x2 - x1) {
                                        if (scores[step + i][next] < scores[i][cur]) {
                                              scores[step + i][next] = scores[i][cur];
                                              bs.x = x; bs.status = cur;
                                              pred[step + i][next] = bs;
                                              }
                                        }
                                  }
                            }
                      }
                }
          return(scores[x2][0]);
          }

    //---------------------------------------------------------
    //   identifySystems
    //---------------------------------------------------------
    void OmrPage::identifySystems()
          {
          int numStaves = staves.size();
          if(numStaves == 0) return;

          //
          //memory allocation
          //
          float **temp_scores = new float*[numStaves];
          for (int i = 0; i < numStaves; i++)
                temp_scores[i] = new float[numStaves];

          int **hashed = new int*[numStaves];
          for (int i = 0; i < numStaves; i++)
                hashed[i] = new int[numStaves];

          float **scores = new float*[numStaves];
          for (int i = 0; i < numStaves; i++)
                scores[i] = new float[2];
          SysState **pred = new SysState*[numStaves];
          for (int i = 0; i < numStaves; i++)
                pred[i] = new SysState[2];

          //
          //initialization
          //
          for (int i = 0; i < numStaves; i++) {
                for (int j = 0; j < numStaves; j++) {
                      hashed[i][j] = 0;
                      }
                }

          SysState ss;
          ss.index = -1; ss.status = -1;

          for (int i = 0; i < numStaves; ++i) {
                for (int j = 0; j < 2; j++) {
                      scores[i][j] = -HUGE_VAL;
                      pred[i][j] = ss;
                      }
                }
          scores[0][0] = 0;

          //
          //identify solid note heads
          //
          int **note_labels = new int*[numStaves];
          for (int i = 0; i < numStaves; i++)
                note_labels[i] = new int[width()];
          for (int i = 0; i < numStaves; i++){
                for (int j = 0; j < width(); j++){
                      note_labels[i][j] = 0;
                      }
                }

          //
          // search notes for each system
          //
    //      int note_ran = spatium()/4.0;
    //      for (int i = 0; i < numStaves; i++){
    //            OmrSystem omrSystem(this);
    //            omrSystem.staves().append(staves[i]);
    //            omrSystem.searchNotes(note_labels[i], note_ran);
    //            }

          //
          // System Identification
          //
          int status;
          float cur_score;
          int cur_staff,next_staff;
          for (cur_staff = 0; cur_staff < numStaves; ++cur_staff) {
                status = 0;//start_staff
                for (next_staff = cur_staff; next_staff < numStaves; ++next_staff) { //connects to end_staff
                      if (hashed[cur_staff][next_staff]) {
                            cur_score = temp_scores[cur_staff][next_staff];
                            }
                      else {
                            //evaluate staff segment [c,n], dp here
                            OmrSystem omrSystem(this);
                            for (int i = cur_staff; i <= next_staff; ++i) {
                                  omrSystem.staves().append(staves[i]);
                                  }

                            cur_score = omrSystem.searchBarLinesvar(next_staff - cur_staff + 1, note_labels + cur_staff);
                            temp_scores[cur_staff][next_staff] = cur_score;
                            hashed[cur_staff][next_staff] = 1;
                            }

                      //forward pass
                      if (scores[cur_staff][status] + cur_score > scores[next_staff][1 - status]) {
                            scores[next_staff][1 - status] = scores[cur_staff][status] + cur_score;
                            ss.status = status; ss.index = cur_staff;
                            pred[next_staff][1 - status] = ss;
                            }
                      }

                status = 1;//end_staff
                next_staff = cur_staff + 1;
                if (next_staff < numStaves) {
                      //forward pass
                      if (scores[cur_staff][status] > scores[next_staff][1 - status]) {
                            scores[next_staff][1 - status] = scores[cur_staff][status];
                            ss.status = status; ss.index = cur_staff;
                            pred[next_staff][1 - status] = ss;
                            }
                      }
                }

          //backtrack
          status = 1;//end_staff
          cur_staff = numStaves - 1;//last staff index
          while (cur_staff > 0 || status != 0) {
                ss = pred[cur_staff][status];

                if (!ss.status) {
                      //add system here
                      OmrSystem omrSystem(this);
                      for (int i = ss.index; i <= cur_staff; ++i) {
                            omrSystem.staves().append(staves[i]);
                            }
                      omrSystem.searchBarLinesvar(cur_staff - ss.index + 1, note_labels + ss.index);
                      _systems.append(omrSystem);
                      }
                cur_staff = ss.index;
                status = ss.status;
                }

          int systems = _systems.size();
          for (int i = 0; i < systems; ++i) {
                OmrSystem* system = &_systems[i];
                int n = system->barLines.size();
                for (int k = 0; k < n - 1; ++k) {
                      const QLine& l1 = system->barLines[k];
                      const QLine& l2 = system->barLines[k + 1];
                      OmrMeasure m(l1.x1(), l2.x1());
                      system->measures().append(m);
                      }
                }

          //delete allocated space
          for (int i = 0; i < numStaves; i++) {
                delete[] scores[i];
                delete[] pred[i];
                delete[] temp_scores[i];
                delete[] hashed[i];
                delete[] note_labels[i];
                }

          delete[] note_labels;
          delete[] scores;
          delete[] pred;
          delete[] temp_scores;
          delete[] hashed;
          }

    //---------------------------------------------------------
    //   readBarLines
    //---------------------------------------------------------

    void OmrPage::readBarLines()
          {
          //QFuture<void> bl = QtConcurrent::run(_systems, &OmrSystem::searchSysBarLines());
          //bl.waitForFinished();

          //int numStaves    = staves.size();
          int systems = _systems.size();

          for (int i = 0; i < systems; ++i) {
                OmrSystem* system = &_systems[i];
                int n = system->barLines.size();
                for (int k = 0; k < n - 1; ++k) {
                      const QLine& l1 = system->barLines[k];
                      const QLine& l2 = system->barLines[k + 1];
                      OmrMeasure m(l1.x1(), l2.x1());
                      for (int ss = 0; ss < system->staves().size(); ++ss) {
                            OmrStaff& staff = system->staves()[ss];
                            QList<OmrChord> chords;
                            int nx = 0;
                            SymId nsym = SymId::noSym;
                            OmrChord chord;
                            foreach(OmrNote* n, staff.notes()) {
                                  int x = n->x();
                                  if (x >= m.x2())
                                        break;
                                  if (x >= m.x1() && x < m.x2()) {
                                        if (qAbs(x - nx) > int(_spatium / 2) || (nsym != n->sym)) {
                                              if (!chord.notes.isEmpty()) {
                                                    SymId sym = chord.notes.front()->sym;
                                                    if (sym == SymId::noteheadBlack)
                                                          chord.duration.setType(TDuration::DurationType::V_QUARTER);
                                                    else if (sym == SymId::noteheadHalf)
                                                          chord.duration.setType(TDuration::DurationType::V_HALF);
                                                    chords.append(chord);
                                                    chord.notes.clear();
                                                    }
                                              }
                                        nx = x;
                                        nsym = n->sym;
                                        chord.notes.append(n);
                                        }
                                  }
                            if (!chord.notes.isEmpty()) {
                                  SymId sym = chord.notes.front()->sym;
                                  if (sym == SymId::noteheadBlack)
                                        chord.duration.setType(TDuration::DurationType::V_QUARTER);
                                  else if (sym == SymId::noteheadHalf)
                                        chord.duration.setType(TDuration::DurationType::V_HALF);
                                  chords.append(chord);
                                  }
                            m.chords().append(chords);
                            }
                      system->measures().append(m);
                      }
                }

          }

    //---------------------------------------------------------
    //   searchClef
    //---------------------------------------------------------

    OmrClef OmrPage::searchClef(OmrSystem* system, OmrStaff* staff)
          {
          std::vector<Pattern*> pl = {
                Omr::trebleclefPattern,
                Omr::bassclefPattern
                };
          const OmrMeasure& m = system->measures().front();
          printf("search clef %d   %d-%d\n", staff->y(), m.x1(), m.x2());

          int x1 = m.x1() + 2;
          int x2 = x1 + (m.x2() - x1) / 2;
          OmrPattern p = searchPattern(pl, staff->y(), x1, x2);

          OmrClef clef(p);
          if (p.sym == SymId::gClef)
                clef.type = ClefType::G;
          else if (p.sym == SymId::fClef)
                clef.type = ClefType::F;
          else
                clef.type = ClefType::G;
          return clef;
          }

    //---------------------------------------------------------
    //   searchPattern
    //---------------------------------------------------------

    OmrPattern OmrPage::searchPattern(const std::vector<Pattern*>& pl, int y, int x1, int x2)
          {
          OmrPattern p;
          p.sym = SymId::noSym;
          p.prob = 0.0;
          for (Pattern* pattern : pl) {
                double val = 0.0;
                int xx = 0;
                int hw = pattern->w();

                for (int x = x1; x < (x2 - hw); ++x) {
                      double val1 = pattern->match(&image(), x - pattern->base().x(), y - pattern->base().y());
                      if (val1 > val) {
                            val = val1;
                            xx = x;
                            }
                      }

                if (val > p.prob) {
                      p.setRect(xx, y, pattern->w(), pattern->h());
                      p.sym = pattern->id();
                      p.prob = val;
                      }
                printf("Pattern found %d %f %d\n", int(pattern->id()), val, xx);
                }
          return p;
          }

    //---------------------------------------------------------
    //   searchTimeSig
    //---------------------------------------------------------

    OmrTimesig* OmrPage::searchTimeSig(OmrSystem* system)
          {
          int z = -1;
          int n = -1;
          double zval = 0;
          double nval = 0;
          QRect rz, rn;

          int y         = system->staves().front().y();
          OmrMeasure* m = &system->measures().front();
          int x1        = m->x1();

          for (int i = 0; i < 10; ++i) {
                Pattern* pattern = Omr::timesigPattern[i];
                double val = 0.0;
                int hh = pattern->h();
                int hw = pattern->w();
                int x2 = m->x2() - hw;
                QRect r;

                for (int x = x1; x < x2; ++x) {
                      double val1 = pattern->match(&image(), x, y);
                      if (val1 > val) {
                            val = val1;
                            r = QRect(x, y, hw, hh);
                            }
                      }

                if (val > timesigTH && val > zval) {
                      z = i;
                      zval = val;
                      rz = r;
                      }
                //            printf("   found %d %f\n", i, val);
                }

          if (z < 0)
                return 0;
          y = system->staves().front().y() + lrint(_spatium * 2);

          x1 = rz.x();
          int x2 = x1 + 1;
          OmrTimesig* ts = 0;
          for (int i = 0; i < 10; ++i) {
                Pattern* pattern = Omr::timesigPattern[i];
                double val = 0.0;
                int hh = pattern->h();
                int hw = pattern->w();
                QRect r;

                for (int x = x1; x < x2; ++x) {
                      double val1 = pattern->match(&image(), x, y);
                      if (val1 > val) {
                            val = val1;
                            r = QRect(x, y, hw, hh);
                            }
                      }

                if (val > timesigTH && val > nval) {
                      n = i;
                      nval = val;
                      rn = r;
                      }
                //            printf("   found %d %f\n", i, val);
                }
          if (n > 0) {
                ts = new OmrTimesig(rz | rn);
                ts->timesig = Fraction(z, n);
                printf("timesig  %d/%d\n", z, n);
                }
          return ts;
          }

    //---------------------------------------------------------
    //   searchKeySig
    //---------------------------------------------------------

    void OmrPage::searchKeySig(OmrSystem* system, OmrStaff* staff)
          {
          Pattern* pl[2];
          pl[0] = Omr::sharpPattern;
          pl[1] = Omr::flatPattern;

          double zval = 0;

          int y = system->staves().front().y();
          OmrMeasure* m = &system->measures().front();
          int x1 = m->x1();

          for (int i = 0; i < 2; ++i) {
                Pattern* pattern = pl[i];
                double val = 0.0;
                int hh = pattern->h();
                int hw = pattern->w();
                int x2 = m->x2() - hw;
                QRect r;

                for (int x = x1; x < x2; ++x) {
                      double val1 = pattern->match(&image(), x, y - hh / 2);
                      if (val1 > val) {
                            val = val1;
                            r = QRect(x, y, hw, hh);
                            }
                      }

                if (val > keysigTH && val > zval) {
                      zval = val;
                      OmrKeySig key(r);
                      key.type = i == 0 ? 1 : -1;
                      staff->setKeySig(key);
                      }
                printf(" key found %d %f\n", i, val);
                }
          }

    //---------------------------------------------------------
    //   maxP
    //---------------------------------------------------------

    int maxP(int* projection, int x1, int x2)
          {
          int xx = x1;
          int max = 0;
          for (int x = x1; x < x2; ++x) {
                if (projection[x] > max) {
                      max = projection[x];
                      xx = x;
                      }
                }
          return xx;
          }

    //---------------------------------------------------------
    //   BSTATE
    //---------------------------------------------------------
    struct BSTATE {
          int x;
          int status;//0 represents white space, 1 represents bar
          };

    //---------------------------------------------------------
    //   searchSysBarLines
    //---------------------------------------------------------

    void OmrSystem::searchSysBarLines()
          {
          OmrStaff& r1 = _staves[0];
          OmrStaff& r2 = _staves[_staves.size() - 1];//[1];

          int x1 = r1.x();
          int x2 = x1 + r1.width();
          int y1 = r1.y();
          int y2 = r2.y() + r2.height();
          int h = y2 - y1 + 1;
          int th = /*r1.height() + r2.height() - 5;*/ h / 2;     // threshold, data score for null model

          int vpw = x2 - x1;
          float vp[vpw];
          memset(vp, 0, sizeof(float) * vpw);

          //using note constraints
          searchNotes();

          int note_constraints[x2 - x1];
          for (int i = 0; i < _staves.size(); i++) {
                OmrStaff& r = _staves[i];
                for (OmrNote* n : r.notes()) {
                      for (int x = n->x(); x <= n->x() + n->width(); ++x)
                            note_constraints[x - x1] = 1;
                      }
                }

          //
          // compute vertical projections
          //
          for (int x = x1; x < x2; ++x) {
                int dots = 0;
                for (int y = y1; y < y2; ++y) {
                      if (_page->dot(x, y))
                            ++dots;
                      }
                if (!note_constraints[x - x1])
                      vp[x - x1] = dots - th;
                else
                      vp[x - x1] = -HUGE_VAL;
                }

          float scores[x2 - x1 + 1][2];
          BSTATE pred[x2 - x1 + 1][2];
          BSTATE bs;

          //initialization
          bs.x = -1; bs.status = -1;
          for (int x = x1; x <= x2; ++x) {
                int i = x - x1;
                for (int status = 0; status <= 1; ++status) {
                      scores[i][status] = -HUGE_VAL;
                      pred[i][status] = bs;
                      }
                }
          scores[0][0] = 0;

          //forward pass
          for (int x = x1; x < x2; ++x) {
                int i = x - x1;
                for (int cur = 0; cur <= 1; ++cur) {
                      //current state
                      if (scores[i][cur] < -1000) continue;

                      if (cur) {
                            scores[i][cur] += vp[i];
                            int next = 0;
                            int step = 3 * _page->spatium();//constraints between adjacent barlines
                            if (step + i <= x2 - x1) {
                                  if (scores[step + i][next] < scores[i][cur]) {
                                        scores[step + i][next] = scores[i][cur];
                                        bs.x = x; bs.status = cur;
                                        pred[step + i][next] = bs;
                                        }
                                  }
                            }
                      else {
                            for (int next = 0; next <= 1; ++next) {
                                  int step = 1;
                                  if (step + i <= x2 - x1) {
                                        if(scores[step + i][next] < scores[i][cur]) {
                                              scores[step + i][next] = scores[i][cur];
                                              bs.x = x; bs.status = cur;
                                              pred[step + i][next] = bs;
                                              }
                                        }
                                  }
                            }
                      }
                }

          //trace back
          int state = 0;
          for (int x = x2; x > x1; ) {
                int i = x - x1;
                bs = pred[i][state];
                state = bs.status;
                x = bs.x;

                if(state)
                      barLines.append(QLine(x, y1, x, y2));
                }
          }

    //-------------------------------------------------------------------
    //   searchBarLinesvar: dynamic programming for system identification
    //-------------------------------------------------------------------
    float OmrSystem::searchBarLinesvar(int n_staff, int **note_labels)
          {
          OmrStaff& r1 = _staves[0];
          OmrStaff& r2 = _staves[n_staff - 1];

          int x1 = r1.left();
          int x2 = r1.right();
          int y1 = r1.top();
          int y2 = r2.bottom();

          int th = (y2 - y1) / 2;//threshold

          //
          //compute note constraints
          //
          int *note_constraints = new int[x2 - x1 + 1];
          memset(note_constraints, 0, sizeof(int) * (x2 - x1 + 1));
          for (int i = 0; i < n_staff; ++i) {
                for (int j = x1; j <= x2; ++j) {
                      if (note_labels[i][j] == 1)
                            note_constraints[j - x1] = 1;
                      }
                }

          //
          // compute vertical projections
          //

          int vpw = x2 - x1 + 1;
          float vp[vpw];

          for(int i = 0; i < vpw; i++) vp[i] = -HUGE_VAL;

          int bar_max_width = 3;
          for (int x = x1 + bar_max_width; x < x2 - bar_max_width; ++x) {
                int dots = 0;
                for (int y = y1; y < y2; ++y) {
                      for (int w = -bar_max_width; w <= bar_max_width; w++) {
                            if (_page->dot(x + w, y)) {
                                  ++dots;
                                  break;
                                  }
                            }
                      }

                if (!note_constraints[x - x1])
                      vp[x - x1] = dots - th;
                else
                      vp[x - x1] = -HUGE_VAL;
                }

          float **scores = new float *[x2 - x1 + 1];
          BSTATE **pred = new BSTATE *[x2 - x1 + 1];
          for (int i = 0; i< x2 - x1 + 1; i++) {
                scores[i] = new float[2];
                pred[i] = new BSTATE[2];
                }

          BSTATE bs;

          //initialization
          bs.x = -1; bs.status = -1;
          for (int x = x1; x <= x2; ++x) {
                int i = x - x1;
                for (int status = 0; status <= 1; ++status) {
                      scores[i][status] = -HUGE_VAL;
                      pred[i][status] = bs;
                      }
                }
          scores[0][0] = 0;

          //forward pass

          for (int x = x1; x <= x2; ++x) {
                int i = x - x1;

                for(int cur = 0; cur <= 1; ++cur) {
                      //current state
                      if (scores[i][cur] < -1000)
                            continue;

                      if (cur) {
                            scores[i][cur] += vp[i];
                            int next = 0;
                            int step = 8 * _page->spatium();//minimum distance between adjacent barlines
                            if (step + i <= x2 - x1) {
                                  if (scores[step + i][next] < scores[i][cur]) {
                                        scores[step + i][next] = scores[i][cur];
                                        bs.x = x; bs.status = cur;
                                        pred[step + i][next] = bs;
                                        }
                                  }
                            else {
                                  if (scores[x2 - x1][next] < scores[i][cur]) {
                                        scores[x2 - x1][next] = scores[i][cur];
                                        bs.x = x; bs.status = cur;
                                        pred[x2 - x1][next] = bs;
                                        }
                                  }

                            }
                      else {
                            for (int next = 0; next <= 1; ++next) {
                                  int step = 1;
                                  if (step + i <= x2 - x1) {
                                        if (scores[step + i][next] < scores[i][cur]) {
                                              scores[step + i][next] = scores[i][cur];
                                              bs.x = x; bs.status = cur;
                                              pred[step + i][next] = bs;
                                              }
                                        }
                                  }
                            }
                      }
                }

          //trace back
          int state = 0;

          for (int x = x2; x > x1; ) {
                int i = x - x1;
                bs = pred[i][state];
                state = bs.status;
                x = bs.x;

                if(state) {
                      barLines.append(QLine(x, y1, x, y2));
                      }
                }

          float best_score = scores[x2 - x1][0];

          delete []note_constraints;
          for (int i = 0; i< x2 - x1 + 1; i++) {
                delete []scores[i];
                delete []pred[i];
                }
          delete []scores;
          delete []pred;

          return(best_score);
          }

    //---------------------------------------------------------
    //   noteCompare
    //---------------------------------------------------------

    static bool noteCompare(OmrNote* n1, OmrNote* n2)
          {
          return n1->x() < n2->x();
          }

    static bool intersectFuzz(const QRect& a, const QRect& b, int fuzz)
          {
          int xfuzz = fuzz;
          int yfuzz = fuzz;
          if (a.x() > b.x())
                xfuzz = -xfuzz;
          if (a.y() > b.y())
                yfuzz = -yfuzz;

          return (a.intersects(b.translated(xfuzz, yfuzz)));
          }

    //---------------------------------------------------------
    //   searchNotes
    //---------------------------------------------------------

    void OmrSystem::searchNotes()
          {
          for (int i = 0; i < _staves.size(); ++i) {
                OmrStaff* r = &_staves[i];
                int x1 = r->x();
                int x2 = x1 + r->width();

                //
                // search notes on a range of vertical line position
                //
                for (int line = 0; line < 8; ++line)
                      searchNotes(&r->notes(), x1, x2, r->y(), line);

                //
                // detect collisions
                //
                int fuzz = int(_page->spatium()) / 2;
                foreach(OmrNote* n, r->notes()) {
                      foreach(OmrNote* m, r->notes()) {
                            if (m == n)
                                  continue;
                            if (intersectFuzz(*m, *n, fuzz)) {
                                  if (m->prob > n->prob)
                                        r->notes().removeOne(n);
                                  else
                                        r->notes().removeOne(m);
                                  }
                            }
                      }
                qSort(r->notes().begin(), r->notes().end(), noteCompare);
                }
          }

    //---------------------------------------------------------
    //   searchNotes
    //---------------------------------------------------------

    void OmrSystem::searchNotes(int *note_labels, int ran)
          {
          for (int i = 0; i < _staves.size(); ++i) {
                OmrStaff* r = &_staves[i];
                int x1 = r->x();
                int x2 = x1 + r->width();

                //
                // search notes on a range of vertical line position
                //
                for (int line = -5; line < 14; ++line)
                      searchNotes(&r->notes(), x1, x2, r->y(), line);

                //
                // save detected note horizontal positions into note_labels
                //
                foreach(OmrNote* n, r->notes()) {
                      QPoint p = n->center();
                      int h_cent = p.x();
                      for(int h = h_cent - ran; h <= h_cent + ran; h++){
                            if(h >= x1 && h < x2) note_labels[h] = 1;
                            }
                      }
                }
          }

    //---------------------------------------------------------
    //   addText
    //---------------------------------------------------------

    #ifdef OCR
    static void addText(Score* score, int subtype, const QString& s)
    {
    #if 0 //TODO-1
          MeasureBase* measure = score->first();
          if (measure == 0 || measure->type() != Element::VBOX) {
                measure = new VBox(score);
                measure->setNext(score->first());
                measure->setTick(0);
                score->add(measure);
                }
          Text* text = new Text(score);
          switch(subtype) {
                case TEXT_TITLE:    text->setTextStyle(TEXT_STYLE_TITLE);    break;
                case TEXT_SUBTITLE: text->setTextStyle(TEXT_STYLE_SUBTITLE); break;
                case TEXT_COMPOSER: text->setTextStyle(TEXT_STYLE_COMPOSER); break;
                case TEXT_POET:     text->setTextStyle(TEXT_STYLE_POET);     break;
                }
          text->setSubtype(subtype);
          text->setText(s);
          measure->add(text);
    #endif
          }
    #endif

    //---------------------------------------------------------
    //   readHeader
    //---------------------------------------------------------

    void OmrPage::readHeader(Score*)
          {
          if (_slices.isEmpty())
                return;
          double maxHeight = _spatium * 4 * 2;

          int slice = 0;
          double maxH = 0.0;
          //      int maxIdx;
          for (; slice < _slices.size(); ++slice) {
                double h = _slices[slice].height();

                if (h > maxHeight)
                      break;
                if (h > maxH) {
                      maxH = h;
                      //                  maxIdx = slice;
                      }
                }
    #ifdef OCR
          //
          // assume that highest slice contains header text
          //
          OcrImage img = OcrImage(_image.bits(), _slices[maxIdx], (_image.width() + 31) / 32);
          QString s = _omr->ocr()->readLine(img).trimmed();
          if (!s.isEmpty())
                score->addText("title", s);

          QString subTitle;
          for (int i = maxIdx + 1; i < slice; ++i) {
                OcrImage img = OcrImage(_image.bits(), _slices[i], (_image.width() + 31) / 32);
                QString s = _omr->ocr()->readLine(img).trimmed();
                if (!s.isEmpty()) {
                      if (!subTitle.isEmpty())
                            subTitle += "\n";
                      subTitle += s;
                      }
                }
          if (!subTitle.isEmpty())
                score->addText("subtitle", subTitle);
    #endif
    #if 0
          OcrImage img = OcrImage(_image.bits(), _slices[0], (_image.width() + 31) / 32);
          QString s = _omr->ocr()->readLine(img).trimmed();
          if (!s.isEmpty())
                addText(score, TEXT_TITLE, s);

          img = OcrImage(_image.bits(), _slices[1], (_image.width() + 31) / 32);
          s = _omr->ocr()->readLine(img).trimmed();
          if (!s.isEmpty())
                addText(score, TEXT_SUBTITLE, s);

          img = OcrImage(_image.bits(), _slices[2], (_image.width() + 31) / 32);
          s = _omr->ocr()->readLine(img).trimmed();
          if (!s.isEmpty())
                addText(score, TEXT_COMPOSER, s);
    #endif
          }

    //---------------------------------------------------------
    //   removeBorder
    //---------------------------------------------------------

    void OmrPage::removeBorder()
          {
    #if 0
          cropT = cropB = cropL = cropR = 0;
          for (int y = 0; y < height(); ++y) {
                for (int x = 0; x < width(); ++x) {
                      QRgb c = _image.pixel(x,y);
                      if(qGray(c) > 0){
                            cropT = y;
                            break;
                            }
                      }
                if (cropT)
                      break;
                }
          for (int y = height() - 1; y >= cropT; --y) {
                for (int x = 0; x < width(); ++x) {
                      QRgb c = _image.pixel(x,y);
                      if(qGray(c) > 0){
                            cropB = y;
                            break;
                            }
                      }
                if (cropB)
                      break;
                }
          int y1 = cropT;
          int y2 = cropB;
          for (int x = 0; x < width(); ++x) {
                for (int y = y1; y < y2; ++y) {
                      QRgb c = _image.pixel(x,y);
                      if(qGray(c) > 0){
                            cropL = x;
                            break;
                            }
                      }
                if (cropL)
                      break;
                }
          for (int x = width() - 1; x >= cropL; --x) {
                for (int y = y1; y < y2; ++y) {
                      QRgb c = _image.pixel(x,y);
                      if(qGray(c) > 0){
                            cropR = x;
                            break;
                            }
                      }
                if (cropR)
                      break;
                }

          _image = _image.copy(cropL, cropT, cropR - cropL + 1, cropB - cropT + 1);
    #endif
          }

    //---------------------------------------------------------
    //   crop
    //---------------------------------------------------------

    void OmrPage::crop()
          {
          int wl = wordsPerLine();
          cropT = cropB = cropL = cropR = 0;
          for (int y = 0; y < height(); ++y) {
                const uint* p = scanLine(y);
                for (int k = 0; k < wl; ++k) {
                      if (*p++) {
                            cropT = y;
                            break;
                            }
                      }
                if (cropT)
                      break;
                }
          for (int y = height() - 1; y >= cropT; --y) {
                const uint* p = scanLine(y);
                for (int k = 0; k < wl; ++k) {
                      if (*p++) {
                            cropB = height() - y - 1;
                            break;
                            }
                      }
                if (cropB)
                      break;
                }
          int y1 = cropT;
          int y2 = height() - cropT - cropB;
          for (int x = 0; x < wl; ++x) {
                for (int y = y1; y < y2; ++y) {
                      if (*(scanLine(y) + x)) {
                            cropL = x;
                            break;
                            }
                      }
                if (cropL)
                      break;
                }
          for (int x = wl - 1; x >= cropL; --x) {
                for (int y = y1; y < y2; ++y) {
                      if (*(scanLine(y) + x)) {
                            cropR = wl - x - 1;
                            break;
                            }
                      }
                if (cropR)
                      break;
                }
          //      printf("*** crop: T%d B%d L%d R:%d\n", cropT, cropB, cropL, cropR);
          }

    //---------------------------------------------------------
    //   slice
    //---------------------------------------------------------

    void OmrPage::slice()
          {
          _slices.clear();
          int y1 = cropT;
          int y2 = height() - cropB;
          int x1 = cropL;
          int x2 = wordsPerLine() - cropR;
          int xbits = (x2 - x1) * 32;

          //      _slices.append(QRect(cropL*32, y1, xbits, y2-y1));
    #if 1
          for (int y = y1; y < y2;) {
                //
                // skip contents
                //
                int ys = y;
                for (; y < y2; ++y) {
                      const uint* p = scanLine(y) + cropL;
                      bool bits = false;
                      for (int x = cropL; x < x2; ++x) {
                            if (*p) {
                                  bits = true;
                                  break;
                                  }
                            ++p;
                            }
                      if (!bits)
                            break;
                      }
                if (y - ys)
                      _slices.append(QRect(cropL * 32, ys, xbits, y - ys));
                //
                // skip space
                //
                for (; y < y2; ++y) {
                      const uint* p = scanLine(y) + cropL;
                      bool bits = false;
                      for (int x = cropL; x < x2; ++x) {
                            if (*p) {
                                  bits = true;
                                  break;
                                  }
                            ++p;
                            }
                      if (bits)
                            break;
                      }
                }
    #endif
          }

    //---------------------------------------------------------
    //    deSkew
    //---------------------------------------------------------

    void OmrPage::deSkew()
          {
          int wl = wordsPerLine();
          int h = height();
          uint* db = new uint[wl * h];
          memset(db, 0, wl * h * sizeof(uint));

          foreach(const QRect& r, _slices) {
                double rot = skew(r);
                if (qAbs(rot) < 0.1) {
                      memcpy(db + wl * r.y(), scanLine(r.y()), wl * r.height() * sizeof(uint));
                      continue;
                      }

                QTransform t;
                t.rotate(rot);
                QTransform tt = QImage::trueMatrix(t, width(), r.height());

                double m11 = tt.m11();
                double m12 = tt.m12();
                double m21 = tt.m21();
                double m22 = tt.m22();
                double dx = tt.m31();
                double dy = tt.m32();

                double m21y = r.y() * m21;
                double m22y = r.y() * m22;
                int y2 = r.y() + r.height();

                for (int y = r.y(); y < y2; ++y) {
                      const uint* s = scanLine(y);

                      m21y += m21;
                      m22y += m22;
                      for (int x = 0; x < wl; ++x) {
                            uint c = *s++;
                            if (c == 0)
                                  continue;
                            uint mask = 1;
                            for (int xx = 0; xx < 32; ++xx) {
                                  if (c & mask) {
                                        int xs = x * 32 + xx;
                                        int xd = lrint(m11 * xs + m21y + dx);
                                        int yd = lrint(m22y + m12 * xs + dy);

                                        int wxd = xd / 32;
                                        if ((xd >= 0) && (wxd < wl) && (yd >= 0) && (yd < h)) {
                                              uint* d = db + wl * yd + wxd;
                                              if( d < db + wl * h) //check that we are in the bounds.
                                                    *d |= (0x1 << (xd % 32));
                                              }
                                        }
                                  mask <<= 1;
                                  }
                            }
                      }
                }
          memcpy(_image.bits(), db, wl * h * sizeof(uint));
          delete[] db;
          }

    struct ScanLine {
          int run;
          int x1, x2;
          ScanLine() { run = 0; x1 = 100000; x2 = 0; }
          };

    struct H {
          int y;
          int bits;
          H(int a, int b) : y(a), bits(b) {}
          };

    //---------------------------------------------------------
    //   xproject
    //---------------------------------------------------------

    int OmrPage::xproject(const uint* p, int wl)
          {
          int run = 0;
          int w = wl - cropL - cropR;
          int x1 = cropL + w / 4;         // only look at part of page
          int x2 = x1 + w / 2;
          for (int x = cropL; x < x2; ++x) {
                uint v = *p++;
                run += Omr::bitsSetTable[v & 0xff]
                   + Omr::bitsSetTable[(v >> 8) & 0xff]
                   + Omr::bitsSetTable[(v >> 16) & 0xff]
                   + Omr::bitsSetTable[v >> 24];
                }
          return run;
          }

    //---------------------------------------------------------
    //   xproject2
    //---------------------------------------------------------

    double OmrPage::xproject2(int y1)
          {
          int wl = wordsPerLine();
          const uint* db = bits();
          double val = 0.0;

          int w = wl - cropL - cropR;
          int x1 = (cropL + w / 4) * 32;       // only look at part of page
          int x2 = x1 + (w / 2 * 32);

          int ddx = x2 - x1;
          for (int dy = -12; dy < 12; ++dy) {
                int onRun = 0;
                int offRun = 0;
                int on = 0;
                int off = 0;
                bool onFlag = false;
                int incy = (dy > 0) ? 1 : (dy < 0) ? -1 : 0;
                int ddy = dy < 0 ? -dy : dy;
                int y = y1;
                if (y < 1)
                      y = 0;
                int err = ddx / 2;
                for (int x = x1; x < x2;) {
                      const uint* d = db + wl * y + (x / 32);
                      if (d < db + wl)
                            break;
                      if (d >= db + (wl - 1) * height()) //check that we are in the bounds.
                            break;
                      bool bit = ((*d) & (0x1 << (x % 32)));
                      bit = bit || ((*(d + wl)) & (0x1 << (x % 32)));
                      bit = bit || ((*(d - wl)) & (0x1 << (x % 32)));
                      if (bit != onFlag) {
                            if (!onFlag) {
                                  //
                                  // end of offrun:
                                  //
                                  if (offRun > 20) {
                                        off += offRun * offRun;
                                        on += onRun * onRun;
                                        onRun = 0;
                                        offRun = 0;
                                        }
                                  else {
                                        onRun += offRun;
                                        offRun = 0;
                                        }
                                  }
                            onFlag = bit;
                            }
                      (bit ? onRun : offRun)++;
                      if (offRun > 100) {
                            offRun = 0;
                            off = 1;
                            on = 0;
                            onRun = 0;
                            break;
                            }
                      err -= ddy;
                      if (err < 0) {
                            err += ddx;
                            y += incy;
                            if (y < 1)
                                  y = 1;
                            else if (y >= height())
                                  y = height() - 1;
                            }
                      ++x;
                      }
                if (offRun > 20)
                      off += offRun * offRun;
                else
                      onRun += offRun;
                on += onRun * onRun;
                if (off == 0)
                      off = 1;
                double nval = double(on) / double(off);
                if (nval > val)
                      val = nval;
                }
          return val;
          }

    static bool sortLvStaves(const Lv& a, const Lv& b)
          {
          return a.line < b.line;
          }

    //---------------------------------------------------------
    //   getStaffLines
    //---------------------------------------------------------

    void OmrPage::getStaffLines()
          {
          staves.clear();
          int h = height();
          int wl = wordsPerLine();
          // printf("getStaffLines %d %d  crop %d %d\n", h, wl, cropT, cropB);
          if (h < 1)
                return;

          int y1 = cropT;
          int y2 = h - cropB;
          if (y2 >= h)
                --y2;

          double projection[h];
          for (int y = 0; y < y1; ++y)
                projection[y] = 0;
          for (int y = y2; y < h; ++y)
                projection[y] = 0;
          for (int y = y1; y < y2; ++y)
                projection[y] = xproject2(y);
          int autoTableSize = (wl * 32) / 20;       // 1/20 page width
          if (autoTableSize > y2 - y1)
                autoTableSize = y2 - y1;
          double autoTable[autoTableSize];
          memset(autoTable, 0, sizeof(autoTable));
          for (int i = 0; i < autoTableSize; ++i) {
                autoTable[i] = covariance(projection + y1, projection + i + y1, y2 - y1 - i);
                }

          //
          // search for first maximum in covariance starting at 10 to skip
          // line width. Staff line distance (spatium) must be at least 10 dots
          //
          double maxCorrelation = 0;
          _spatium = 0;
          for (int i = 5; i < autoTableSize; ++i) {
                if (autoTable[i] > maxCorrelation) {
                      maxCorrelation = autoTable[i];
                      _spatium = i;
                      }
                }
          if (_spatium == 0) {
                printf("*** no staff lines found\n");
                return;
                }

          //---------------------------------------------------
          //    look for staves
          //---------------------------------------------------

          QList<Lv> lv;
          int ly = 0;
          int lval = -1000.0;
          for (int y = y1; y < (y2 - _spatium * 4); ++y) {
                double val = 0.0;
                for (int i = 0; i < 5; ++i)
                      val += projection[y + i * int(_spatium)];
                if (val < lval) {
                      lv.append(Lv(ly, lval));
                      }
                lval = val;
                ly = y;
                }
          qSort(lv);

          QList<Lv> staveTop;
          int staffHeight = _spatium * 6;
          foreach (Lv a, lv) {
                if (a.val < 500)   // MAGIC to avoid false positives
                      continue;
                int line = a.line;
                bool ok = true;
                foreach (Lv b, staveTop) {
                      if ((line > (b.line - staffHeight)) && (line < (b.line + staffHeight))) {
                            ok = false;
                            break;
                            }
                      }
                if (ok)
                      staveTop.append(a);
                }
          qSort(staveTop.begin(), staveTop.end(), sortLvStaves);
          foreach (Lv a, staveTop) {
                staves.append(OmrStaff(cropL * 32, a.line, (wordsPerLine() - cropL - cropR) * 32, _spatium * 4));
                }
          }

    //---------------------------------------------------------
    //   getRatio
    //---------------------------------------------------------

    void OmrPage::getRatio()
          {
          double num_black;
          double num_white;
          _ratio = 0.0;

          num_black = 1;
          num_white = 1;
          for (int x = 0; x < width(); ++x) {
                for (int y = 0; y < height(); ++y) {
                      if(isBlack(x,y)) num_black++;
                      else num_white++;
                      }
                }
          _ratio = num_black / (num_black + num_white);
          }

    //---------------------------------------------------------
    //   searchNotes
    //---------------------------------------------------------

    void OmrSystem::searchNotes(QList<OmrNote*>* noteList, int x1, int x2, int y, int line)
          {
          //a simple and cheap note detector (heuristic approach)
          double _spatium = _page->spatium();
          y += line * _spatium * .5;

          QList<Peak> notePeaks;
          Pattern* pattern = Omr::quartheadPattern;
          int hh = pattern->h();
          int hw = pattern->w();
          double val;
          int step_size = 2;
          int note_thresh = 50;

          for (int x = x1; x < (x2 - hw); x += step_size) {
                val = pattern->match(&_page->image(), x, y - hh / 2, _page->ratio());
                if (val > note_thresh) {
                      notePeaks.append(Peak(x, val, 0));
                      }
                }

          int n = notePeaks.size();
          for (int i = 0; i < n; ++i) {
                OmrNote* note = new OmrNote;
                int hh = pattern->h();
                int hw = pattern->w();
                note->setRect(notePeaks[i].x, y - hh / 2, hw, hh);
                note->line = line;
                note->sym = pattern->id();
                note->prob = notePeaks[i].val;
                noteList->append(note);
                }
          return;
          }

    //---------------------------------------------------------
    //   staffDistance
    //---------------------------------------------------------

    double OmrPage::staffDistance() const
          {
          if (staves.size() < 2)
                return 5;
          return ((staves[1].y() - staves[0].y()) / _spatium) - 4.0;
          }

    //---------------------------------------------------------
    //   systemDistance
    //---------------------------------------------------------

    double OmrPage::systemDistance() const
          {
          if (staves.size() < 3)
                return 6.0;
          return ((staves[2].y() - staves[1].y()) / _spatium) - 4.0;
          }

    //---------------------------------------------------------
    //   write
    //---------------------------------------------------------

    void OmrPage::write(XmlWriter& xml) const
          {
          xml.stag("OmrPage");
          xml.tag("cropL", cropL);
          xml.tag("cropR", cropR);
          xml.tag("cropT", cropT);
          xml.tag("cropB", cropB);
          foreach(const QRect& r, staves)
                xml.tag("staff", QRectF(r));
          xml.etag();
          }

    //---------------------------------------------------------
    //   read
    //---------------------------------------------------------

    void OmrPage::read(XmlReader& e)
          {
          while (e.readNextStartElement()) {
                const QStringRef& tag(e.name());

                if (tag == "cropL")
                      cropL = e.readInt();
                else if (tag == "cropR")
                      cropR = e.readInt();
                else if (tag == "cropT")
                      cropT = e.readInt();
                else if (tag == "cropB")
                      cropB = e.readInt();
                else if (tag == "staff") {
                      OmrStaff r(e.readRect().toRect());
                      staves.append(r);
                      }
                else
                      e.unknown();
                }
          }
    }

↑ V519 The '_ratio' variable is assigned values twice successively. Perhaps this is a mistake. Check lines: 1551, 1561.

↑ V730 Not all members of a class are initialized inside the constructor. Consider inspecting: _spatium, _ratio.

↑ V730 Not all members of a class are initialized inside the constructor. Consider inspecting: sym.