modules/stdEncoder.cpp

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