rawVolModel.cpp

/* Copyright (c) 2007       Maxim Makhinya
 *
 * This library is free software; you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License version 2.1 as published
 * by the Free Software Foundation.
 *  
 * This library 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 Lesser General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */

#include "rawVolModel.h"
#include "hlp.h"

namespace eVolve
{


using hlpFuncs::clip;
using hlpFuncs::hFile;


static GLuint createPreintegrationTable( const uint8_t* Table );

static bool readTransferFunction( FILE* file, std::vector<uint8_t>& TF );

static bool readDimensionsAndScaling
( 
    FILE* file, 
    uint32_t& w, uint32_t& h, uint32_t& d, 
    VolumeScaling& volScaling 
);


// Read volume dimensions, scaling and transfer function
RawVolumeModel::RawVolumeModel( const std::string& filename  )
        : _headerLoaded( false )
        , _filename( filename )
        , _preintName  ( 0 )
        , _glewContext( 0 )
{}

bool RawVolumeModel::loadHeader( const float brightness, const float alpha )
{
    EQASSERT( !_headerLoaded );
    EQASSERT( brightness > 0.0f );
    EQASSERT( alpha > 0.0f );
    EQLOG( eq::LOG_CUSTOM ) << "-------------------------------------------"
                            << std::endl << "model: " << _filename;

    hFile header( fopen( ( _filename + std::string( ".vhf" ) ).c_str(), "rb" ));

    if( header.f==0 )
    {
        EQERROR << "Can't open header file" << std::endl;
        return false;
    }

    if( !readDimensionsAndScaling( header.f, _w, _h, _d, _volScaling ) )
        return false;

    _resolution = EQ_MAX( _w, EQ_MAX( _h, _d ) );

    if( !readTransferFunction( header.f, _TF ))
        return false;

    _headerLoaded = true;

    if( brightness != 1.0f )
    {
        for( size_t i = 0; i < _TF.size(); i+=4 )
        {
            _TF[i+0] = static_cast< uint8_t >( _TF[i+0] * brightness );
            _TF[i+1] = static_cast< uint8_t >( _TF[i+1] * brightness );
            _TF[i+2] = static_cast< uint8_t >( _TF[i+2] * brightness );
        }
    }

    if( alpha != 1.0f )
        for( size_t i = 3; i < _TF.size(); i+=4 )
            _TF[i] = static_cast< uint8_t >( _TF[i] * alpha );

    return true;
}


static int32_t calcHashKey( const eq::Range& range )
{
    return static_cast<int32_t>(( range.start*10000.f + range.end )*10000.f );
}


bool RawVolumeModel::getVolumeInfo( VolumeInfo& info, const eq::Range& range )
{
    if( !_headerLoaded && !loadHeader( 1.0f, 1.0f ))
        return false;

    if( _preintName == 0 )
    {
        EQLOG( eq::LOG_CUSTOM ) << "Creating preint" << std::endl;
        _preintName = createPreintegrationTable( &_TF[0] );
    }

          VolumePart* volumePart = 0;
    const int32_t     key        = calcHashKey( range );

    if( _volumeHash.find( key ) == _volumeHash.end( ) )
    {
        // new key
        volumePart = &_volumeHash[ key ];
        if( !_createVolumeTexture( volumePart->volume, volumePart->TD, range ))
            return false;
    }
    else
    {   // old key
        volumePart = &_volumeHash[ key ];
    }

    info.volume     = volumePart->volume;
    info.TD         = volumePart->TD;
    info.preint     = _preintName;
    info.volScaling  = _volScaling;

    return true;
}


void RawVolumeModel::releaseVolumeInfo( const eq::Range& range )
{
    const int32_t key = calcHashKey( range );
    if( _volumeHash.find( key ) == _volumeHash.end() )
        return;

    _volumeHash.erase( key );
}


static uint32_t calcMinPow2( uint32_t size )
{
    if( size == 0 )
        return 0;
    
    size--;
    uint32_t res = 1;

    while( size > 0 )
    {
        res  <<= 1;
        size >>= 1;
    }
    return res;
}


bool RawVolumeModel::_createVolumeTexture(        GLuint&    volume,
                                                  DataInTextureDimensions& TD,
                                            const eq::Range& range    )
{
    const uint32_t w = _w;
    const uint32_t h = _h;
    const uint32_t d = _d;

    const int32_t bwStart = 2; //border width from left
    const int32_t bwEnd   = 2; //border width from right

    const int32_t s =
            clip<int32_t>( static_cast< int32_t >( d*range.start ), 0, d-1 );

    const int32_t e =
            clip<int32_t>( static_cast< int32_t >( d*range.end-1 ), 0, d-1 );

    const uint32_t start =
                static_cast<uint32_t>( clip<int32_t>( s-bwStart, 0, d-1 ) );

    const uint32_t end   =
                static_cast<uint32_t>( clip<int32_t>( e+bwEnd  , 0, d-1 ) );

    const uint32_t depth = end-start+1;

    const uint32_t tW = calcMinPow2( w );
    const uint32_t tH = calcMinPow2( h );
    const uint32_t tD = calcMinPow2( depth );

    //texture scaling coefficients
    TD.W  = static_cast<float>( w     ) / static_cast<float>( tW );
    TD.H  = static_cast<float>( h     ) / static_cast<float>( tH );
    TD.D  = static_cast<float>( e-s+1 ) / static_cast<float>( tD );
    TD.D /= range.end>range.start ? (range.end-range.start) : 1.0 ;

    // Shift coefficient and left border in texture for depth
    TD.Do = range.start;
    TD.Db = range.start > 0.0001 ? bwStart / static_cast<float>(tD) : 0;

    EQLOG( eq::LOG_CUSTOM )
            << "==============================================="   << std::endl
            << " w: "  << w << " " << tW
            << " h: "  << h << " " << tH
            << " d: "  << d << " " << depth << " " << tD           << std::endl
            << " r: "  << _resolution                              << std::endl
            << " ws: " << TD.W  << " hs: " << TD.H  << " wd: " << TD.D 
            << " Do: " << TD.Do << " Db: " << TD.Db                << std::endl
            << " s= "  << start << " e= "  << end                  << std::endl;

    // Reading of requested part of a volume
    std::vector<uint8_t> data( tW*tH*tD*4, 0 );
    const uint32_t  wh4 =  w *  h * 4;
    const uint32_t tWH4 = tW * tH * 4;

    std::ifstream file ( _filename.c_str(), std::ifstream::in |
                         std::ifstream::binary | std::ifstream::ate );

    if( !file.is_open() )
    {
        EQERROR << "Can't open model data file";
        return false;
    }

    file.seekg( wh4*start, std::ios::beg );

    if( w==tW && h==tH ) // width and height are power of 2
    {
        file.read( (char*)( &data[0] ), wh4*depth );
    }
    else if( w==tW )     // only width is power of 2
    {
        for( uint32_t i=0; i<depth; i++ )
            file.read( (char*)( &data[i*tWH4] ), wh4 );
    }
    else
    {               // nor width nor heigh is power of 2
        const uint32_t   w4 =  w * 4;
        const uint32_t  tW4 = tW * 4;
        
        for( uint32_t i=0; i<depth; i++ )
            for( uint32_t j=0; j<h; j++ )
                file.read( (char*)( &data[ i*tWH4 + j*tW4] ), w4 );
    }

    file.close();

    EQASSERT( _glewContext );
    // create 3D texture
    glGenTextures( 1, &volume );
    EQLOG( eq::LOG_CUSTOM ) << "generated texture: " << volume << std::endl;
    glBindTexture(GL_TEXTURE_3D, volume);

    glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_WRAP_S    , GL_CLAMP_TO_EDGE );
    glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_WRAP_T    , GL_CLAMP_TO_EDGE );
    glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_WRAP_R    , GL_CLAMP_TO_EDGE );
    glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR        );
    glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR        );

    glTexImage3D(   GL_TEXTURE_3D, 
                    0, GL_RGBA, tW, tH, tD, 
                    0, GL_RGBA, GL_UNSIGNED_BYTE, (GLvoid*)(&data[0]) );
    return true;
}


static void normalizeScaling
(
    const uint32_t       w,
    const uint32_t       h,
    const uint32_t       d,
          VolumeScaling& scaling
)
{
//Correct proportions according to real size of volume
    float maxS = EQ_MAX( w, EQ_MAX( h, d ) );

    scaling.W *= w / maxS;
    scaling.H *= h / maxS;
    scaling.D *= d / maxS;

//Make maximum proportion equal to 1.0
    maxS = EQ_MAX( scaling.W, EQ_MAX( scaling.H, scaling.D ) );

    scaling.W /= maxS;
    scaling.H /= maxS;
    scaling.D /= maxS;
}


static bool readDimensionsAndScaling
( 
    FILE* file, 
    uint32_t& w, uint32_t& h, uint32_t& d, 
    VolumeScaling& volScaling 
)
{
        fscanf( file, "w=%u\n", &w );
        fscanf( file, "h=%u\n", &h );
    if( fscanf( file, "d=%u\n", &d ) != 1 )
    {
        EQERROR << "Can't read dimensions from header file" << std::endl;
        return false;
    }

        fscanf( file, "wScale=%g\n", &volScaling.W );
        fscanf( file, "hScale=%g\n", &volScaling.H );
    if( fscanf( file, "dScale=%g\n", &volScaling.D ) != 1 )
    {
        EQERROR << "Can't read scaling from header file: " << std::endl;
        return false;
    }

    if( w<1 || h<1 || d<1 ||
        volScaling.W<0.001 || 
        volScaling.H<0.001 || 
        volScaling.W<0.001    )
    {
        EQERROR << "volume scaling is incorrect, check header file"<< std::endl;
        return false;
    }

    normalizeScaling( w, h, d, volScaling );

    EQLOG( eq::LOG_CUSTOM ) << " "  << w << "x" << h << "x" << d << std::endl
                            << "( " << volScaling.W << " x "
                                    << volScaling.H << " x "
                                    << volScaling.D << " )"       << std::endl;
    return true;
}


static bool readTransferFunction( FILE* file,  std::vector<uint8_t>& TF )
{
    if( fscanf(file,"TF:\n") !=0 )
    {
        EQERROR << "Error in header file, can't find 'TF:' marker.";
        return false;
    }

    uint32_t TFSize;
    fscanf( file, "size=%u\n", &TFSize );

    if( TFSize!=256  )
        EQWARN << "Wrong size of transfer function, should be 256" << std::endl;

    TFSize = clip<int32_t>( TFSize, 1, 256 );
    TF.resize( TFSize*4 );

    int tmp;
    for( uint32_t i=0; i<TFSize; i++ )
    {
            fscanf( file, "r=%d\n", &tmp ); TF[4*i  ] = tmp;
            fscanf( file, "g=%d\n", &tmp ); TF[4*i+1] = tmp;
            fscanf( file, "b=%d\n", &tmp ); TF[4*i+2] = tmp;
        if( fscanf( file, "a=%d\n", &tmp ) != 1 )
        {
            EQERROR << "Failed to read entity #" << i 
                    << " of TF from header file" << std::endl;
            return i;
        }
        TF[4*i+3] = tmp;
    }

    return true;
}


static GLuint createPreintegrationTable( const uint8_t *Table )
{
    EQLOG( eq::LOG_CUSTOM ) << "Calculating preintegration table" << std::endl;

    double rInt[256]; rInt[0] = 0.;
    double gInt[256]; gInt[0] = 0.;
    double bInt[256]; bInt[0] = 0.;
    double aInt[256]; aInt[0] = 0.;

    for( int i=1; i<256; i++ )
    {
        const double tauc = ( Table[(i-1)*4+3] + Table[i*4+3] ) / 2.;

        rInt[i] = rInt[i-1] + ( Table[(i-1)*4+0] + Table[i*4+0] )/2.*tauc/255.;
        gInt[i] = gInt[i-1] + ( Table[(i-1)*4+1] + Table[i*4+1] )/2.*tauc/255.;
        bInt[i] = bInt[i-1] + ( Table[(i-1)*4+2] + Table[i*4+2] )/2.*tauc/255.;
        aInt[i] = aInt[i-1] + tauc;
    }

    std::vector<GLubyte> lookupImg( 256*256*4, 255 ); // Preint Texture

    int lookupindex=0;

    for( int sb=0; sb<256; sb++ )
    {
        for( int sf=0; sf<256; sf++ )
        {
            int smin, smax;
            if( sb<sf ) { smin=sb; smax=sf; }
            else        { smin=sf; smax=sb; }

            int rcol, gcol, bcol, acol;
            if( smax != smin )
            {
                const double factor = 1. / (double)(smax-smin);
                rcol = static_cast<int>( (rInt[smax]-rInt[smin])*factor );
                gcol = static_cast<int>( (gInt[smax]-gInt[smin])*factor );
                bcol = static_cast<int>( (bInt[smax]-bInt[smin])*factor );
                acol = static_cast<int>( 
                        256.*(    exp(-(aInt[smax]-aInt[smin])*factor/255.)));
            } else
            {
                const int    index  = smin*4;
                const double factor = 1./255.;
                rcol = static_cast<int>( Table[index+0]*Table[index+3]*factor );
                gcol = static_cast<int>( Table[index+1]*Table[index+3]*factor );
                bcol = static_cast<int>( Table[index+2]*Table[index+3]*factor );
                acol = static_cast<int>( 256.*(    exp(-Table[index+3]/255.)) );
            }
            lookupImg[lookupindex++] = clip( rcol, 0, 255 );//MIN( rcol, 255 );
            lookupImg[lookupindex++] = clip( gcol, 0, 255 );//MIN( gcol, 255 );
            lookupImg[lookupindex++] = clip( bcol, 0, 255 );//MIN( bcol, 255 );
            lookupImg[lookupindex++] = clip( acol, 0, 255 );//MIN( acol, 255 );
        }
    }

    GLuint preintName = 0;
    glGenTextures( 1, &preintName );
    glBindTexture( GL_TEXTURE_2D, preintName );
    glTexImage2D(  GL_TEXTURE_2D, 0, GL_RGBA, 256, 256, 0, 
                   GL_RGBA, GL_UNSIGNED_BYTE, &lookupImg[0] );

    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S,     GL_CLAMP_TO_EDGE );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR        );
    glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR        );

    return preintName;
}


}