Equalizer  1.2.1
rawVolModel.cpp
00001 
00002 /* Copyright (c) 2007-2011, Maxim Makhinya  <maxmah@gmail.com>
00003  *
00004  * Redistribution and use in source and binary forms, with or without
00005  * modification, are permitted provided that the following conditions are met:
00006  *
00007  * - Redistributions of source code must retain the above copyright notice, this
00008  *   list of conditions and the following disclaimer.
00009  * - Redistributions in binary form must reproduce the above copyright notice,
00010  *   this list of conditions and the following disclaimer in the documentation
00011  *   and/or other materials provided with the distribution.
00012  * - Neither the name of Eyescale Software GmbH nor the names of its
00013  *   contributors may be used to endorse or promote products derived from this
00014  *   software without specific prior written permission.
00015  *
00016  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00017  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00018  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00019  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
00020  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00021  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
00022  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
00023  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
00024  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
00025  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00026  * POSSIBILITY OF SUCH DAMAGE.
00027  */
00028 
00029 #include "rawVolModel.h"
00030 #include "hlp.h"
00031 
00032 namespace eVolve
00033 {
00034 
00035 
00036 using hlpFuncs::clip;
00037 using hlpFuncs::hFile;
00038 
00039 
00040 static GLuint createPreintegrationTable( const uint8_t* Table );
00041 
00042 static bool readTransferFunction( FILE* file, std::vector<uint8_t>& TF );
00043 
00044 static bool readDimensionsAndScaling
00045 ( 
00046     FILE* file, 
00047     uint32_t& w, uint32_t& h, uint32_t& d, 
00048     VolumeScaling& volScaling 
00049 );
00050 
00051 
00052 // Read volume dimensions, scaling and transfer function
00053 RawVolumeModel::RawVolumeModel( const std::string& filename  )
00054         : _headerLoaded( false )
00055         , _filename( filename )
00056         , _preintName  ( 0 )
00057         , _w( 0 )
00058         , _h( 0 )
00059         , _d( 0 )
00060         , _tW( 0 )
00061         , _tH( 0 )
00062         , _tD( 0 )
00063         , _hasDerivatives( true )
00064         , _glewContext( 0 )
00065 {}
00066 
00067 bool RawVolumeModel::loadHeader( const float brightness, const float alpha )
00068 {
00069     EQASSERT( !_headerLoaded );
00070     EQASSERT( brightness > 0.0f );
00071     EQASSERT( alpha > 0.0f );
00072     EQLOG( eq::LOG_CUSTOM ) << "-------------------------------------------"
00073                             << std::endl << "model: " << _filename;
00074 
00075     hFile header( fopen( ( _filename + std::string( ".vhf" ) ).c_str(), "rb" ));
00076 
00077     if( header.f==0 )
00078     {
00079         EQERROR << "Can't open header file" << std::endl;
00080         return false;
00081     }
00082 
00083     // test for raw+der or raw
00084     const size_t fNameLen = _filename.length();
00085     _hasDerivatives = 
00086         ( fNameLen >= 6 && _filename.substr( fNameLen-6, 6 ) == "_d.raw" ) ?
00087         true : false;
00088 
00089     if( !readDimensionsAndScaling( header.f, _w, _h, _d, _volScaling ) )
00090         return false;
00091 
00092     _resolution = EQ_MAX( _w, EQ_MAX( _h, _d ) );
00093 
00094     if( !readTransferFunction( header.f, _TF ))
00095         return false;
00096 
00097     _headerLoaded = true;
00098 
00099     if( brightness != 1.0f )
00100     {
00101         for( size_t i = 0; i < _TF.size(); i+=4 )
00102         {
00103             _TF[i+0] = static_cast< uint8_t >( _TF[i+0] * brightness );
00104             _TF[i+1] = static_cast< uint8_t >( _TF[i+1] * brightness );
00105             _TF[i+2] = static_cast< uint8_t >( _TF[i+2] * brightness );
00106         }
00107     }
00108 
00109     if( alpha != 1.0f )
00110         for( size_t i = 3; i < _TF.size(); i+=4 )
00111             _TF[i] = static_cast< uint8_t >( _TF[i] * alpha );
00112 
00113     return true;
00114 }
00115 
00116 
00117 static int32_t calcHashKey( const eq::Range& range )
00118 {
00119     return static_cast<int32_t>(( range.start*10000.f + range.end )*10000.f );
00120 }
00121 
00122 
00123 bool RawVolumeModel::getVolumeInfo( VolumeInfo& info, const eq::Range& range )
00124 {
00125     if( !_headerLoaded && !loadHeader( 1.0f, 1.0f ))
00126         return false;
00127 
00128     if( _preintName == 0 )
00129     {
00130         EQLOG( eq::LOG_CUSTOM ) << "Creating preint" << std::endl;
00131         _preintName = createPreintegrationTable( &_TF[0] );
00132     }
00133 
00134           VolumePart* volumePart = 0;
00135     const int32_t     key        = calcHashKey( range );
00136 
00137     if( _volumeHash.find( key ) == _volumeHash.end( ) )
00138     {
00139         // new key
00140         volumePart = &_volumeHash[ key ];
00141         if( !_createVolumeTexture( volumePart->volume, volumePart->TD, range ))
00142             return false;
00143     }
00144     else
00145     {   // old key
00146         volumePart = &_volumeHash[ key ];
00147     }
00148 
00149     info.volume     = volumePart->volume;
00150     info.TD         = volumePart->TD;
00151     info.preint     = _preintName;
00152     info.volScaling = _volScaling;
00153     info.hasDerivatives = _hasDerivatives;
00154     if( _hasDerivatives )
00155     {
00156         info.voxelSize.W  = 1.f;
00157         info.voxelSize.H  = 1.f;
00158         info.voxelSize.D  = 1.f;
00159     }else
00160     {
00161         info.voxelSize.W  = 1.f / _tW;
00162         info.voxelSize.H  = 1.f / _tH;
00163         info.voxelSize.D  = 1.f / _tD;
00164     }
00165     return true;
00166 }
00167 
00168 
00169 void RawVolumeModel::releaseVolumeInfo( const eq::Range& range )
00170 {
00171     const int32_t key = calcHashKey( range );
00172     if( _volumeHash.find( key ) == _volumeHash.end() )
00173         return;
00174 
00175     _volumeHash.erase( key );
00176 }
00177 
00178 
00181 static uint32_t calcMinPow2( uint32_t size )
00182 {
00183     if( size == 0 )
00184         return 0;
00185     
00186     size--;
00187     uint32_t res = 1;
00188 
00189     while( size > 0 )
00190     {
00191         res  <<= 1;
00192         size >>= 1;
00193     }
00194     return res;
00195 }
00196 
00197 
00200 bool RawVolumeModel::_createVolumeTexture(        GLuint&    volume,
00201                                                   DataInTextureDimensions& TD,
00202                                             const eq::Range& range    )
00203 {
00204     const uint32_t w = _w;
00205     const uint32_t h = _h;
00206     const uint32_t d = _d;
00207     const uint32_t bytes = _hasDerivatives ? 4 : 1;
00208 
00209     const int32_t bwStart = 2; //border width from left
00210     const int32_t bwEnd   = 2; //border width from right
00211 
00212     const int32_t s =
00213             clip<int32_t>( static_cast< int32_t >( d*range.start ), 0, d-1 );
00214 
00215     const int32_t e =
00216             clip<int32_t>( static_cast< int32_t >( d*range.end-1 ), 0, d-1 );
00217 
00218     const uint32_t start =
00219                 static_cast<uint32_t>( clip<int32_t>( s-bwStart, 0, d-1 ) );
00220 
00221     const uint32_t end   =
00222                 static_cast<uint32_t>( clip<int32_t>( e+bwEnd  , 0, d-1 ) );
00223 
00224     const uint32_t depth = end-start+1;
00225 
00226     _tW = calcMinPow2( w );
00227     _tH = calcMinPow2( h );
00228     _tD = calcMinPow2( depth );
00229 
00230     //texture scaling coefficients
00231     TD.W  = static_cast<float>( w     ) / static_cast<float>( _tW );
00232     TD.H  = static_cast<float>( h     ) / static_cast<float>( _tH );
00233     TD.D  = static_cast<float>( e-s+1 ) / static_cast<float>( _tD );
00234     TD.D /= range.end>range.start ? (range.end-range.start) : 1.0f;
00235 
00236     // Shift coefficient and left border in texture for depth
00237     TD.Do = range.start;
00238     TD.Db = range.start > 0.0001 ? bwStart / static_cast<float>(_tD) : 0;
00239 
00240     EQLOG( eq::LOG_CUSTOM )
00241             << "==============================================="   << std::endl
00242             << " w: "  << w << " " << _tW
00243             << " h: "  << h << " " << _tH
00244             << " d: "  << d << " " << depth << " " << _tD           << std::endl
00245             << " r: "  << _resolution                              << std::endl
00246             << " ws: " << TD.W  << " hs: " << TD.H  << " wd: " << TD.D 
00247             << " Do: " << TD.Do << " Db: " << TD.Db                << std::endl
00248             << " s= "  << start << " e= "  << end                  << std::endl;
00249 
00250     // Reading of requested part of a volume
00251     std::vector<uint8_t> data( _tW*_tH*_tD*bytes, 0 );
00252     const uint32_t  wh4 =   w *   h * bytes;
00253     const uint32_t tWH4 = _tW * _tH * bytes;
00254 
00255     std::ifstream file ( _filename.c_str(), std::ifstream::in |
00256                          std::ifstream::binary | std::ifstream::ate );
00257 
00258     if( !file.is_open() )
00259     {
00260         EQERROR << "Can't open model data file";
00261         return false;
00262     }
00263 
00264     file.seekg( wh4*start, std::ios::beg );
00265 
00266     if( w==_tW && h==_tH ) // width and height are power of 2
00267     {
00268         file.read( (char*)( &data[0] ), wh4*depth );
00269     }
00270     else if( w==_tW )     // only width is power of 2
00271     {
00272         for( uint32_t i=0; i<depth; i++ )
00273             file.read( (char*)( &data[i*tWH4] ), wh4 );
00274     }
00275     else
00276     {               // nor width nor heigh is power of 2
00277         const uint32_t   w4 =   w * bytes;
00278         const uint32_t  tW4 = _tW * bytes;
00279 
00280         for( uint32_t i=0; i<depth; i++ )
00281             for( uint32_t j=0; j<h; j++ )
00282                 file.read( (char*)( &data[ i*tWH4 + j*tW4] ), w4 );
00283     }
00284 
00285     file.close();
00286 
00287     EQASSERT( _glewContext );
00288     // create 3D texture
00289     glGenTextures( 1, &volume );
00290     EQLOG( eq::LOG_CUSTOM ) << "generated texture: " << volume << std::endl;
00291     glBindTexture(GL_TEXTURE_3D, volume);
00292 
00293     glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_WRAP_S    , GL_CLAMP_TO_EDGE );
00294     glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_WRAP_T    , GL_CLAMP_TO_EDGE );
00295     glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_WRAP_R    , GL_CLAMP_TO_EDGE );
00296     glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR        );
00297     glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR        );
00298 
00299     if( _hasDerivatives )
00300     {
00301         glTexImage3D(   GL_TEXTURE_3D,
00302                         0, GL_RGBA, _tW, _tH, _tD,
00303                         0, GL_RGBA, GL_UNSIGNED_BYTE, (GLvoid*)(&data[0]) );
00304     }else
00305     {
00306         glTexImage3D(   GL_TEXTURE_3D,
00307                         0, GL_ALPHA, _tW, _tH, _tD,
00308                         0, GL_ALPHA, GL_UNSIGNED_BYTE, (GLvoid*)(&data[0]) );
00309     }
00310 
00311     return true;
00312 }
00313 
00314 
00319 static void normalizeScaling
00320 (
00321     const uint32_t       w,
00322     const uint32_t       h,
00323     const uint32_t       d,
00324           VolumeScaling& scaling
00325 )
00326 {
00327 //Correct proportions according to real size of volume
00328     float maxS = float( EQ_MAX( w, EQ_MAX( h, d ) ));
00329 
00330     scaling.W *= w / maxS;
00331     scaling.H *= h / maxS;
00332     scaling.D *= d / maxS;
00333 
00334 //Make maximum proportion equal to 1.0
00335     maxS = EQ_MAX( scaling.W, EQ_MAX( scaling.H, scaling.D ) );
00336 
00337     scaling.W /= maxS;
00338     scaling.H /= maxS;
00339     scaling.D /= maxS;
00340 }
00341 
00342 
00343 static bool readDimensionsAndScaling
00344 ( 
00345     FILE* file, 
00346     uint32_t& w, uint32_t& h, uint32_t& d, 
00347     VolumeScaling& volScaling 
00348 )
00349 {
00350     if( fscanf( file, "w=%u\n", &w ) == EOF )
00351         return false;
00352     if( fscanf( file, "h=%u\n", &h ) == EOF )
00353         return false;
00354     if( fscanf( file, "d=%u\n", &d ) != 1 )
00355     {
00356         EQERROR << "Can't read dimensions from header file" << std::endl;
00357         return false;
00358     }
00359 
00360     if( fscanf( file, "wScale=%g\n", &volScaling.W ) == EOF )
00361         return false;
00362     if( fscanf( file, "hScale=%g\n", &volScaling.H ) == EOF )
00363         return false;
00364     if( fscanf( file, "dScale=%g\n", &volScaling.D ) != 1 )
00365     {
00366         EQERROR << "Can't read scaling from header file: " << std::endl;
00367         return false;
00368     }
00369 
00370     if( w<1 || h<1 || d<1 ||
00371         volScaling.W<0.001 || 
00372         volScaling.H<0.001 || 
00373         volScaling.W<0.001    )
00374     {
00375         EQERROR << "volume scaling is incorrect, check header file"<< std::endl;
00376         return false;
00377     }
00378 
00379     normalizeScaling( w, h, d, volScaling );
00380 
00381     EQLOG( eq::LOG_CUSTOM ) << " "  << w << "x" << h << "x" << d << std::endl
00382                             << "( " << volScaling.W << " x "
00383                                     << volScaling.H << " x "
00384                                     << volScaling.D << " )"       << std::endl;
00385     return true;
00386 }
00387 
00388 
00389 static bool readTransferFunction( FILE* file,  std::vector<uint8_t>& TF )
00390 {
00391     if( fscanf(file,"TF:\n") !=0 )
00392     {
00393         EQERROR << "Error in header file, can't find 'TF:' marker.";
00394         return false;
00395     }
00396 
00397     uint32_t TFSize;
00398     if( fscanf( file, "size=%u\n", &TFSize ) == EOF )
00399         return false;
00400 
00401     if( TFSize!=256  )
00402         EQWARN << "Wrong size of transfer function, should be 256" << std::endl;
00403 
00404     TFSize = clip<int32_t>( TFSize, 1, 256 );
00405     TF.resize( TFSize*4 );
00406 
00407     int tmp;
00408     for( uint32_t i=0; i<TFSize; i++ )
00409     {
00410         if( fscanf( file, "r=%d\n", &tmp ) == EOF )
00411             return false;
00412         TF[4*i  ] = tmp;
00413         if( fscanf( file, "g=%d\n", &tmp ) == EOF )
00414             return false;
00415         TF[4*i+1] = tmp;
00416         if( fscanf( file, "b=%d\n", &tmp ) == EOF )
00417             return false;
00418         TF[4*i+2] = tmp;
00419         if( fscanf( file, "a=%d\n", &tmp ) != 1 )
00420         {
00421             EQERROR << "Failed to read entity #" << i 
00422                     << " of TF from header file" << std::endl;
00423             return i;
00424         }
00425         TF[4*i+3] = tmp;
00426     }
00427     return true;
00428 }
00429 
00430 
00431 static GLuint createPreintegrationTable( const uint8_t *Table )
00432 {
00433     EQLOG( eq::LOG_CUSTOM ) << "Calculating preintegration table" << std::endl;
00434 
00435     double rInt[256]; rInt[0] = 0.;
00436     double gInt[256]; gInt[0] = 0.;
00437     double bInt[256]; bInt[0] = 0.;
00438     double aInt[256]; aInt[0] = 0.;
00439 
00440     // Creating SAT (Summed Area Tables) from averaged neighbouring RGBA values
00441     for( int i=1; i<256; i++ )
00442     {
00443         // average Alpha from two neighbouring TF values
00444         const double tauc =   ( Table[(i-1)*4+3] + Table[i*4+3] ) / 2. / 255.;
00445 
00446         // SAT of average RGBs from two neighbouring TF values 
00447         // multiplied with Alpha
00448         rInt[i] = rInt[i-1] + ( Table[(i-1)*4+0] + Table[i*4+0] )/2.*tauc;
00449         gInt[i] = gInt[i-1] + ( Table[(i-1)*4+1] + Table[i*4+1] )/2.*tauc;
00450         bInt[i] = bInt[i-1] + ( Table[(i-1)*4+2] + Table[i*4+2] )/2.*tauc;
00451 
00452         // SAT of average Alpha values
00453         aInt[i] = aInt[i-1] + tauc;
00454     }
00455 
00456     std::vector<GLubyte> lookupImg( 256*256*4, 255 ); // Preint Texture
00457 
00458     int lookupindex=0;
00459 
00460     for( int sb=0; sb<256; sb++ )
00461     {
00462         for( int sf=0; sf<256; sf++ )
00463         {
00464             int smin, smax;
00465             if( sb<sf ) { smin=sb; smax=sf; }
00466             else        { smin=sf; smax=sb; }
00467 
00468             double rcol, gcol, bcol, acol;
00469             if( smax != smin )
00470             {
00471                 const double factor = 1. / (double)(smax-smin);
00472                 rcol =       ( rInt[smax] - rInt[smin] ) * factor;
00473                 gcol =       ( gInt[smax] - gInt[smin] ) * factor;
00474                 bcol =       ( bInt[smax] - bInt[smin] ) * factor;
00475                 acol = exp( -( aInt[smax] - aInt[smin] ) * factor) * 255.;
00476             } else
00477             {
00478                 const int    index  = smin*4;
00479                 const double factor = 1./255.;
00480                 rcol =       Table[index+0] * Table[index+3] * factor;
00481                 gcol =       Table[index+1] * Table[index+3] * factor;
00482                 bcol =       Table[index+2] * Table[index+3] * factor;
00483                 acol = exp(                  -Table[index+3] * factor ) * 256.;
00484             }
00485             lookupImg[lookupindex++] = clip( static_cast<int>(rcol), 0, 255 );
00486             lookupImg[lookupindex++] = clip( static_cast<int>(gcol), 0, 255 );
00487             lookupImg[lookupindex++] = clip( static_cast<int>(bcol), 0, 255 );
00488             lookupImg[lookupindex++] = clip( static_cast<int>(acol), 0, 255 );
00489         }
00490     }
00491 
00492     GLuint preintName = 0;
00493     glGenTextures( 1, &preintName );
00494     glBindTexture( GL_TEXTURE_2D, preintName );
00495     glTexImage2D(  GL_TEXTURE_2D, 0, GL_RGBA, 256, 256, 0, 
00496                    GL_RGBA, GL_UNSIGNED_BYTE, &lookupImg[0] );
00497 
00498     glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S,     GL_CLAMP_TO_EDGE );
00499     glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE );
00500     glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR        );
00501     glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR        );
00502 
00503     return preintName;
00504 }
00505 
00506 
00507 }
Generated on Fri Jun 8 2012 15:44:32 for Equalizer 1.2.1 by  doxygen 1.8.0