aligner.cpp

本程序是计算机视觉稠密匹配的程序

/*
 *
 *  All Contents Copyright 2000 by Jared Samet. All Rights Reserved.
 *
 */

#include <stdafx.h>
#include "sequence.h"
#include "aligner.h"
#include "SteppedArray.h"

CAligner::CAligner(int iGapStart, int iGapContinue, MATCH Match, 
                   int iMinDisp, int iMaxDisp)
    : m_iGapStart(iGapStart), m_iGapContinue(iGapContinue),
      m_Match(Match), m_iMinDisp(iMinDisp), m_iMaxDisp(iMaxDisp),
      m_fAligned(false), m_fHaveSeqs(false), m_pSeqLeft(NULL), m_pSeqRight(NULL)
{
}

int CAligner::AlignmentCount()
{
    return m_fAligned ? 1 : 0;
}

bool CAligner::isAligned()
{
    return m_fAligned;
}

bool CAligner::SetSequences(CSequence *pSeq1, CSequence *pSeq2)
{
    if (m_fHaveSeqs && 
        (m_iLeftLen != (int) pSeq1->Length() || m_iRightLen != (int) pSeq2->Length()))
    {
        return false;
    }

    m_pSeqLeft = pSeq1;
    m_pSeqRight = pSeq2;

    if (!m_fHaveSeqs)
    {
        m_iLeftLen = (int) pSeq1->Length();
        m_iRightLen = (int) pSeq2->Length();
        m_pTable = new CSteppedArray<SDPBox>(m_iLeftLen+1, m_iRightLen+1, m_iMinDisp, m_iMaxDisp);     
        if (m_pTable->isValid(0,0))
        {
            m_pTable->ElementAt(0,0).A.iScore = 0;
        }
        int i;
        for (i=1; i <= m_iLeftLen; i++)
        {
            if (m_pTable->isValid(0,i))
            {
                m_pTable->ElementAt(0,i).B.iScore = -(m_iGapStart + i*m_iGapContinue);
            }
        }
        for (i=1; i <= m_iRightLen; i++)
        {
            if (m_pTable->isValid(i,0))
            {
                m_pTable->ElementAt(i,0).C.iScore = -(m_iGapStart + i*m_iGapContinue);
            }
        }
    }

    m_fHaveSeqs = true;
    return true;
}

bool CAligner::isValid(C3D_ARRAY arPrev, int iRow, int iCol)
{
    if (C3D_B == arPrev && 0 == iRow)
        return true;
    if (C3D_C == arPrev && 0 == iCol)
        return true;
    if  ((iRow > m_iLeftLen) || (iCol > m_iRightLen) ||
        (0 > iRow) || (0 > iCol) ||
        (m_iMinDisp > (iRow - iCol)) || ((iRow - iCol) > m_iMaxDisp))
        return false;
    switch(arPrev)
    {
    case C3D_A:
        return (0 == iRow && 0 == iCol) || (0 != iRow && 0 != iCol);
        break;
    case C3D_B:
        if (0 == iCol)
            return false;
        if (m_iMaxDisp == (iRow - iCol))
            return false;
        return true;
        break;
    case C3D_C:
        if (0 == iRow)
            return false;
        if (m_iMinDisp == (iRow - iCol))
            return false;
        return true;
        break;
    default:
        ASSERT(false);
        return false;
        break;
    }
}

void CAligner::Align()
{
    int i, d, j, iMatch;
    ASSERT(m_pSeqLeft && m_pSeqRight);
    
    for (i=1; i <= m_iLeftLen; i++)
    {
        for (d = m_iMaxDisp; d >= m_iMinDisp; d--)
        {
            j = i - d;
            if (1 <= j && j <= m_iRightLen)
            {
                // compute all three boxes of the table, if they are valid.
                if (isValid(C3D_A, i, j))
                {
                    iMatch = m_Match(m_pSeqLeft->LetterAt(i-1), m_pSeqRight->LetterAt(j-1));
                    if (isValid(C3D_A, i-1, j-1))
                    {
                        m_pTable->ElementAt(i,j).A.iScore = iMatch + m_pTable->ElementAt(i-1,j-1).A.iScore;
                        m_pTable->ElementAt(i,j).A.arPrev = C3D_A;
                        m_pTable->ElementAt(i,j).A.dirPrev = C3D_NORTHWEST;
                    }
                    else
                    {
                        m_pTable->ElementAt(i,j).A.iScore = INT_MIN;
                        m_pTable->ElementAt(i,j).A.arPrev = C3D_INVALID_AR;
                        m_pTable->ElementAt(i,j).A.dirPrev = C3D_INVALID_DIR;
                    }
                    if (isValid(C3D_B, i-1, j-1) &&
                        (m_pTable->ElementAt(i,j).A.iScore < iMatch + m_pTable->ElementAt(i-1,j-1).B.iScore))
                    {
                        m_pTable->ElementAt(i,j).A.iScore = iMatch + m_pTable->ElementAt(i-1,j-1).B.iScore;
                        m_pTable->ElementAt(i,j).A.arPrev = C3D_B;
                        m_pTable->ElementAt(i,j).A.dirPrev = C3D_NORTHWEST;
                    }
                    if (isValid(C3D_C, i-1, j-1) &&
                        (m_pTable->ElementAt(i,j).A.iScore < iMatch + m_pTable->ElementAt(i-1,j-1).C.iScore))
                    {
                        m_pTable->ElementAt(i,j).A.iScore = iMatch + m_pTable->ElementAt(i-1,j-1).C.iScore;
                        m_pTable->ElementAt(i,j).A.arPrev = C3D_C;
                        m_pTable->ElementAt(i,j).A.dirPrev = C3D_NORTHWEST;
                    }
                }
                if (isValid(C3D_B, i, j))
                {
                    if (isValid(C3D_A, i, j-1))
                    {
                        m_pTable->ElementAt(i,j).B.iScore = -(m_iGapStart + m_iGapContinue) + m_pTable->ElementAt(i,j-1).A.iScore;
                        m_pTable->ElementAt(i,j).B.arPrev = C3D_A;
                        m_pTable->ElementAt(i,j).B.dirPrev = C3D_WEST;
                    }
                    else
                    {
                        m_pTable->ElementAt(i,j).B.iScore = INT_MIN;
                        m_pTable->ElementAt(i,j).B.arPrev = C3D_INVALID_AR;
                        m_pTable->ElementAt(i,j).B.dirPrev = C3D_INVALID_DIR;
                    }
                    if (isValid(C3D_B, i, j-1) &&
                        (m_pTable->ElementAt(i,j).B.iScore < -(m_iGapContinue) + m_pTable->ElementAt(i,j-1).B.iScore))
                    {
                        m_pTable->ElementAt(i,j).B.iScore = -(m_iGapContinue) + m_pTable->ElementAt(i,j-1).B.iScore;
                        m_pTable->ElementAt(i,j).B.arPrev = C3D_B;
                        m_pTable->ElementAt(i,j).B.dirPrev = C3D_WEST;
                    }
                    if (isValid(C3D_C, i, j-1) &&
                        (m_pTable->ElementAt(i,j).B.iScore < -(m_iGapStart + m_iGapContinue) + m_pTable->ElementAt(i,j-1).C.iScore))
                    {
                        m_pTable->ElementAt(i,j).B.iScore = -(m_iGapStart + m_iGapContinue) + m_pTable->ElementAt(i,j-1).C.iScore;
                        m_pTable->ElementAt(i,j).B.arPrev = C3D_C;
                        m_pTable->ElementAt(i,j).B.dirPrev = C3D_WEST;
                    }
                }
                if (isValid(C3D_C, i, j))
                {
                    if (isValid(C3D_A, i-1, j))
                    {
                        m_pTable->ElementAt(i,j).C.iScore = -(m_iGapStart + m_iGapContinue) + m_pTable->ElementAt(i-1,j).A.iScore;
                        m_pTable->ElementAt(i,j).C.arPrev = C3D_A;
                        m_pTable->ElementAt(i,j).C.dirPrev = C3D_NORTH;
                    }
                    else
                    {
                        m_pTable->ElementAt(i,j).C.iScore = INT_MIN;
                        m_pTable->ElementAt(i,j).C.arPrev = C3D_INVALID_AR;
                        m_pTable->ElementAt(i,j).C.dirPrev = C3D_INVALID_DIR;
                    }
                    if (isValid(C3D_B, i-1, j) &&
                        (m_pTable->ElementAt(i,j).C.iScore < -(m_iGapStart + m_iGapContinue) + m_pTable->ElementAt(i-1,j).B.iScore))
                    {
                        m_pTable->ElementAt(i,j).C.iScore = -(m_iGapStart + m_iGapContinue) + m_pTable->ElementAt(i-1,j).B.iScore;
                        m_pTable->ElementAt(i,j).C.arPrev = C3D_B;
                        m_pTable->ElementAt(i,j).C.dirPrev = C3D_NORTH;
                    }
                    if (isValid(C3D_C, i-1, j) &&
                        (m_pTable->ElementAt(i,j).C.iScore < -(m_iGapContinue) + m_pTable->ElementAt(i-1,j).C.iScore))
                    {
                        m_pTable->ElementAt(i,j).C.iScore = -(m_iGapContinue) + m_pTable->ElementAt(i-1,j).C.iScore;
                        m_pTable->ElementAt(i,j).C.arPrev = C3D_C;
                        m_pTable->ElementAt(i,j).C.dirPrev = C3D_NORTH;
                    }
                }
            }
        }
    }
    m_fAligned = true;
}

void CAligner::GetAlignment(int nAlignment, CSequence **ppSeq1, CSequence **ppSeq2)
{
    ASSERT(0 == nAlignment && m_fAligned);

    int iRow = m_iLeftLen, jCol = m_iRightLen;
    int cLeftSpc=0, cRightSpc=0;
    int iMaxScore;
    C3D_ARRAY arCur = C3D_A, arStart;
    
    while ((iRow - jCol) < m_iMinDisp)
    {
        cLeftSpc++;
        jCol--;
    }
    while ((iRow - jCol) > m_iMaxDisp)
    {
        cRightSpc++;
        iRow--;
    }

    ASSERT(isValid(arCur, iRow, jCol));
    iMaxScore = m_pTable->ElementAt(iRow, jCol).A.iScore;
    if (isValid(C3D_B, iRow, jCol) && m_pTable->ElementAt(iRow, jCol).B.iScore > iMaxScore)
    {
        iMaxScore = m_pTable->ElementAt(iRow, jCol).B.iScore;
        arCur = C3D_B;
    }
    if (isValid(C3D_C, iRow, jCol) && m_pTable->ElementAt(iRow, jCol).C.iScore > iMaxScore)
    {
        iMaxScore = m_pTable->ElementAt(iRow, jCol).C.iScore;
        arCur = C3D_C;
    }
    arStart = arCur;

    while (iRow > 0 || jCol > 0)
    {
        ASSERT(isValid(arCur, iRow, jCol));
        if (0 == iRow)
        {
            cLeftSpc++;
            jCol--;
            continue;
        }
        if (0 == jCol)
        {
            cRightSpc++;
            iRow--;
            continue;
        }
        switch(arCur)
        {
        case C3D_A:
            arCur = m_pTable->ElementAt(iRow, jCol).A.arPrev;
            iRow--;
            jCol--;
            break;
        case C3D_B:
            arCur = m_pTable->ElementAt(iRow, jCol).B.arPrev;
            cLeftSpc++;
            jCol--;
            break;
        case C3D_C:
            arCur = m_pTable->ElementAt(iRow, jCol).C.arPrev;
            cRightSpc++;
            iRow--;
            break;
        default:
            ASSERT(false);
        }
    }

    ASSERT(m_iLeftLen + cLeftSpc == m_iRightLen + cRightSpc);
    CSequence *pAlignLeft = new CSequence(NULL, m_iLeftLen + cLeftSpc);
    CSequence *pAlignRight = new CSequence(NULL, m_iRightLen + cRightSpc);
    int iCurPos = m_iLeftLen + cLeftSpc - 1;

    iRow = m_iLeftLen;
    jCol = m_iRightLen;    
    arCur = arStart;
    
    while ((iRow - jCol) < m_iMinDisp)
    {
        pAlignLeft->LetterAt(iCurPos) = 0xFFFFFFFF;
        pAlignRight->LetterAt(iCurPos) = m_pSeqRight->LetterAt(jCol - 1);
        jCol--;
        iCurPos--;
    }
    while ((iRow - jCol) > m_iMaxDisp)
    {
        pAlignLeft->LetterAt(iCurPos) = m_pSeqLeft->LetterAt(iRow - 1);
        pAlignRight->LetterAt(iCurPos) = 0xFFFFFFFF;
        iRow--;
        iCurPos--;
    }

    while (iRow > 0 || jCol > 0)
    {
        ASSERT(isValid(arCur, iRow, jCol));
        if (0 == iRow)
        {
            pAlignLeft->LetterAt(iCurPos) = 0xFFFFFFFF;
            pAlignRight->LetterAt(iCurPos) = m_pSeqRight->LetterAt(jCol - 1);
            iCurPos--;
            jCol--;
            continue;
        }
        if (0 == jCol)
        {
            pAlignLeft->LetterAt(iCurPos) = m_pSeqLeft->LetterAt(iRow - 1);
            pAlignRight->LetterAt(iCurPos) = 0xFFFFFFFF;
            iRow--;
            iCurPos--;
            continue;
        }
        switch(arCur)
        {
        case C3D_A:
            arCur = m_pTable->ElementAt(iRow, jCol).A.arPrev;
            pAlignLeft->LetterAt(iCurPos) = m_pSeqLeft->LetterAt(iRow - 1);
            pAlignRight->LetterAt(iCurPos) = m_pSeqRight->LetterAt(jCol - 1);
            iRow--;
            jCol--;
            iCurPos--;
            break;
        case C3D_B:
            arCur = m_pTable->ElementAt(iRow, jCol).B.arPrev;
            pAlignLeft->LetterAt(iCurPos) = 0xFFFFFFFF;
            pAlignRight->LetterAt(iCurPos) = m_pSeqRight->LetterAt(jCol - 1);
            iCurPos--;
            jCol--;
            break;
        case C3D_C:
            arCur = m_pTable->ElementAt(iRow, jCol).C.arPrev;
            pAlignLeft->LetterAt(iCurPos) = m_pSeqLeft->LetterAt(iRow - 1);
            pAlignRight->LetterAt(iCurPos) = 0xFFFFFFFF;
            iRow--;
            iCurPos--;
            break;
        default:
            ASSERT(false);
        }
    }
    ASSERT(-1 == iCurPos);
    *ppSeq1 = pAlignLeft;
    *ppSeq2 = pAlignRight;
}

CAligner::~CAligner()
{
    delete m_pTable;
}