trace_graph.cpp

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

                double aln_to = trans.GetAlnPosFromSeqPos(to + (b_neg ? 0 : 1));
                glVertex2d(aln_to - offset_x, bottom_y - v); 
                //glVertex2d(aln_to - offset_x, bottom_y);  
            } else {
               TModelUnit seqpos = positions[sm_end];
                double v = values[sm_end];
                v *= amp;
        
                double aln_pos = trans.GetAlnPosFromSeqPos(seqpos);
                glVertex2d(aln_pos - offset_x, bottom_y - v); 
                //glVertex2d(aln_pos - offset_x, bottom_y);                 
            }
            glEnd();
        }
    }
}


/// Render signals for all channels as gradient-color bands with color intensity
/// proprotional to signal strength.
void    CTraceGraph::x_RenderIntensityGraphs(CGlPane& pane, int y, int h, 
                                             const TChunkVec& chunks)
{   
    //_TRACE("\nx_RenderIntensityGraphs");
    const CTraceData::TPositions& positions = m_pData->GetPositions();
    
    const TVPRect& rc_vp = pane.GetViewport();
    const TModelRect rc_vis = pane.GetVisibleRect();
  
    // calculate layout
    int av_h = h - 2 * kGraphOffsetY; // height available for intensity bands
    int band_h = av_h / 4;
    int off = (av_h - 4 * band_h) / 2; // rounding error compensation
    int top_y = rc_vp.Bottom() - (y + kGraphOffsetY + off);
    
    TModelUnit  left = rc_vis.Left();
    TModelUnit  scale_x = pane.GetScaleX();

    TVPUnit range_pix_start = rc_vp.Left();
    TVPUnit range_pix_end = rc_vp.Right();

    SChunkTranslator trans;

    glBegin(GL_LINES);
    
    int i_chunk = 0, n_chunks = chunks->size();            
    CConstRef<CAlnMap::CAlnChunk> chunk;
  
    // iterate by pixels,samples and chunks
    for( int pix = range_pix_start; 
         pix <= range_pix_end  &&  i_chunk < n_chunks;  i_chunk++ )
    {
        _ASSERT(i_chunk >=0  && i_chunk < n_chunks);

        chunk = (*chunks)[i_chunk];            
        if(chunk->GetType() & CAlnMap::fSeq)    { // if chunk is aligned

            //calc samples range by chunk 
            const TSignedSeqRange& ch_range = chunk->GetRange();  
            const TSignedSeqRange& ch_aln_range = chunk->GetAlnRange();
            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);
            
            // [sm_start, sm_end] - samples range being rendered
            int sm_start = x_FindSampleToLeft(from);
            int sm_end = x_FindSampleToRight(to + 1);
            sm_start = max(sm_start, 0);
            sm_end = min(sm_end, m_pData->GetSamplesCount() - 1);
            
            // calculate pixels range to render        
            double aln_from = trans.GetAlnPosFromSeqPos(from);
            double aln_to = trans.GetAlnPosFromSeqPos(to);
            if(m_pData->IsNegative())   {
                swap(aln_from, aln_to);
            }
            aln_to += 1;
            
            TVPUnit pix_start = range_pix_start + (TVPUnit) floor((aln_from - left) / scale_x);
            TVPUnit pix_end = range_pix_start + (TVPUnit) ceil((aln_to - left) / scale_x);                        
            pix_start = max(pix_start, range_pix_start);
            pix_end = min(pix_end, range_pix_end);
            
            int sm_inc = m_pData->IsNegative() ? -1 : +1;

            int band_y = top_y;
            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::TValues& values = m_pData->GetValues((CTraceData::EChannel) ch);
                CTraceData::TSignalValue MaxSignal = m_pData->GetMax((CTraceData::EChannel) ch);

                double x1, x2, pix_x1, pix_x2;
                double v1, v2, s, dx, sum, sum_pix_x;

                int sample = m_pData->IsNegative() ? sm_end : sm_start;                
                _ASSERT(sample >= 0);

                for( TVPUnit pix = pix_start; pix <= pix_end;  pix++ )  { // for each pixel
                    // calculate average value for "pix" pixel by integrating values
                    // over the range [pix, pix+1]
                    sum = 0; // integral from values by pix_x
                    sum_pix_x = 0; // length of integrated range

                    x1 = positions[sample];       
                    if(m_pData->IsNegative()) 
                        x1 -= 1.0;
                    pix_x1 = range_pix_start + (trans.GetAlnPosFromSeqPos(x1) - left) / scale_x; 
                    v1 = v2 = values[sample]; //#####

                    if(pix_x1 < pix + 1)    {

                        bool b_next_point = m_pData->IsNegative() ? (sample > sm_start) : (sample < sm_end);
                        if(b_next_point) { // there is second point available
                            x2 = positions[sample + sm_inc];
                            if(m_pData->IsNegative()) 
                                x2 -= 1.0;
                            pix_x2 = range_pix_start + (trans.GetAlnPosFromSeqPos(x2) - left) / scale_x; 
                            v2 = values[sample + sm_inc];
                        
                            if(pix_x1 < pix)    { // fisrt sample is to the left of this pixel                        
                                // replace it fake interpolated sample at x = "pix"
                                v1 += (v2 - v1) * (pix - pix_x1) / (pix_x2 - pix_x1);
                                pix_x1 = pix;
                            }
                        }    

                        while(b_next_point  &&  pix_x2 <= pix + 1) // while second point is inside pixel
                        {
                            dx = pix_x2 - pix_x1;
                            s = 0.5 * (v1 + v2) * dx;
                            _ASSERT(s >=0  &&  dx >= 0);
                            
                            sum += s;
                            sum_pix_x += dx;

                            sample += sm_inc; // advance, x2 becomes x1
                            pix_x1 = pix_x2;
                            v1 = v2;                        
                        
                            b_next_point = m_pData->IsNegative() ? (sample > sm_start) : (sample < sm_end);
                            if(b_next_point) {
                                x2 = positions[sample + sm_inc];
                                if(m_pData->IsNegative()) 
                                    x2 -= 1.0;
                                pix_x2 = range_pix_start + (trans.GetAlnPosFromSeqPos(x2) - left) / scale_x;
                                v2 = values[sample + sm_inc];  
                            } else break;
                        }   
                        _ASSERT(pix_x1 <= pix + 1);
                    
                        if(b_next_point  &&  pix_x2 > pix + 1)    { // second point is outside pixel
                            dx = pix + 1 - pix_x1;
                            _ASSERT(dx >= 0);
                            
                            double v = v1 + (v2 - v1) * dx / (pix_x2 - pix_x1);
                            s = 0.5 * (v1 + v) * dx;
                        
                            _ASSERT(s >=0  &&  dx >= 0);

                            sum += s;
                            sum_pix_x += dx;
                        }
                        double av_v = (sum_pix_x) > 0  ? sum / sum_pix_x : 0;
                    
                        // render pixel
                        double norm = (MaxSignal == 0) ? 0 : (av_v / MaxSignal);
                        const CGlColor& col = GetColorByValue(norm, i_ch);
                        glColorC(col);
                        glVertex2d(pix, band_y);
                        glVertex2d(pix, band_y - band_h - 1);
                    } // if(pix < pix + 1)
                }
                band_y -= band_h;
            }
        }
    }
    glEnd();
}

// returns gradient color corresponding to normalized [0, 1.0] value "value" 
// for channel specified by "signal" 
const CGlColor&     CTraceGraph::GetColorByValue(double value, int signal) const
{
    _ASSERT(value >= 0  &&  value <= 1.0);
    _ASSERT(signal >= 0  &&  signal <= 3);

    int i = (int) (value * kGradColors);
    i = min(i, kGradColors);
    int index = signal * kGradColors + i;
    return m_vSignalColors[index];
}

/// returns index of rightmost sample having m_SeqPos less then "pos".
/// if "pos" is to the left of the trace range function returns -1,
/// if "pos" is to the right of the trace range functions returns "n_samples"
int CTraceGraph::x_FindSampleToLeft(double pos) const
{
    int n_samples = m_pData->GetSamplesCount();        
    if(pos < m_pData->GetSeqFrom()  || n_samples == 0)  {
        return -1;
    } else if(pos > m_pData->GetSeqTo())   {
        return n_samples;
    } else {
        const CTraceData::TPositions& positions = m_pData->GetPositions();
        double scale = ((double) n_samples) / m_pData->GetSeqLength();
        
        // calculate approximate sample index
        int i = (int) (scale * (pos - m_pData->GetSeqFrom()));
        i = min(i, n_samples - 1);
        i = max(i, 0);
        
        if(positions[i] > pos)  {
            for(  ; i > 0  &&  positions[i] > pos;  i-- )  {          
            }
        } else {
            for(  ; ++i < n_samples  &&  positions[i] < pos;  )  {          
            }
            i--;
        }
        return i;
    }
}

/// returns index of the leftmost sample having m_SeqPos greater than "pos"
/// if "pos" is to the left of the trace range function returns -1,
/// if "pos" is to the right of the trace range functions returns "n_samples"
int CTraceGraph::x_FindSampleToRight(double pos) const
{
    int n_samples = m_pData->GetSamplesCount();        
    if(pos < m_pData->GetSeqFrom()  || n_samples == 0)  {
        return -1;
    } else if(pos > m_pData->GetSeqTo())   {
        return n_samples;
    } else {
        const CTraceData::TPositions& positions = m_pData->GetPositions();        
        double scale = ((double) n_samples) / m_pData->GetSeqLength();
        
        // calculate approximate sample index
        int i = (int) (scale * (pos - m_pData->GetSeqFrom()));
        i = min(i, n_samples - 1);
        i = max(i, 0);
        
        if(positions[i] > pos)  {
            for(  ; i > 0  &&  positions[i] > pos;  i-- ) {          
            }
            i++;
        } else {
            for(  ; ++i < n_samples  &&  positions[i] < pos;  ) {          
            }
        }
        return i;
    }
}

END_NCBI_SCOPE

/*
 * ===========================================================================
 * $Log: trace_graph.cpp,v $
 * Revision 1000.1  2004/06/01 21:07:29  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.9
 *
 * Revision 1.9  2004/05/21 22:27:52  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.8  2004/05/14 13:20:51  yazhuk
 * Implemented color reversing for traces on negative strand
 *
 * Revision 1.7  2004/05/10 17:46:35  yazhuk
 * Addressed GCC warnings
 *
 * Revision 1.6  2004/03/29 18:59:47  yazhuk
 * Added support for data loading and Properties management
 *
 * Revision 1.5  2004/03/26 14:59:49  yazhuk
 * Renamed EMode to ESignalStyle; a number of rendering improvements
 *
 * Revision 1.4  2004/03/25 13:05:50  dicuccio
 * Use _TRACE() instead of cout
 *
 * Revision 1.3  2004/03/11 17:50:41  dicuccio
 * Updated typedefs: dropped TDimension, TPosition, TIndex, TColor; use TSeqRange
 * instead of TRange
 *
 * Revision 1.2  2004/03/08 15:38:16  yazhuk
 * Implemented state management expanded/collapsed/hidden
 *
 * Revision 1.1  2004/03/02 15:13:51  yazhuk
 * Initial revision
 *
 * ===========================================================================
 */