modules/stdDomains.cpp

   1 #include "stdDomains.h"
   2 #include "../fileUtil.h"
   3
   4 enum { MinDomSize=8, MinRngSize=4 };
   5
   6 using namespace std;
   7
   8
   9 /** Implementations of shrinking routines by using MatrixWalkers */
  10 namespace NOSPACE {
  11         using namespace MatrixWalkers;
  12         /** Performs a simple 50\%^2 image shrink (dimensions belong to the destination)
  13          *      \relates MStdDomains */
  14         void shrinkToHalf( CSMatrix src, SMatrix dest, int width, int height ) {
  15                 walkOperate( Checked<SReal>(dest,Block(0,0,width,height))
  16                 , HalfShrinker<const SReal>(src), ReverseAssigner() );
  17         }
  18         /** Performs a simple 33\%x66\% image horizontal shrink
  19          *      (dimensions belong to the destination) \relates MStdDomains */
  20         void shrinkHorizontally( CSMatrix src, SMatrix dest, int width, int height ) {
  21                 walkOperate( Checked<SReal>(dest,Block(0,0,width,height))
  22                 , HorizShrinker<const SReal>(src), ReverseAssigner() );
  23         }
  24         /** Performs a simple 66\%x33\% image vertical shrink
  25          *      (dimensions belong to the destination) \relates MStdDomains */
  26         void shrinkVertically( CSMatrix src, SMatrix dest, int width, int height ) {
  27                 walkOperate( Checked<SReal>(dest,Block(0,0,width,height))
  28                 , VertShrinker<const SReal>(src), ReverseAssigner() );
  29         }
  30         /** Performs 50\% shrink with 45-degree anticlockwise rotation
  31          *      (\p side - the length of the destination square;
  32          *      \p sx0, \p sy0 - the top-left of the enclosing source square) \relates MStdDomains */
  33         void shrinkToDiamond( CSMatrix src, SMatrix dest, int side ) {
  34                 walkOperate( Checked<SReal>(dest,Block(0,0,side,side))
  35                 , DiamShrinker<const SReal>(src,side), ReverseAssigner() );
  36         }//     shrinkToDiamond
  37 }
  38
  39
  40 namespace NOSPACE {
  41         /** Pool ordering according to Pool::type (primary key) and Pool::level (secondary) */
  42         struct PoolTypeLevelComparator {
  43                 typedef MStdDomains::Pool Pool;
  44                 bool operator()(const Pool &a,const Pool &b) {
  45                         if (a.type!=b.type)
  46                                 return a.type<b.type;
  47                         else
  48                                 return a.level<b.level;
  49                 }
  50         };
  51         /* All sizes are unzoomed */
  52         inline static int minSizeNeededForDiamond()
  53                 { return 2*MinDomSize-1; }
  54         inline static int getDiamondSize(int fromSize)
  55                 { return fromSize<minSizeNeededForDiamond() ? 0 : fromSize/2; }
  56         inline static int getDiamondShift(int fromSize)
  57                 { return ( getDiamondSize(fromSize) - MinRngSize + 1 )*2; }
  58 }
  59 void MStdDomains::initPools(const PlaneBlock &planeBlock) {
  60         zoom=   planeBlock.settings->zoom;
  61         width=  rShift( planeBlock.width, zoom );
  62         height= rShift( planeBlock.height, zoom );
  63 //      checks some things
  64         ASSERT( width>0 && height>0 && this && settings && pools.empty() );
  65         if ( min(width,height)/2 < MinDomSize )
  66         //      no domains (too small)
  67                 return;
  68 //      create the first set of domain pools for standard, horizontal and vertical domains
  69         if ( settingsInt(DomPortion_Standard) )
  70                 pools.push_back(Pool( width/2, height/2, DomPortion_Standard, 1, 0.25, zoom ));
  71         if ( settingsInt(DomPortion_Horiz) && width/3>=MinDomSize )
  72                 pools.push_back(Pool( width/3, height/3*2, DomPortion_Horiz, 1, 2.0/9.0, zoom ));
  73         if ( settingsInt(DomPortion_Vert) && height/3>=MinDomSize )
  74                 pools.push_back(Pool( width/3*2, height/3, DomPortion_Vert, 1, 2.0/9.0, zoom ));
  75 //      create the first set of domain pools for diamond domains
  76         if ( settingsInt(DomPortion_Diamond) ) {
  77         //      get longer and shorter dimension
  78                 int shorter= min(width,height);
  79                 int shift= getDiamondShift(shorter);
  80         //      generate the pool sizes
  81                 for (int longer=max(width,height); longer>minSizeNeededForDiamond(); longer-=shift) {
  82                         int side= getDiamondSize( min(longer,shorter) );
  83                         ASSERT(side>=0);
  84                         if (side)
  85                                 pools.push_back(Pool( side, side, DomPortion_Diamond, 1, 0.5, zoom ));
  86                         else
  87                                 break;
  88                 }
  89         }
  90 //      if allowed, create more downscaled domains (at most 256 pools)
  91         if ( settingsInt(MultiDownScaling) )
  92                 for (Uint i=0; i<pools.size() && i<256; ++i) {
  93                         const Pool &pool= pools[i];
  94                 //      compute new dimensions and add the pool if it's big enough
  95                         int w= rShift<int>(pool.width,zoom+1);
  96                         int h= rShift<int>(pool.height,zoom+1);
  97                         float cf= ldexp(pool.contrFactor,-2);
  98                         if ( min(w,h) >= MinDomSize )
  99                                 pools.push_back( Pool( w, h, pool.type, pool.level+1, cf, zoom ) );
 100                 }
 101 //      sort the pools according to their types and levels (stable so the diamonds can't be swapped)
 102         stable_sort( pools.begin(), pools.end(), PoolTypeLevelComparator() );
 103 }
 104
 105 namespace NOSPACE {
 106         typedef MStdDomains::PoolList::const_iterator PoolIt;
 107         /** To be called before creating shrinked domains to check the shrink is OK (debug only) */
 108         static inline bool halfShrinkOK(PoolIt src,PoolIt dest,int zoom) {
 109                 return src->level+1 == dest->level && src->type == dest->type
 110                         && rShift<int>(src->width,zoom+1) == rShift<int>(dest->width,zoom)
 111                         && rShift<int>(src->height,zoom+1) == rShift<int>(dest->height,zoom);
 112         }
 113 }
 114 void MStdDomains::fillPixelsInPools(PlaneBlock &planeBlock) {
 115         ASSERT( !pools.empty() ); // assuming initPools has already been called
 116 //      iterate over pool types
 117         PoolList::iterator end= pools.begin();
 118         while ( end != pools.end() ) {
 119                 PoolList::iterator begin= end;
 120                 char type= begin->type;
 121         //      find the end of the same-pool-type block and invalidate the summers on the way
 122                 while ( end!=pools.end() && end->type==type ) {
 123                         end->summers_invalidate();
 124                         ++end;
 125                 }
 126
 127         //      we've got the interval, find out about the type
 128                 if (type!=DomPortion_Diamond) {
 129                 //      non-diamond domains all behave similarly
 130                         void (*shrinkProc)( CSMatrix , SMatrix , int, int );
 131                         switch (type) {
 132                                 case DomPortion_Standard:       shrinkProc= &shrinkToHalf;              break;
 133                                 case DomPortion_Horiz:          shrinkProc= &shrinkHorizontally;break;
 134                                 case DomPortion_Vert:           shrinkProc= &shrinkVertically;  break;
 135                                 default: ASSERT(false), shrinkProc=0;
 136                         }
 137                 //      we have the right procedure -> fill the first pool
 138                         ASSERT( begin->level == 1 );
 139                         shrinkProc( planeBlock.pixels, begin->pixels, begin->width, begin->height );
 140                 //      fill the rest (in the same-type interval)
 141                         while (++begin != end) {
 142                                 ASSERT( halfShrinkOK(begin-1,begin,zoom) );
 143                                 shrinkToHalf( (begin-1)->pixels, begin->pixels, begin->width, begin->height );
 144                         }
 145
 146                 } else { //     handle diamond-type domains
 147                         PoolList::iterator it= begin; //< the currently filled domain pool
 148                 //      fill the first set of diamond-type domain pools
 149                         bool horiz= width>=height;
 150                         int shift= lShift( getDiamondShift(min(width,height)), zoom );
 151                         int longerEnd= lShift( max(width,height)-minSizeNeededForDiamond(), zoom );
 152                         CSMatrix source= planeBlock.pixels;
 153
 154                         for (int l=0; l<=longerEnd; ++it,l+=shift) {
 155                                 ASSERT( it!=end && it->level==1 && it->width==it->height );
 156                                 shrinkToDiamond( planeBlock.pixels, it->pixels, it->width );
 157                                 source.shiftMatrix( (horiz?shift:0), (horiz?0:shift) );
 158                         }
 159                 //      now fill the multiscaled diamond pools
 160                         while (it!=end) {
 161                         //      too small pools are skipped
 162                                 while ( min(begin->width,begin->height) < 2*MinDomSize )
 163                                         ++begin;
 164                                 ASSERT( halfShrinkOK(begin,it,zoom) );
 165                                 shrinkToHalf( begin->pixels, it->pixels, it->width, it->height );
 166                         //      move on
 167                                 ++it;
 168                                 ++begin;
 169                         }
 170                 }//     if non-diamond else diamond
 171
 172         //      we just handled the whole interval (of diamond or other type)
 173
 174         }//     for (iterate over single-type intervals)
 175
 176 //      (cancelled) we filled all pools, let's prepare the summers
 177         //for_each( pools, mem_fun_ref(&Pool::summers_makeValid) );
 178 }//     ::fillPixelsInPools
 179
 180 namespace NOSPACE {
 181         /** Computes the ideal domain density for pool, level and max.\ domain count,
 182          *      the density is push_back-ed, returns the generated domain count (used once)
 183          *      \relates MStdDomains */
 184         inline static int bestDomainDensity( PoolIt pool, int level, int maxCount, int zoom
 185         , vector<short> &result) {
 186                 level-= zoom;
 187                 ASSERT( maxCount>=0 && level>0 && zoom>=0 );
 188                 int wms= rShift<int>(pool->width,zoom) -powers[level];
 189                 int hms= rShift<int>(pool->height,zoom) -powers[level];
 190         //      check whether any domain can fit and whether we should generate any more
 191                 if ( wms<0 || hms<0 || !maxCount ) {
 192                         result.push_back(0);
 193                         return 0;
 194                 }
 195                 int dens;
 196                 if (maxCount>1) {
 197                         Real temp= (  wms+hms+sqrt( sqr<Real>(wms+hms) +Real(4*wms*hms)*(maxCount-1) )  )
 198                                 / ( 2*(maxCount-1) );
 199                         dens= (int)ceil(temp);
 200                 } else
 201                         dens= 1+max(wms,hms);
 202
 203                 if (!dens)
 204                         dens= 1;
 205                 ASSERT(dens>0);
 206                 int count= (wms/dens+1)*(hms/dens+1);
 207                 ASSERT(count<=maxCount);
 208
 209                 result.push_back(dens);
 210                 return count;
 211         }
 212         /** Generates (\p results.push_back) densities for domains on level \p level for
 213          *      all pools of one type. It distributes at most \p maxCount domains among
 214          *      the pools' scale-levels in one of three ways (\p divType)
 215          *      and returns the number of generated domains \relates MStdDomains */
 216         static int divideDomsInType( PoolIt begin, PoolIt end, int maxCount, int level
 217         , char divType, int zoom, vector<short> &results ) {
 218                 ASSERT( begin!=end && divType>=0 && divType<=2 );
 219                 int scaleLevels= (end-1)->level - begin->level + 1;
 220                 int genCount= 0;
 221         //      iterate over same-scaleLevel intervals
 222                 for (; begin!=end; --scaleLevels) {
 223                         PoolIt it= begin+1;
 224                         while ( it!=end && it->level==begin->level )
 225                                 ++it;
 226                 //      we have the same-scale interval, find out how many domains to generate for it
 227                         int toGenerate= divType==2
 228                         //      half per scale level
 229                                 ? (maxCount-genCount)/2
 230                         //      uniform dividing or no multiscaling (then scaleLevels==1)
 231                                 : (maxCount-genCount)/scaleLevels ;
 232                 //      distribute it uniformly among the interval
 233                 //      (there are more than one for diamond-type only)
 234                         genCount+= toGenerate;  // genCount: "assume" we generate exactly toGenerate
 235                         for (; begin!=it; ++begin)
 236                                 toGenerate-=
 237                                         bestDomainDensity( begin, level, toGenerate/(it-begin), zoom, results );
 238                         genCount-= toGenerate;  // genCount: correct the count
 239                 }
 240                 return genCount;
 241         }
 242 }
 243 vector<short> MStdDomains::getLevelDensities(int level,int stdDomCountLog2) {
 244         ASSERT(level>=2);
 245 //      compute the sum of shares, check for no-domain situations
 246         int totalShares= settingsInt(DomPortion_Standard) + settingsInt(DomPortion_Horiz)
 247         + settingsInt(DomPortion_Vert) + settingsInt(DomPortion_Diamond);
 248         stdDomCountLog2-= (level-zoom-2)*settingsInt(MaxDomCountLevelDivisor);
 249
 250         if ( pools.empty() || !totalShares || stdDomCountLog2<0 )
 251                 return vector<short>( pools.size(), 0 );
 252 //      get the real domain count for this level
 253         int domCountLeft= powers[stdDomCountLog2];
 254
 255         vector<short> result;
 256         result.reserve(pools.size());
 257         PoolList::iterator begin, end;
 258         end= pools.begin();
 259 //      iterate over single-type domain-pool intervals
 260         while ( end!=pools.end() ) {
 261                 begin= end;
 262         //      find the end of the current interval
 263                 while ( end!=pools.end() && begin->type==end->type )
 264                         ++end;
 265         //      we've got a single-type interval
 266                 int share= settingsInt( (Settings)begin->type );
 267                 domCountLeft-= divideDomsInType( begin, end, domCountLeft*share/totalShares
 268                         , level, settingsInt(MultiDownScaling), zoom, result );
 269                 totalShares-= share;
 270         }
 271 //      check we created the correct number of densities
 272         ASSERT( result.size() == pools.size() );
 273         return result;
 274 }
 275
 276 void MStdDomains::writeSettings(ostream &file) {
 277         ASSERT( this && settings );
 278 //      all settings are small integers a need to be preserved
 279         int setLength= info().setLength;
 280         for (int i=0; i<setLength; ++i)
 281                 put<Uchar>( file, settingsInt(i) );
 282 }
 283
 284 void MStdDomains::readSettings(istream &file) {
 285         ASSERT( this && settings );
 286 //      all settings are small integers a need to be preserved
 287         int setLength= info().setLength;
 288         for (int i=0; i<setLength; ++i)
 289                 settingsInt(i)= get<Uchar>(file);
 290 }