trace_graph.cpp

来自「ncbi源码」· C++ 代码 · 共 926 行 · 第 1/3 页

    TModelUnit y3 = y + top_h + 1;
    TModelUnit y4 = y + total_h - 1;

    TSignedSeqRange gr_range = TSignedSeqRange(m_pData->GetSeqFrom(), m_pData->GetSeqTo());
    SChunkTranslator    trans;
    
    for( int i = 0; i < chunks->size(); i++  )   { // iterate by chunks
        CConstRef<CAlnMap::CAlnChunk> ch = (*chunks)[i];
        
        if(ch->GetType() & CAlnMap::fSeq) { // if chunk aligned
            TSignedSeqRange ch_range = ch->GetRange();  
            const TSignedSeqRange& ch_aln_range = ch->GetAlnRange();
            trans.Init(ch_range, ch_aln_range, m_pData->IsNegative());
            
            // intersect 
            ch_range.IntersectWith(gr_range);
            
            if(ch_range.NotEmpty())   {           
                
                double x1 = trans.GetAlnPosFromSeqPos(ch_range.GetFrom()) - offset_x;
                double x2 = trans.GetAlnPosFromSeqPos(ch_range.GetTo()) - offset_x;
                if(x1 > x2) {
                    swap(x1, x2);
                }
                x2 += 1.0;
                glRectd(x1, y1, x2, y2);
                glRectd(x1, y3, x2, y4);
            }
        }
    }
}


const static float kEps = 0.000001f; // precision for float point comparisions

// renders confidence graph
void    CTraceGraph::x_RenderConfGraph(CGlPane& pane, int y, int h, 
                              const TChunkVec& chunks)
{
    const TVPRect& rc_vp = pane.GetViewport();
    const TModelRect rc_vis = pane.GetVisibleRect();

    glDisable(GL_BLEND);   
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
    
    TModelUnit amp  = ((TModelUnit) h - 2 * kGraphOffsetY -1)
                        / m_pData->GetMaxConfidence();
    TModelUnit base_y = y + kGraphOffsetY;
    TModelUnit scale_x = pane.GetScaleX();

    TSignedSeqPos gr_from = m_pData->GetSeqFrom();
    TSignedSeqPos gr_to = m_pData->GetSeqTo();
    
    SChunkTranslator trans;

    bool b_average = (scale_x > 0.3);
    if(b_average)    {
        pane.Close();
        pane.OpenPixels();

        //calculate visible range in seq coords
        TVPUnit range_pix_start = rc_vp.Left();
        TVPUnit range_pix_end = rc_vp.Right();
        
        TModelUnit  left = rc_vis.Left();
        TModelUnit  scale_x = pane.GetScaleX();
        TModelUnit  top_y = rc_vp.Bottom() - (y + kGraphOffsetY);

        glBegin(GL_LINES);

        int i_prev_chunk = -1, i_chunk = 0, n_chunks = chunks->size();
        CConstRef<CAlnMap::CAlnChunk> ch;
        CAlnMap::TSegTypeFlags type;
        TModelUnit ch_pix_end;
        TModelUnit pos_pix_left, pos_pix_right;

        TSignedSeqPos from = -1, to = -2;
        TSignedSeqPos pos = gr_from, pos_inc = 1;
        double v = 0, v_min = 0, v_max = 0;
            
        for( int pix = range_pix_start; 
             pix <= range_pix_end  &&  i_chunk < n_chunks  &&  pos <= gr_to; )   { // iterate by pixels
            bool b_first = true;

            // if chunks ends before pixel 
            if(i_prev_chunk != i_chunk) {
                // advance chunk
                ch = (*chunks)[i_chunk];
                type = ch->GetType();

                if(type & CAlnMap::fSeq) {             
                    const TSignedSeqRange& ch_range = ch->GetRange();  
                    const TSignedSeqRange& ch_aln_range = ch->GetAlnRange();
                
                    trans.Init(ch_range, ch_aln_range, m_pData->IsNegative());
                
                    // calculate intersection of chunk in seq_coords with graph range
                    from = max(ch_range.GetFrom(), gr_from);
                    to = min(ch_range.GetTo(), gr_to);
                
                    // clip with visible align range (expanding integer clip)
                    TSignedSeqPos vis_from = trans.GetSeqPosFromAlnPos((TSignedSeqPos) floor(rc_vis.Left()));
                    TSignedSeqPos vis_to = trans.GetSeqPosFromAlnPos((TSignedSeqPos) ceil(rc_vis.Right()));
                    if(vis_to < vis_from)   {
                        _ASSERT(m_pData->IsNegative());
                        swap(vis_from, vis_to);
                    }
                    from = max(from, vis_from);
                    to = min(to, vis_to);
                    if(m_pData->IsNegative()) {
                        ch_pix_end = range_pix_start + (trans.GetAlnPosFromSeqPos(from) + 1 - left) / scale_x;
                        pos = to;
                        pos_inc = -1;
                    } else {
                        ch_pix_end = range_pix_start + (trans.GetAlnPosFromSeqPos(to) + 1 - left) / scale_x;
                        pos = from;
                        pos_inc = 1;
                    }
                }
                i_prev_chunk = i_chunk;
            }
            
            if(type & CAlnMap::fSeq) {         
                TSignedSeqPos aln_pos = trans.GetAlnPosFromSeqPos(pos);
                pos_pix_left = range_pix_start + (aln_pos - left) / scale_x;
                pos_pix_right = pos_pix_left + 1 / scale_x;
                
                pos_pix_left = max((TModelUnit) pix, pos_pix_left);                

                _ASSERT(pos_pix_left >= pix);

                while(pos >= from  &&  pos <= to  &&  pos_pix_left < pix + 1)   {                
                    // calculate overlap with pixel and integrate
                    double v = amp * m_pData->GetConfidence(pos);
                    
                    v_min = b_first ? v : min(v_min, v);
                    v_max = b_first ? v : max(v_max, v);
                    b_first = false;
                    
                    if(pos_pix_right < pix +1)  {
                        pos +=pos_inc; // advance to next pos
                        aln_pos = trans.GetAlnPosFromSeqPos(pos);
                        pos_pix_left = range_pix_start + (aln_pos - left) / scale_x;
                        pos_pix_right = pos_pix_left + 1 / scale_x;
                    } else break;                    
                }
            } 
            if(ch_pix_end < pix + 1) {
                i_chunk++;
            } else  {           
                if(type & CAlnMap::fSeq) { 
                    glColor3d(0.0f, 0.5f, 0.0f);                
                    glVertex2d(pix, top_y);
                    glVertex2d(pix, top_y - v_min);

                    glColor3d(0.0f, 0.75f, 0.25f);
                    glVertex2d(pix, top_y - v_min);
                    glVertex2d(pix, top_y - v_max);
                
                }
                pix++;
                v = v_min = v_max = 0;
            }
        }    
        glEnd();
    } else { // render without averaging
        _ASSERT(pane.GetProjMode() == CGlPane::eOrtho);
        
        TModelUnit offset_x = pane.GetOffsetX();
        glColor3d(0.0f, 0.5f, 0.0f); //###

        for( int i = 0; i < chunks->size(); i++  )   { // iterate by chunks
            CConstRef<CAlnMap::CAlnChunk> ch = (*chunks)[i];
            CAlnMap::TSegTypeFlags type = ch->GetType();

            if(type & CAlnMap::fSeq) { // if chunk aligned
                const TSignedSeqRange& ch_range = ch->GetRange();  
                const TSignedSeqRange& ch_aln_range = ch->GetAlnRange();
                trans.Init(ch_range, ch_aln_range, m_pData->IsNegative());

                TSignedSeqPos from = max(gr_from, ch_range.GetFrom());
                TSignedSeqPos to = min(gr_to, ch_range.GetTo());
           
                for( TSignedSeqPos pos = from; pos <= to; pos ++ ) {
                    double v = m_pData->GetConfidence(pos);
                    v *= amp;
                    TSignedSeqPos aln_pos = trans.GetAlnPosFromSeqPos(pos);
                    double x = aln_pos - offset_x;
                    glRectd(x , base_y, x + 1.0, base_y + v);
                }
            }
        }
    }
}

const static int kIntBandSpace = 1;

void    CTraceGraph::x_RenderSignalGraph(CGlPane& pane, int y, int h, 
                                const TChunkVec& chunks)
{
    CTraceData::TSignalValue MaxSignal = 0;
    for( int ch = CTraceData::eA; ch <= CTraceData::eG; ch++ )   {
        MaxSignal = max(MaxSignal, m_pData->GetMax( (CTraceData::EChannel) ch));
    }
    TModelUnit amp  = ((TModelUnit) (h - 2 * kGraphOffsetY)) / MaxSignal; 
    int bottom_y = y + (h - 1) - kGraphOffsetY;
    
    for( int i_ch = 0; i_ch < 4; i_ch++ )   {   // for every channel (i_ch - channel index)
        int ch_index = (m_Props.m_bReverseColors  &&  m_pData->IsNegative()) ? (i_ch ^ 2) : i_ch;
        int ch = CTraceData::eA + ch_index;
    
        const CTraceData::TPositions& positions = m_pData->GetPositions();
        
        CTraceData::EChannel channel = (CTraceData::EChannel) ch;
        const CTraceData::TValues& values = m_pData->GetValues(channel);

        glColorC(m_vSignalColors[kGradColors * (i_ch + 1) -1]);
        glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);

        _ASSERT(m_Props.m_SignalStyle == CTraceGraphProperties::eCurve);

        for( int i = 0; i < chunks->size(); i++  )   {
            CConstRef<CAlnMap::CAlnChunk> ch = (*chunks)[i];
            CAlnMap::TSegTypeFlags type = ch->GetType();

            if(type & CAlnMap::fSeq) {
                x_RenderCurveChunk(pane, ch, positions, values, bottom_y, h, (int) amp);            
            }
        }
    }
}

// Renders chromatogram corresponding to given chunk of sequence.
// This function renders data as a curve for a single channel specified by "values".
void    CTraceGraph::x_RenderCurveChunk(CGlPane& pane, 
                                          CConstRef<CAlnMap::CAlnChunk> ch, 
                                          const CTraceData::TPositions& positions, 
                                          const CTraceData::TValues& values,
                                          int bottom_y, int h, int amp)
{
    const TModelRect rc_vis = pane.GetVisibleRect();

    const TSignedSeqRange& ch_range = ch->GetRange();  
    const TSignedSeqRange& ch_aln_range = ch->GetAlnRange();

    SChunkTranslator trans;
    bool b_neg = m_pData->IsNegative();
    trans.Init(ch_range, ch_aln_range, m_pData->IsNegative());
    
    // [from, to] - is sequence range for which graph is rendered
    double from = max(m_pData->GetSeqFrom(), ch_range.GetFrom());
    double to = min(m_pData->GetSeqTo(), ch_range.GetTo());
    
    double vis_from = trans.GetSeqPosFromAlnPos(rc_vis.Left());
    double vis_to = trans.GetSeqPosFromAlnPos(rc_vis.Right());
    if(vis_to < vis_from)   {
        _ASSERT(m_pData->IsNegative());
        swap(vis_from, vis_to);
    }

    from = max(from, vis_from);
    to = min(to, vis_to);

    if(from <= to)  {
        int sm_start = x_FindSampleToLeft(from);
        int sm_end = x_FindSampleToRight(to + 1);
        sm_start = max(0, sm_start);
        sm_end = min(sm_end, m_pData->GetSamplesCount() - 1);
   
        if(sm_start <= sm_end)  {
            glBegin(GL_QUAD_STRIP);
            TModelUnit offset_x = pane.GetOffsetX();
         
            // check if start interpolation is needed
            CTraceData::TFloatSeqPos sm_start_seqpos = positions[sm_start];
            if(sm_start_seqpos < from)   {
                if(sm_start + 1 < sm_end)   {
                    double v1 = values[sm_start];
                    double v2 = values[sm_start + 1];
                    double x1 = sm_start_seqpos;
                    double x2 = positions[sm_start + 1];
                    double v = v1 + ((from - x1) * (v2 - v1) / (x2 - x1));
                    v *= amp;
                    
                    double aln_from = trans.GetAlnPosFromSeqPos(from + (b_neg ? -1 : 0));
                    glVertex2d(aln_from - offset_x, bottom_y - v); 
                    //glVertex2d(aln_from - offset_x, bottom_y);                 
                }
                sm_start++;
            }
            // render regular samples
            for( int i_sm = sm_start; i_sm < sm_end; i_sm++  ) {
                TModelUnit seqpos = positions[i_sm];
                double v = values[i_sm];
                v *= amp;
        
                double aln_pos = trans.GetAlnPosFromSeqPos(seqpos + (b_neg ? -1 : 0));
                glVertex2d(aln_pos - offset_x, bottom_y - v); 
                //glVertex2d(aln_pos - offset_x, bottom_y);                 
            }
            // render end point
            if( sm_end - 1 > sm_start) { // interpolate end point
                double v1 = values[sm_end -1];
                double v2 = values[sm_end];
                double x1 = positions[sm_end - 1];
                double x2 = positions[sm_end];
                double v = v1 + ((to + 1 - x1) * (v2 - v1) / (x2 - x1));
                v *= amp;