生物技术

开发平台：
C/C++

trace_graph.cpp：源码内容
							/*
 * ===========================================================================
 * PRODUCTION $Log: trace_graph.cpp,v $
 * PRODUCTION Revision 1000.1  2004/06/01 21:07:29  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.9
 * PRODUCTION
 * ===========================================================================
 */
/*  $Id: trace_graph.cpp,v 1000.1 2004/06/01 21:07:29 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a 'United States Government Work' under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Authors:  Andrey Yazhuk
 *
 * File Description:
 *
 */
#include <ncbi_pch.hpp>
#include <corelib/ncbistd.hpp>
#include <gui/widgets/aln_multiple/trace_graph.hpp>
#include <gui/opengl/glhelpers.hpp>
#include <gui/types.hpp>
#include <math.h>
BEGIN_NCBI_SCOPE
CTraceGraphProperties::CTraceGraphProperties()
:   m_SignalStyle(eCurve),
    m_ConfGraphState(eExpanded),
    m_SignalGraphState(eExpanded),
    m_bReverseColors(true)
{
}
CTraceGraphProperties&  
    CTraceGraphProperties::operator=(const CTraceGraphProperties& orig)
{
    m_SignalStyle = orig.m_SignalStyle;
    m_ConfGraphState = orig.m_ConfGraphState;
    m_SignalGraphState = orig.m_SignalGraphState;
    return *this;
}
////////////////////////////////////////////////////////////////////////////////
/// CTraceGraph
const static int kGradColors = 32;
CTraceGraph::CTraceGraph(const CBioseq_Handle& handle, bool b_neg_strand)
:   m_DataProxy(handle, b_neg_strand),
    m_pData(NULL)
{
}
CTraceGraph::~CTraceGraph()
{
    Destroy();
}
bool    CTraceGraph::HasData() const
{
    return m_DataProxy.HasData();
}
bool    CTraceGraph::IsCreated() const
{
    return (m_pData != NULL);
}
bool    CTraceGraph::Create()
{
    m_pData = m_DataProxy.LoadData();
    if(m_pData) {
        SetConfGraphState(CTraceGraphProperties::eCollapsed);
    
        bool b_ch = m_pData->GetSamplesCount() > 0;
        SetSignalGraphState(b_ch    ? CTraceGraphProperties::eExpanded 
                                    : CTraceGraphProperties::eHidden);
        
        m_vSignalColors.resize(4 * kGradColors);
        double k = 1.0 / kGradColors;
    
        for( int j = 0; j < kGradColors; j++ )  {
            float v = 1.0f - k * j;
            m_vSignalColors[j] = CGlColor(1.0f, v, v); //red
            m_vSignalColors[kGradColors + j] = CGlColor(v, 1.0f, v); //green
            m_vSignalColors[2 * kGradColors + j] = CGlColor(v, v, 1.0f); //blue
            m_vSignalColors[3 * kGradColors + j] = CGlColor(0.5f * (1 + v), v, 0.5f * (1 + v)); //purple        
        }
        m_pData->CalculateMax();
        return true;
    } else return false;
}
void    CTraceGraph::Destroy()
{
    if(m_pData) {
        delete m_pData;
        m_pData = NULL;
        m_vSignalColors.clear();
    }
}
const IAlnRowGraphProperties*     CTraceGraph::GetProperties() const
{
    return &m_Props;
}
void    CTraceGraph::SetProperties(IAlnRowGraphProperties* props)
{
    CTraceGraphProperties* trace_props = 
        dynamic_cast<CTraceGraphProperties*>(props);
    if(trace_props) {
        m_Props = *trace_props;
    }
}
void    CTraceGraph::SetConfGraphState(EGraphState state)
{
    m_Props.m_ConfGraphState = state;
}
void    CTraceGraph::SetSignalGraphState(EGraphState state)
{
    m_Props.m_SignalGraphState = state;
}
// vertical spacing between graph area border and graph
static const int kGraphOffsetY = 1;
static const int kConfGraphPrefH = 24;
static const int kSignalGraphPrefH = 40;
static const int kCollapsedGraphH = 11;
int     CTraceGraph::GetPreferredHeight()    const
{
    int h = x_GetConfGraphH();
    h += x_GetSignalGraphH();
    h += 2;    
    return h;
}
int     CTraceGraph::x_GetConfGraphH() const
{
    switch(m_Props.m_ConfGraphState)    {
    case CTraceGraphProperties::eCollapsed: return kCollapsedGraphH; 
    case CTraceGraphProperties::eExpanded: return kConfGraphPrefH;
    case CTraceGraphProperties::eHidden: return 0;
    }
    return 0;
}
int     CTraceGraph::x_GetSignalGraphH() const
{
    switch(m_Props.m_SignalGraphState)    {
    case CTraceGraphProperties::eCollapsed: return kCollapsedGraphH; 
    case CTraceGraphProperties::eExpanded: return kSignalGraphPrefH;
    case CTraceGraphProperties::eHidden: return 0;
    }    
    return 0;
}
// renders both Confidence graph and Chromatograms
void    CTraceGraph::Render(CGlPane& pane, int y, int h, const TChunkVec& chunks)
{
    CGlAttrGuard AttrGuard(GL_ENABLE_BIT  | GL_COLOR_BUFFER_BIT | GL_HINT_BIT 
                            | GL_LINE_SMOOTH | GL_POLYGON_MODE | GL_LINE_BIT);
         
    int conf_h = x_GetConfGraphH();
    int sig_h = x_GetSignalGraphH();
    int top_y = y;
    
    pane.OpenOrtho();
    x_RenderContour(pane, top_y, sig_h, conf_h + sig_h, chunks);  // render background
        
    if(m_Props.m_SignalGraphState == CTraceGraphProperties::eExpanded)    {
        if(pane.GetScaleX() < 1.0)   { // render signal graph
            if(m_Props.m_SignalStyle == CTraceGraphProperties::eCurve)   { // render curves
                glEnable(GL_BLEND);
                glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    
                glEnable(GL_LINE_SMOOTH);
                glHint(GL_LINE_SMOOTH_HINT, GL_NICEST);
                glLineWidth(0.5);
                x_RenderSignalGraph(pane, top_y, sig_h, chunks);
            } else {            // render intensity bands
                pane.Close();
                pane.OpenPixels();
                glDisable(GL_BLEND);
                glDisable(GL_LINE_SMOOTH);
                glLineWidth(1.0);
                x_RenderIntensityGraphs(pane, top_y, sig_h, chunks);                            
            
                pane.Close();
                pane.OpenOrtho();
            }        
        }
        top_y += sig_h;
    }
        
    if(m_Props.m_ConfGraphState == CTraceGraphProperties::eExpanded)   {
        // render confidence graph
        x_RenderConfGraph(pane, top_y, conf_h, chunks);        
    }
    
    pane.Close();    
}
// Functor for coordinate translation align <->seq
struct SChunkTranslator
{
public:
    SChunkTranslator()
    : m_pSeqRange(NULL), m_pAlnRange(NULL), m_bNegative(false)
    {
    }
    
    void Init(const TSignedSeqRange& seq_range, 
              const TSignedSeqRange& aln_range, bool negative)
    {        
        m_pSeqRange = &seq_range;
        m_pAlnRange = &aln_range;
        m_bNegative = negative;
        _ASSERT(m_pSeqRange->GetLength() == m_pAlnRange->GetLength());    
    }
    inline TSignedSeqPos GetAlnPosFromSeqPos(TSignedSeqPos seq_pos)
    {
        _ASSERT(m_pSeqRange  && m_pAlnRange);
        TSignedSeqPos seq_off = seq_pos - m_pSeqRange->GetFrom();
        if(m_bNegative)    {
          return m_pAlnRange->GetTo() - seq_off;
        } else return seq_off + m_pAlnRange->GetFrom();
    }
    inline double GetAlnPosFromSeqPos(double seq_pos)
    {
        _ASSERT(m_pSeqRange  && m_pAlnRange);
        double seq_off = seq_pos - m_pSeqRange->GetFrom();
        if(m_bNegative)    {
          return m_pAlnRange->GetTo() - seq_off;
        } else return seq_off + m_pAlnRange->GetFrom();
    }
    inline TSignedSeqPos GetSeqPosFromAlnPos(TSignedSeqPos aln_pos)
    {
        _ASSERT(m_pSeqRange  && m_pAlnRange);
        TSignedSeqPos aln_off = aln_pos - m_pAlnRange->GetFrom();
        if(m_bNegative)    {
          return m_pSeqRange->GetTo() - aln_off;
        } else return aln_off + m_pSeqRange->GetFrom();
    }
    inline double GetSeqPosFromAlnPos(double aln_pos)
    {
        _ASSERT(m_pSeqRange  && m_pAlnRange);
        double aln_off = aln_pos - m_pAlnRange->GetFrom();
        if(m_bNegative)    {
          return m_pSeqRange->GetTo() - aln_off;
        } else return aln_off + m_pSeqRange->GetFrom();
    }
protected:
    const TSignedSeqRange* m_pSeqRange;
    const TSignedSeqRange* m_pAlnRange;
    bool    m_bNegative;
};
void    CTraceGraph::x_RenderContour(CGlPane& pane, int y, int top_h, int total_h,
                                      const TChunkVec& chunks)
{
    glDisable(GL_BLEND);
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
    glColor3d(0.9, 0.9, 0.9);
    
    TModelUnit offset_x = pane.GetOffsetX();
    TModelUnit y1 = y;
    TModelUnit y2 = y + top_h - 1;
    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;
    
                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
 *
 * ===========================================================================
 */

						