modules/stdEncoder.cpp

   1 #include <memory> // auto_ptr
   2
   3 #include "stdEncoder.h"
   4 #include "../fileUtil.h"
   5
   6 using namespace std;
   7
   8 /// \file
   9
  10 /** Quantizing stuff used by MStdEncoder \relates MStdEncoder */
  11 namespace Quantizer {
  12         /** Quantizes f that belongs to [0;possib/2^scale] into [0;possib-1] */
  13         inline int quantizeByPower(Real f,int scale,int possib) {
  14                 ASSERT( f>=0 && f<=possib/(Real)powers[scale] );
  15                 int result= (int)trunc(ldexp(f,scale));
  16                 ASSERT( result>=0 && result<=possib );
  17                 return result<possib ? result : --result;
  18         }
  19         /** Performs the opposite to ::quantizeByPower */
  20         inline Real dequantizeByPower(int i,int scale,int DEBUG_ONLY(possib)) {
  21                 ASSERT( i>=0 && i<possib );
  22                 Real result= ldexp(i+Real(0.5),-scale);
  23                 ASSERT( result>=0 && result<= possib/(Real)powers[scale] );
  24                 return result;
  25         }
  26
  27         class QuantBase {
  28         protected:
  29                 int scale       ///  how much should the values be scaled (like in ::quantizeByPower)
  30                 , possib;       ///< the number of possibilities for quantization
  31         public:
  32                 /** Quantizes a value */
  33                 int quant(Real val) const
  34                         { return quantizeByPower(val,scale,possib); }
  35                 /** Dequantizes a value */
  36                 Real dequant(int i) const
  37                         { return dequantizeByPower(i,scale,possib); }
  38                 /** Quant-rounds a value - returns its state after quantization and dequantization */
  39                 Real qRound(Real val) const
  40                         { return dequant(quant(val)); }
  41         };
  42
  43         /** (De)%Quantizer for range-block averages, only initializes QuantBase correctly */
  44         class Average: public QuantBase {
  45         public:
  46                 Average(int possibLog2) {
  47                         ASSERT(possibLog2>0);
  48                 //      the average is from [0,1]
  49                         scale= possibLog2;
  50                         possib= powers[possibLog2];
  51                 }
  52         };
  53
  54         /** (De)%Quantizer for range-block deviations, only initializes QuantBase correctly */
  55         class Deviation: public QuantBase {
  56         public:
  57                 Deviation(int possibLog2) {
  58                         ASSERT(possibLog2>0);
  59                 //      the deviation is from [0,0.5] -> we have to scale twice more (log2+=1)
  60                         scale= possibLog2+1;
  61                         possib= powers[possibLog2];
  62                 }
  63         };
  64 } // Quantizer namespace
  65
  66
  67 /** Adjust block of a domain to be mapped equally on an incomplete range (depends on rotation) */
  68 inline static Block adjustDomainForIncompleteRange( Block range, int rotation, Block domain ) {
  69         // 0, 0T: top left
  70         // 1, 1T: top right
  71         // 2, 2T: bottom right
  72         // 3, 3T: bottom left
  73
  74         ASSERT( 0<=rotation && rotation<8 );
  75         int rotID= rotation/2;
  76
  77         int xSize= range.width();
  78         int ySize= range.height();
  79
  80         int magic= rotID+rotation%2;
  81         if ( magic==1 || magic==3 ) // for rotations 0T, 1, 2T, 3 -> swapped x and y
  82                 swap(xSize,ySize);
  83
  84         if ( rotID>1 )  // domain rotations 2, 3, 2T, 3T -> aligned to the bottom
  85                 domain.y0= domain.yend-ySize;
  86         else                    // other rotations are aligned to the top
  87                 domain.yend= domain.y0+ySize;
  88
  89         if ( rotID==1 || rotID==2 )     // rotations 1, 2, 1T, 2T -> aligned to the right
  90                 domain.x0= domain.xend-xSize;
  91         else                                            // other rotations are aligned to the left
  92                 domain.xend= domain.x0+xSize;
  93
  94         return domain;
  95 }
  96
  97 ////    Member methods
  98
  99 void MStdEncoder::initialize( IRoot::Mode mode, PlaneBlock &planeBlock_ ) {
 100         planeBlock= &planeBlock_;
 101 //      prepare the domains-module
 102         planeBlock->domains->initPools(*planeBlock);
 103
 104         int maxLevel= 1 +log2ceil(max( planeBlock->width, planeBlock->height ));
 105 //      prepare levelPoolInfos
 106         levelPoolInfos.resize(maxLevel);
 107
 108         if (mode==IRoot::Encode) {
 109                 typedef ISquareDomains::PoolList PoolList;
 110         //      prepare the domains
 111                 planeBlock->domains->fillPixelsInPools(*planeBlock);
 112                 for_each( planeBlock->domains->getPools(), mem_fun_ref(&Pool::summers_makeValid) );
 113         //      initialize the range summers
 114                 planeBlock->summers_makeValid();
 115
 116         //      prepare maximum SquareErrors allowed for regular range blocks
 117                 stdRangeSEs.resize(maxLevel+1);
 118                 planeBlock->settings->moduleQ2SE->regularRangeErrors
 119                         ( planeBlock->settings->quality, maxLevel+1, &stdRangeSEs.front() );
 120         }
 121 }
 122
 123
 124 namespace NOSPACE {
 125         struct LevelPoolInfo_indexComparator {
 126                 typedef ISquareEncoder::LevelPoolInfo LevelPoolInfo;
 127                 bool operator()( const LevelPoolInfo &a, const LevelPoolInfo &b )
 128                         { return a.indexBegin < b.indexBegin; }
 129         };
 130 }
 131 MStdEncoder::PoolInfos::const_iterator MStdEncoder
 132 ::getPoolFromDomID( int domID, const PoolInfos &poolInfos ) {
 133         ASSERT( domID>=0 && domID<poolInfos.back().indexBegin );
 134 //      get the right pool
 135         PoolInfos::value_type toFind;
 136         toFind.indexBegin= domID;
 137         return upper_bound( poolInfos.begin(), poolInfos.end(), toFind
 138                 , LevelPoolInfo_indexComparator() ) -1;
 139 }
 140
 141 const ISquareDomains::Pool& MStdEncoder::getDomainData
 142 ( const RangeNode &rangeBlock, const ISquareDomains::PoolList &pools
 143 , const PoolInfos &poolInfos, int domIndex, int zoom, Block &block ) {
 144
 145 //      get the pool
 146         PoolInfos::const_iterator it= getPoolFromDomID( domIndex, poolInfos );
 147         const Pool &pool= pools[ it - poolInfos.begin() ];
 148 //      get the coordinates
 149         int indexInPool= domIndex - it->indexBegin;
 150         ASSERT( indexInPool>=0 && indexInPool<(it+1)->indexBegin );
 151         int sizeNZ= powers[rangeBlock.level-zoom];
 152         int zoomFactor= powers[zoom];
 153 //      changed: the domains are mapped along columns and not rows
 154         int domsInCol= getCountForDensity( pool.height/zoomFactor, it->density, sizeNZ );
 155         block.x0= (indexInPool/domsInCol)*it->density*zoomFactor;
 156         block.y0= (indexInPool%domsInCol)*it->density*zoomFactor;
 157         block.xend= block.x0+powers[rangeBlock.level];
 158         block.yend= block.y0+powers[rangeBlock.level];
 159
 160         return pool;
 161 }
 162
 163
 164 namespace NOSPACE {
 165         typedef IStdEncPredictor::NewPredictorData StableInfo;
 166
 167         /** Information about the best (found so far) domain block for a range block */
 168         struct BestInfo: public IStdEncPredictor::Prediction {
 169                 float error;    ///< The square error of the mapping
 170                 bool inverted;  ///< Whether the colours are mapped with a negative linear coefficient
 171                 #ifndef NDEBUG
 172                 Real rdSum              ///
 173                 , dSum                  ///  Sum of all pixels in best domain
 174                 , d2Sum;                ///< Sum of squares of all pixesl in the best domain
 175                 #endif
 176
 177                 /** Only initializes ::error to the maximum float value */
 178                 BestInfo()
 179                 : Prediction(), error( numeric_limits<float>::infinity() ) {}
 180                 /** Only returns reference to *this typed as IStdEncPredictor::Prediciton */
 181                 IStdEncPredictor::Prediction& prediction()
 182                         { return *this; }
 183         };
 184 } // namespace NOSPACE
 185
 186 /** Structure constructed for a range to try domains */
 187 struct MStdEncoder::EncodingInfo {
 188         /** The type for exact-comparing methods */
 189         typedef bool (EncodingInfo::*ExactCompareProc)( Prediction prediction );
 190
 191 private:
 192         static const ExactCompareProc exactCompareArray[];      ///< Possible comparing methods
 193
 194         ExactCompareProc selectedProc;                                          ///< Selected comparing method
 195 public:
 196         StableInfo stable;      ///< Information only depending on the range (and not domain) block
 197         BestInfo best;          ///< Information about the best (at the moment) mapping found
 198         float targetSE;         ///< The target SE (square error) for the range
 199
 200 public:
 201         /** Only nulls ::selectedProc */
 202         EncodingInfo()
 203         : selectedProc(0) {}
 204
 205         /** Initializes a RangeInfo object from \p this */
 206         RangeInfo* initRangeInfo(RangeInfo *ri) const {
 207                 ri->bestSE=             best.error;
 208                 ri->domainID=   best.domainID;
 209                 ri->rotation=   best.rotation;
 210                 ri->qrAvg=              stable.qrAvg;
 211                 ri->qrDev2=             stable.qrDev2;
 212                 ri->inverted=   best.inverted;
 213                 #ifndef NDEBUG
 214                 ri->exact.avg=  stable.rSum/stable.pixCount;
 215                 ri->exact.dev2= stable.r2Sum/stable.pixCount - sqr(ri->exact.avg);
 216                 ri->exact.linCoeff= ( stable.pixCount*best.rdSum - stable.rSum*best.dSum )
 217                 / ( stable.pixCount*best.d2Sum - sqr(best.dSum) );
 218                 ri->exact.constCoeff= ri->exact.avg - ri->exact.linCoeff*best.dSum/stable.pixCount;
 219                 #endif
 220                 return ri;
 221         }
 222
 223         /** Selects the right comparing method according to ::stable */
 224         void selectExactCompareProc() {
 225                 selectedProc= exactCompareArray
 226                 [(((  stable.quantError *2
 227                         + stable.allowInversion ) *2
 228                         + stable.isRegular ) *2
 229                         + (stable.maxLinCoeff2 >= 0) ) *2
 230                         + (stable.bigScaleCoeff != 0)
 231                 ];
 232         }
 233
 234         /** Uses the selected comparing method */
 235         bool exactCompare(Prediction prediction)
 236                 { return (this->* selectedProc)(prediction); }
 237
 238 private:
 239         /** Template for all the comparing methods */
 240         template < bool quantErrors, bool allowInversion, bool isRegular
 241         , bool restrictMaxLinCoeff, bool bigScalePenalty >
 242         bool exactCompareProc( Prediction prediction );
 243 }; // EncodingInfo struct
 244
 245 #define ALTERNATE_0(name,params...) name<params>
 246 #define ALTERNATE_1(name,params...) ALTERNATE_0(name,##params,false), ALTERNATE_0(name,##params,true)
 247 #define ALTERNATE_2(name,params...) ALTERNATE_1(name,##params,false), ALTERNATE_1(name,##params,true)
 248 #define ALTERNATE_3(name,params...) ALTERNATE_2(name,##params,false), ALTERNATE_2(name,##params,true)
 249 #define ALTERNATE_4(name,params...) ALTERNATE_3(name,##params,false), ALTERNATE_3(name,##params,true)
 250 #define ALTERNATE_5(name,params...) ALTERNATE_4(name,##params,false), ALTERNATE_4(name,##params,true)
 251 const MStdEncoder::EncodingInfo::ExactCompareProc MStdEncoder::EncodingInfo::exactCompareArray[]
 252         = {ALTERNATE_5(&MStdEncoder::EncodingInfo::exactCompareProc)};
 253 #undef ALTERNATE_0
 254 #undef ALTERNATE_1
 255 #undef ALTERNATE_2
 256 #undef ALTERNATE_3
 257 #undef ALTERNATE_4
 258 #undef ALTERNATE_5
 259
 260 template< bool quantErrors, bool allowInversion, bool isRegular
 261 , bool restrictMaxLinCoeff, bool bigScalePenalty >
 262 bool MStdEncoder::EncodingInfo::exactCompareProc( Prediction prediction ) {
 263         using namespace MatrixWalkers;
 264 //      find out which domain was predicted (pixel matrix and position within it)
 265         Block domBlock;
 266         const Pool &pool= getDomainData( *stable.rangeBlock, *stable.pools, *stable.poolInfos
 267         , prediction.domainID, 0/*zoom*/, domBlock );
 268 //      compute domain sums
 269         if (!isRegular)
 270                 domBlock= adjustDomainForIncompleteRange
 271                         ( *stable.rangeBlock, prediction.rotation, domBlock );
 272         Real dSum, d2Sum;
 273         pool.getSums(domBlock).unpack(dSum,d2Sum);
 274 //      compute the denominator common to most formulas
 275         Real test= stable.pixCount*d2Sum - sqr(dSum);
 276         if (test<=0) // skip too flat domains
 277                 return false;
 278         Real denom= 1/test;
 279
 280 //      compute the sum of products of pixels
 281         Real rdSum= walkOperateCheckRotate
 282         ( Checked<const SReal>(stable.rangePixels, *stable.rangeBlock), RDSummer<Real,SReal>()
 283         , pool.pixels, domBlock, prediction.rotation ) .result();
 284
 285         Real nRDs_RsDs= stable.pixCount*rdSum - stable.rSum*dSum;
 286 //      check for negative linear coefficients (if needed)
 287         if (!allowInversion)
 288                 if (nRDs_RsDs<0)
 289                         return false;
 290 //      compute the square of linear coefficient if needed (for restricting or penalty)
 291         Real linCoeff2 DEBUG_ONLY(= numeric_limits<Real>::quiet_NaN() );
 292         if ( restrictMaxLinCoeff || bigScalePenalty )
 293                 linCoeff2= stable.rnDev2 * denom;
 294         if (restrictMaxLinCoeff)
 295                 if ( linCoeff2 > stable.maxLinCoeff2 )
 296                         return false;
 297
 298         float optSE DEBUG_ONLY(= numeric_limits<Real>::quiet_NaN() );
 299
 300         if (quantErrors) {
 301                 optSE= stable.qrAvg * ( stable.pixCount*stable.qrAvg - ldexp(stable.rSum,1) )
 302                         + stable.r2Sum + stable.qrDev
 303                                 * ( stable.pixCount*stable.qrDev - ldexp( abs(nRDs_RsDs)*sqrt(denom), 1 ) );
 304         } else { // !quantErrors
 305         //      assuming different linear coeffitient
 306                 Real inner= stable.rnDev2 - stable.rnDev*abs(nRDs_RsDs)*sqrt(denom);
 307                 optSE= ldexp( inner, 1 ) / stable.pixCount;
 308         }
 309
 310 //      add big-scaling penalty if needed
 311         if (bigScalePenalty)
 312                 optSE+= linCoeff2 * targetSE * pool.contrFactor * stable.bigScaleCoeff;
 313
 314 //      test if the error is the best so far
 315         if ( optSE < best.error ) {
 316                 best.prediction()= prediction;
 317                 best.error= optSE;
 318                 best.inverted= nRDs_RsDs<0;
 319
 320                 #ifndef NDEBUG
 321                 best.rdSum= rdSum;
 322                 best.dSum= dSum;
 323                 best.d2Sum= d2Sum;
 324                 #endif
 325
 326                 return true;
 327         } else
 328                 return false;
 329 } // EncodingInfo::exactCompareProc method
 330
 331
 332 float MStdEncoder::findBestSE(const RangeNode &range,bool allowHigherSE) {
 333         ASSERT( planeBlock && !stdRangeSEs.empty() && !range.encoderData );
 334         const IColorTransformer::PlaneSettings *plSet= planeBlock->settings;
 335
 336 //      initialize an encoding-info object
 337         EncodingInfo info;
 338
 339         info.stable.rangeBlock=         &range;
 340         info.stable.rangePixels=        planeBlock->pixels;
 341         info.stable.pools=                      &planeBlock->domains->getPools();
 342
 343         ASSERT( range.level < (int)levelPoolInfos.size() );
 344         info.stable.poolInfos=          &levelPoolInfos[range.level];
 345 //      check the level has been initialized
 346         if ( info.stable.poolInfos->empty() ) {
 347                 buildPoolInfos4aLevel(range.level);
 348                 ASSERT( !info.stable.poolInfos->empty() );
 349         }
 350
 351         info.stable.allowRotations=     settingsInt(AllowedRotations);
 352         info.stable.quantError=         settingsInt(AllowedQuantError);
 353         info.stable.allowInversion=     settingsInt(AllowedInversion);
 354         info.stable.isRegular=          range.isRegular();
 355         {
 356                 Real coeff= settingsFloat(MaxLinCoeff);
 357                 info.stable.maxLinCoeff2=       coeff==MaxLinCoeff_none ? -1 : sqr(coeff);
 358         }
 359         info.stable.bigScaleCoeff=      settingsFloat(BigScaleCoeff);
 360
 361         planeBlock->getSums(range).unpack( info.stable.rSum, info.stable.r2Sum );
 362         info.stable.pixCount=   range.size();
 363         info.stable.rnDev2=             info.stable.pixCount*info.stable.r2Sum - sqr(info.stable.rSum);
 364         if (info.stable.rnDev2<0)
 365                 info.stable.rnDev2= 0;
 366         info.stable.rnDev=              sqrt(info.stable.rnDev2);
 367
 368         Real variance;
 369         {
 370                 Quantizer::Average quantAvg( settingsInt(QuantStepLog_avg) );
 371                 Quantizer::Deviation quantDev( settingsInt(QuantStepLog_dev) );
 372                 Real average= info.stable.rSum / info.stable.pixCount;
 373                 info.stable.qrAvg= quantAvg.qRound(average);
 374
 375                 variance= info.stable.r2Sum/info.stable.pixCount - sqr(average);
 376                 Real deviance= variance>0 ? sqrt(variance) : 0;
 377                 int qrDev= quantDev.quant(deviance);
 378         //      if we have too little deviance or no domain pool for that big level or no domain in the pool
 379                 if ( !qrDev || range.level >= (int)levelPoolInfos.size()
 380                 || info.stable.poolInfos->back().indexBegin <= 0 ) // -> no domain block, only average
 381                         goto returning; // skips to the end, assigning a constant block
 382                 else { // the regular case, with nonzero quantized deviance
 383                         info.stable.qrDev= quantDev.dequant(qrDev);
 384                         info.stable.qrDev2= sqr(info.stable.qrDev);
 385                 }
 386         }
 387
 388         info.targetSE= info.best.error= info.stable.isRegular ? stdRangeSEs[range.level]
 389                 : plSet->moduleQ2SE->rangeSE( plSet->quality, range.size() );
 390         ASSERT(info.targetSE>=0);
 391         if (allowHigherSE)
 392                 info.best.error= numeric_limits<float>::infinity();
 393         info.selectExactCompareProc();
 394
 395         { // a goto-skippable block
 396         //      create and initialize a new predictor (in auto_ptr because of exceptions)
 397                 auto_ptr<IStdEncPredictor::IOneRangePredictor> predictor
 398                 ( modulePredictor()->newPredictor(info.stable) );
 399                 typedef IStdEncPredictor::Predictions Predictions;
 400                 Predictions predicts;
 401
 402                 float sufficientSE= info.targetSE*settingsFloat(SufficientSEq);
 403         //      get and process prediction chunks until an empty one is returned
 404                 while ( !predictor->getChunk(info.best.error,predicts).empty() )
 405                         for (Predictions::iterator it=predicts.begin(); it!=predicts.end(); ++it) {
 406                                 bool betterSE= info.exactCompare(*it);
 407                                 if ( betterSE && info.best.error<=sufficientSE )
 408                                         goto returning;
 409                         }
 410         }
 411
 412         returning:
 413 //      check the case that the predictor didn't return any domains (or jumped to returning:)
 414         if ( info.best.domainID < 0 ) { // a constant block
 415                 info.stable.qrDev= info.stable.qrDev2= 0;
 416                 info.best.error= info.stable.quantError
 417                         ? info.stable.r2Sum + info.stable.qrAvg
 418                                 *( info.stable.pixCount*info.stable.qrAvg - ldexp(info.stable.rSum,1) )
 419                         : variance*info.stable.pixCount;
 420         }
 421 //      store the important info and return the error
 422         range.encoderData= info.initRangeInfo( /*rangeInfoAlloc.make()*/ new RangeInfo );
 423         return info.best.error;
 424 } // ::findBestSE method
 425
 426 void MStdEncoder::buildPoolInfos4aLevel(int level) {
 427 //      get the real maximum domain count (divide by the number of rotations)
 428         int domainCountLog2= planeBlock->settings->domainCountLog2;
 429         if ( settingsInt(AllowedRotations) )
 430                 domainCountLog2-= 3;
 431 //      get the per-domain-pool densities
 432         vector<short> densities= planeBlock->domains->getLevelDensities(level,domainCountLog2);
 433 //      storing the result in this-level infos with beginIndex values for all pools
 434         vector<LevelPoolInfo> &poolInfos= levelPoolInfos[level];
 435         ASSERT( poolInfos.empty() );
 436         const ISquareDomains::PoolList &pools= planeBlock->domains->getPools();
 437         int zoom= planeBlock->settings->zoom;
 438
 439         int domainSizeNZ= powers[level-zoom];
 440         int poolCount= pools.size();
 441         poolInfos.resize(poolCount+1);
 442
 443         int domCount= poolInfos[0].indexBegin= 0; // accumulated count
 444         for (int i=0; i<poolCount; ++i) {
 445                 int dens= poolInfos[i].density= densities[i];
 446                 ASSERT(dens>=0);
 447                 if (dens) // if dens==0, there are no domains -> no increase
 448                         domCount+= getCountForDensity2D( rShift(pools[i].width,zoom)
 449                                 , rShift(pools[i].height,zoom), dens, domainSizeNZ );
 450                 poolInfos[i+1].indexBegin= domCount;
 451         }
 452         poolInfos[poolCount].density= -1;
 453 }
 454
 455 void MStdEncoder::writeSettings(ostream &file) {
 456         ASSERT( /*modulePredictor() &&*/ moduleCodec(true) && moduleCodec(false) );
 457 //      put settings needed for decoding
 458         put<Uchar>( file, settingsInt(AllowedRotations) );
 459         put<Uchar>( file, settingsInt(AllowedInversion) );
 460         put<Uchar>( file, settingsInt(QuantStepLog_avg) );
 461         put<Uchar>( file, settingsInt(QuantStepLog_dev) );
 462 //      put ID's of connected modules (the predictor module doesn't need to be known)
 463         file_saveModuleType( file, ModuleCodecAvg );
 464         file_saveModuleType( file, ModuleCodecDev );
 465 //      put the settings of connected modules
 466         moduleCodec(true)->writeSettings(file);
 467         moduleCodec(false)->writeSettings(file);
 468 }
 469 void MStdEncoder::readSettings(istream &file) {
 470         ASSERT( !modulePredictor() && !moduleCodec(true) && !moduleCodec(false) );
 471 //      get settings needed for decoding
 472         settingsInt(AllowedRotations)= get<Uchar>(file);
 473         settingsInt(AllowedInversion)= get<Uchar>(file);
 474         settingsInt(QuantStepLog_avg)= get<Uchar>(file);
 475         settingsInt(QuantStepLog_dev)= get<Uchar>(file);
 476 //      create connected modules (the predictor module isn't needed)
 477         file_loadModuleType( file, ModuleCodecAvg );
 478         file_loadModuleType( file, ModuleCodecDev );
 479 //      get the settings of connected modules
 480         moduleCodec(true)->readSettings(file);
 481         moduleCodec(false)->readSettings(file);
 482 }
 483
 484 void MStdEncoder::writeData(ostream &file,int phase) {
 485         typedef RangeList::const_iterator RLcIterator;
 486         ASSERT( moduleCodec(true) && moduleCodec(false) );
 487         ASSERT( planeBlock && planeBlock->ranges && planeBlock->encoder==this );
 488         const RangeList &ranges= planeBlock->ranges->getRangeList();
 489         ASSERT( !ranges.empty() );
 490
 491         switch(phase) {
 492                 case 0: { // save the averages
 493                 //      get the list
 494                         vector<int> averages;
 495                         averages.reserve(ranges.size());
 496                         Quantizer::Average quant( settingsInt(QuantStepLog_avg) );
 497                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 498                                 ASSERT( *it && (*it)->encoderData );
 499                                 STREAM_POS(file);
 500                                 const RangeInfo *info= RangeInfo::get(*it);
 501                                 averages.push_back( quant.quant(info->qrAvg) );
 502                         }
 503                 //      pass it to the codec
 504                         moduleCodec(true)->setPossibilities( powers[settingsInt(QuantStepLog_avg)] );
 505                         moduleCodec(true)->encode( averages, file );
 506                         break;
 507                 }
 508
 509                 case 1: { // save the deviations
 510                 //      get the list
 511                         vector<int> devs;
 512                         devs.reserve(ranges.size());
 513                         Quantizer::Deviation quant( settingsInt(QuantStepLog_dev) );
 514                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 515                                 ASSERT( *it && (*it)->encoderData );
 516                                 STREAM_POS(file);
 517                                 const RangeInfo *info= RangeInfo::get(*it);
 518                                 int dev= ( info->domainID>=0 ? quant.quant(sqrt(info->qrDev2)) : 0 );
 519                                 devs.push_back(dev);
 520                         }
 521                 //      pass it to the codec
 522                         moduleCodec(false)->setPossibilities( powers[settingsInt(QuantStepLog_dev)] );
 523                         moduleCodec(false)->encode( devs, file );
 524                         break;
 525                 }
 526
 527                 case 2: { // save the rest
 528                         BitWriter stream(file);
 529                 //      put the inversion bits if inversion is allowed
 530                         if ( settingsInt(AllowedInversion) )
 531                                 for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 532                                         ASSERT( *it && (*it)->encoderData );
 533                                         STREAM_POS(file);
 534                                         const RangeInfo *info= RangeInfo::get(*it);
 535                                         if ( info->domainID >= 0 )
 536                                                 stream.putBits(info->inverted,1);
 537                                 }
 538                 //      put the rotation bits if rotations are allowed
 539                         if ( settingsInt(AllowedRotations) )
 540                                 for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 541                                         ASSERT( *it && (*it)->encoderData );
 542                                         STREAM_POS(file);
 543                                         const RangeInfo *info= RangeInfo::get(*it);
 544                                         if ( info->domainID >= 0 ) {
 545                                                 ASSERT( 0<=info->rotation && info->rotation<8 );
 546                                                 stream.putBits( info->rotation, 3 );
 547                                         }
 548                                 }
 549                 //      find out bits needed to store IDs of domains for every level
 550                         vector<int> domBits;
 551                         domBits.resize( levelPoolInfos.size() );
 552                         for (int i=0; i<(int)levelPoolInfos.size(); ++i) {
 553                                 int count= levelPoolInfos[i].empty() ? 0 : levelPoolInfos[i].back().indexBegin;
 554                                 ASSERT(count>=0);
 555                                 STREAM_POS(file);
 556                                 domBits[i]= count ? log2ceil(count) : 0;
 557                         }
 558                 //      put the domain bits
 559                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 560                                 ASSERT( *it && (*it)->encoderData );
 561                                 STREAM_POS(file);
 562                                 const RangeInfo *info= RangeInfo::get(*it);
 563                                 int domID= info->domainID;
 564                                 if ( domID >= 0 ) {
 565                                         int bits= domBits[ (*it)->level ];
 566                                         ASSERT( 0<=domID && domID<powers[bits] );
 567                                         if (bits>0)
 568                                                 stream.putBits( domID, bits );
 569                                 }
 570                         }
 571                         break;
 572                 }
 573                 default:
 574                         ASSERT(false);
 575         }
 576 } // ::writeData method
 577 void MStdEncoder::readData(istream &file,int phase) {
 578         typedef RangeList::const_iterator RLcIterator;
 579         ASSERT( moduleCodec(true) && moduleCodec(false) );
 580         ASSERT( planeBlock && planeBlock->ranges && planeBlock->encoder==this );
 581         const RangeList &ranges= planeBlock->ranges->getRangeList();
 582         ASSERT( !ranges.empty() );
 583
 584         switch(phase) {
 585                 case 0: { // create the infos and load the averages
 586                 //      get the list of averages from the codec
 587                         vector<int> averages;
 588                         moduleCodec(true)->setPossibilities( powers[settingsInt(QuantStepLog_avg)] );
 589                         moduleCodec(true)->decode( file, ranges.size(), averages );
 590                 //      iterate the list, create infos and fill them with dequantized averages
 591                         Quantizer::Average quant( settingsInt(QuantStepLog_avg) );
 592                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 593                                 ASSERT( *it && !(*it)->encoderData );
 594                                 STREAM_POS(file);
 595                                 RangeInfo *info= new RangeInfo;//rangeInfoAlloc.make();
 596                                 (*it)->encoderData= info;
 597                                 info->qrAvg= quant.dequant(averages[it-ranges.begin()]);
 598                                 DEBUG_ONLY( info->bestSE= -1; ) // SE not needed for decoding,saving,...
 599                                 if ( !settingsInt(AllowedInversion) )
 600                                         info->inverted= false;
 601                         }
 602                         break;
 603                 }
 604
 605                 case 1: { // load the deviations
 606                 //      get the list of deviations from the codec
 607                         vector<int> devs;
 608                         moduleCodec(false)->setPossibilities( powers[settingsInt(QuantStepLog_dev)] );
 609                         moduleCodec(false)->decode( file, ranges.size(), devs );
 610                 //      iterate the list, create infos and fill the with dequantized deviations
 611                         Quantizer::Deviation quant( settingsInt(QuantStepLog_dev) );
 612                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 613                                 ASSERT( *it && (*it)->encoderData );
 614                                 STREAM_POS(file);
 615                                 RangeInfo *info= RangeInfo::get(*it);
 616                                 int quantDev= devs[it-ranges.begin()];
 617                                 if (quantDev)
 618                                         info->qrDev2= sqr(quant.dequant(quantDev));
 619                                 else { // it's a flat block
 620                                         info->qrDev2= 0;
 621
 622                                         DEBUG_ONLY( info->inverted= false; )
 623                                         ASSERT( info->rotation==-1 && info->domainID==-1 );
 624                                 }
 625                         }
 626                         break;
 627                 }
 628
 629                 case 2: { // load the rest
 630                         BitReader stream(file);
 631                 //      get the inversion bits if inversion is allowed
 632                         if ( settingsInt(AllowedInversion) )
 633                                 for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 634                                         ASSERT( *it && (*it)->encoderData );
 635                                         STREAM_POS(file);
 636                                         RangeInfo *info= RangeInfo::get(*it);
 637                                         if (info->qrDev2)
 638                                                 info->inverted= stream.getBits(1);
 639                                 }
 640                 //      get the rotation bits if rotations are allowed
 641                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 642                                 ASSERT( *it && (*it)->encoderData );
 643                                 STREAM_POS(file);
 644                                 RangeInfo *info= RangeInfo::get(*it);
 645                                 if (info->qrDev2)
 646                                         info->rotation= settingsInt(AllowedRotations) ? stream.getBits(3) : 0;
 647                         }
 648                 //      find out bits needed to store IDs of domains for every level
 649                         vector<int> domBits( levelPoolInfos.size(), -1 );
 650                 //      get the domain bits
 651                         for (RLcIterator it=ranges.begin(); it!=ranges.end(); ++it) {
 652                                 ASSERT( *it && (*it)->encoderData );
 653                                 STREAM_POS(file);
 654                                 RangeInfo *info= RangeInfo::get(*it);
 655                                 if (!info->qrDev2)
 656                                         continue;
 657                                 int level= (*it)->level;
 658                                 int bits= domBits[level];
 659                                 if (bits<0) {
 660                                 //      yet unused level -> initialize it and get the right bits
 661                                         buildPoolInfos4aLevel(level);
 662                                         bits= domBits[level]
 663                                         = log2ceil( levelPoolInfos[level].back().indexBegin );
 664                                         ASSERT( bits>0 );
 665                                 }
 666                         //      get the domain ID, check it's OK and store it in the range's info
 667                                 int domID= stream.getBits(bits);
 668                                 checkThrow( domID < levelPoolInfos[level].back().indexBegin );
 669                                 info->domainID= domID;
 670                         }
 671                 //      initialize decoding accelerators for domain-to-range mappings
 672                         initRangeInfoAccelerators();
 673                         break;
 674                 }
 675                 default:
 676                         ASSERT(false);
 677         }
 678 } // ::readData method
 679
 680 void MStdEncoder::decodeAct( DecodeAct action, int count ) {
 681 //      do some checks
 682         ASSERT( planeBlock && planeBlock->ranges && planeBlock->encoder==this );
 683         const RangeList &ranges= planeBlock->ranges->getRangeList();
 684         ASSERT( !ranges.empty() );
 685
 686         switch (action) {
 687         default:
 688                 ASSERT(false);
 689         case Clear:
 690                 ASSERT(count==1);
 691                 planeBlock->pixels.fillSubMatrix
 692                         ( Block(0,0,planeBlock->width,planeBlock->height), 0.5f );
 693                 planeBlock->summers_invalidate();
 694                 break;
 695         case Iterate:
 696                 ASSERT(count>0);
 697                 do {
 698                 //      prepare the domains, iterate each range block
 699                         planeBlock->domains->fillPixelsInPools(*planeBlock);
 700                         for (RangeList::const_iterator it=ranges.begin(); it!=ranges.end(); ++it) {
 701                                 const RangeInfo &info= *RangeInfo::get(*it);
 702                                 if ( info.domainID < 0 ) { // no domain - constant color
 703                                         planeBlock->pixels.fillSubMatrix( **it, info.qrAvg );
 704                                         continue;
 705                                 }
 706                         //      get domain sums and the pixel count
 707                                 Real dSum, d2Sum;
 708                                 info.decAccel.pool->summers_makeValid();
 709                                 info.decAccel.pool->getSums(info.decAccel.domBlock).unpack(dSum,d2Sum);
 710                                 Real pixCount= (*it)->size();
 711                         //      find out the coefficients and handle constant blocks
 712                                 Real linCoeff= (info.inverted ? -pixCount : pixCount)
 713                                         * sqrt( info.qrDev2 / ( pixCount*d2Sum - sqr(dSum) ) );
 714                                 if ( !isnormal(linCoeff) || !linCoeff ) {
 715                                         planeBlock->pixels.fillSubMatrix( **it, info.qrAvg );
 716                                         continue;
 717                                 }
 718                                 Real constCoeff= info.qrAvg - linCoeff*dSum/pixCount;
 719                         //      map the nonconstant blocks
 720                                 using namespace MatrixWalkers;
 721                                 MulAddCopyChecked<Real> oper( linCoeff, constCoeff, 0, 1 );
 722                                 walkOperateCheckRotate( Checked<SReal>(planeBlock->pixels, **it), oper
 723                                 , info.decAccel.pool->pixels, info.decAccel.domBlock, info.rotation );
 724                         }
 725                 } while (--count);
 726                 break;
 727         } // switch (action)
 728 } // ::decodeAct method
 729
 730 void MStdEncoder::initRangeInfoAccelerators() {
 731 //      get references that are the same for all range blocks
 732         const RangeList &ranges= planeBlock->ranges->getRangeList();
 733         const ISquareDomains::PoolList &pools= planeBlock->domains->getPools();
 734         int zoom= planeBlock->settings->zoom;
 735 //      iterate over the ranges (and their accelerators)
 736         for ( RangeList::const_iterator it=ranges.begin(); it!=ranges.end(); ++it ) {
 737         //      get range level, reference to the range's info and to the level's pool infos
 738                 int level= (*it)->level;
 739                 RangeInfo &info= *RangeInfo::get(*it);
 740                 if ( info.domainID < 0 )
 741                         continue;
 742                 const PoolInfos &poolInfos= levelPoolInfos[level];
 743         //      build pool infos for the level (if neccesary), fill info's accelerators
 744                 if ( poolInfos.empty() )
 745                         buildPoolInfos4aLevel(level), ASSERT( !poolInfos.empty() );
 746                 info.decAccel.pool= &getDomainData
 747                         ( **it, pools, poolInfos, info.domainID, zoom, info.decAccel.domBlock );
 748         //      adjust domain's block if the range isn't regular
 749                 if ( !(*it)->isRegular() )
 750                         info.decAccel.domBlock= adjustDomainForIncompleteRange
 751                                 ( **it, info.rotation, info.decAccel.domBlock );
 752         }
 753 }