- Main Page
- Namespaces
- Classes
- Files
- Directories

- File List
- File Members

# comp_compress.cpp

    00001 /* ***** BEGIN LICENSE BLOCK *****
    00002 *
    00003 * $Id: comp_compress.cpp,v 1.2 2005/01/30 05:11:41 gabest Exp $ $Name:  $
    00004 *
    00005 * Version: MPL 1.1/GPL 2.0/LGPL 2.1
    00006 *
    00007 * The contents of this file are subject to the Mozilla Public License
    00008 * Version 1.1 (the "License"); you may not use this file except in compliance
    00009 * with the License. You may obtain a copy of the License at
    00010 * http://www.mozilla.org/MPL/
    00011 *
    00012 * Software distributed under the License is distributed on an "AS IS" basis,
    00013 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License for
    00014 * the specific language governing rights and limitations under the License.
    00015 *
    00016 * The Original Code is BBC Research and Development code.
    00017 *
    00018 * The Initial Developer of the Original Code is the British Broadcasting
    00019 * Corporation.
    00020 * Portions created by the Initial Developer are Copyright (C) 2004.
    00021 * All Rights Reserved.
    00022 *
    00023 * Contributor(s): Thomas Davies (Original Author),
    00024 *                 Scott R Ladd,
    00025 *                 Anuradha Suraparaju,
    00026 *                 Peter Meerwald
    00027 *
    00028 * Alternatively, the contents of this file may be used under the terms of
    00029 * the GNU General Public License Version 2 (the "GPL"), or the GNU Lesser
    00030 * Public License Version 2.1 (the "LGPL"), in which case the provisions of
    00031 * the GPL or the LGPL are applicable instead of those above. If you wish to
    00032 * allow use of your version of this file only under the terms of the either
    00033 * the GPL or LGPL and not to allow others to use your version of this file
    00034 * under the MPL, indicate your decision by deleting the provisions above
    00035 * and replace them with the notice and other provisions required by the GPL
    00036 * or LGPL. If you do not delete the provisions above, a recipient may use
    00037 * your version of this file under the terms of any one of the MPL, the GPL
    00038 * or the LGPL.
    00039 * ***** END LICENSE BLOCK ***** */
    00040
    00041
    00042 //Compression of an individual component,
    00043 //after motion compensation if appropriate
    00045
    00046 #include <libdirac_encoder/comp_compress.h>
    00047 #include <libdirac_common/band_codec.h>
    00048 #include <libdirac_common/golomb.h>
    00049 using namespace dirac;
    00050
    00051 #include <ctime>
    00052 #include <vector>
    00053 #include <iostream>
    00054
    00055 using std::log;
    00056 using std::floor;
    00057
    00058 static inline double pow4 (double x)
    00059 {
    00060     return x * x * x* x;
    00061 }
    00062
    00063 CompCompressor::CompCompressor( EncoderParams& encp,const FrameParams& fp)
    00064 : m_encparams(encp),
    00065   m_fparams(fp),
    00066   m_fsort( m_fparams.FSort() ),
    00067   m_cformat( m_fparams.CFormat() ),
    00068   m_qflist(60),
    00069   m_qfinvlist(60),
    00070   m_offset(60)
    00071 {}
    00072
    00073 void CompCompressor::Compress(PicArray& pic_data)
    00074 {
    00075
    00076     //need to transform, select quantisers for each band, and then compress each component in turn
    00077     m_csort=pic_data.CSort();
    00078     const int depth=4;
    00079     unsigned int num_band_bits;
    00080
    00081     // A pointer to an object  for coding the subband data
    00082     BandCodec* bcoder;
    00083
    00084     // A pointer to an object for outputting the subband data
    00085     UnitOutputManager* band_op;
    00086
    00087     const size_t CONTEXTS_REQUIRED = 24;
    00088
    00089     Subband node;
    00090
    00091     //set up Lagrangian params
    00092     if (m_fsort == I_frame)
    00093         m_lambda= m_encparams.ILambda();
    00094     else if (m_fsort == L1_frame)
    00095         m_lambda= m_encparams.L1Lambda();
    00096     else
    00097         m_lambda= m_encparams.L2Lambda();
    00098
    00099      if (m_csort == U_COMP)
    00100          m_lambda*= m_encparams.UFactor();
    00101      if (m_csort == V_COMP)
    00102          m_lambda*= m_encparams.VFactor();
    00103
    00104     WaveletTransform wtransform(depth);
    00105
    00106     wtransform.Transform( FORWARD , pic_data );
    00107     wtransform.SetBandWeights( m_encparams.CPD() , m_fparams.FSort() , m_fparams.CFormat(), m_csort);
    00108
    00109     SubbandList& bands=wtransform.BandList();
    00110
    00111     // Generate all the quantisation data
    00112     GenQuantList();
    00113
    00114     // Choose all the quantisers
    00115     OneDArray<unsigned int> estimated_bits( Range( 1 , bands.Length() ) );
    00116     SelectQuantisers( pic_data , bands , estimated_bits );
    00117
    00118     // Loop over all the bands (from DC to HF) quantising and coding them
    00119     for (int b=bands.Length() ; b>=1 ; --b )
    00120     {
    00121         band_op = & m_encparams.BitsOut().FrameOutput().BandOutput( m_csort , b );
    00122
    00123         GolombCode( band_op->Header() , bands(b).Qf(0) );
    00124
    00125         if (bands(b).Qf(0) != -1)
    00126         {   // If not skipped ...
    00127
    00128             bands(b).SetQf( 0 , m_qflist[bands(b).Qf(0)] );
    00129
    00130              // Pick the right codec according to the frame type and subband
    00131             if (b >= bands.Length())
    00132             {
    00133                 if ( m_fsort == I_frame && b == bands.Length() )
    00134                     bcoder=new IntraDCBandCodec( &( band_op->Data() ) , CONTEXTS_REQUIRED , bands);
    00135                 else
    00136                     bcoder=new LFBandCodec( &( band_op->Data() ) ,CONTEXTS_REQUIRED, bands , b);
    00137             }
    00138             else
    00139                 bcoder=new BandCodec( &( band_op->Data() ) , CONTEXTS_REQUIRED , bands , b);
    00140
    00141             num_band_bits = bcoder->Compress(pic_data);
    00142
    00143              // Update the entropy correction factors
    00144             m_encparams.EntropyFactors().Update(b , m_fsort , m_csort , estimated_bits[b] , num_band_bits);
    00145
    00146             // Write the length of the data chunk into the header, and flush everything out to file
    00147             UnsignedGolombCode( band_op->Header() , num_band_bits);
    00148
    00149             delete bcoder;
    00150         }
    00151         else
    00152         {   // ... skipped
    00153
    00154             if (b == bands.Length() && m_fsort == I_frame)
    00155                 SetToVal( pic_data , bands(b) , 2692 );
    00156             else
    00157                 SetToVal( pic_data , bands(b) , 0 );
    00158         }
    00159     }//b
    00160
    00161     // Transform back into the picture domain
    00162     wtransform.Transform( BACKWARD , pic_data );
    00163
    00164 }
    00165
    00166 void CompCompressor::GenQuantList()
    00167 {    //generates the list of quantisers and inverse quantisers
    00168     //there is some repetition in this list but at the moment this is easiest from the perspective of SelectQuant
    00169     //Need to remove this repetition later
    00170
    00171     m_qflist[0]=1;        m_qfinvlist[0]=131072;    m_offset[0]=0;
    00172     m_qflist[1]=1;        m_qfinvlist[1]=131072;    m_offset[1]=0;
    00173     m_qflist[2]=1;        m_qfinvlist[2]=131072;    m_offset[2]=0;
    00174     m_qflist[3]=1;        m_qfinvlist[3]=131072;    m_offset[3]=0;
    00175     m_qflist[4]=2;        m_qfinvlist[4]=65536;        m_offset[4]=1;
    00176     m_qflist[5]=2;        m_qfinvlist[5]=65536;        m_offset[5]=1;
    00177     m_qflist[6]=2;        m_qfinvlist[6]=65536;        m_offset[6]=1;
    00178     m_qflist[7]=3;        m_qfinvlist[7]=43690;        m_offset[7]=1;
    00179     m_qflist[8]=4;        m_qfinvlist[8]=32768;        m_offset[8]=2;
    00180     m_qflist[9]=4;        m_qfinvlist[9]=32768;        m_offset[9]=2;
    00181     m_qflist[10]=5;        m_qfinvlist[10]=26214;    m_offset[10]=2;
    00182     m_qflist[11]=6;        m_qfinvlist[11]=21845;    m_offset[11]=2;
    00183     m_qflist[12]=8;        m_qfinvlist[12]=16384;    m_offset[12]=3;
    00184     m_qflist[13]=9;        m_qfinvlist[13]=14563;    m_offset[13]=3;
    00185     m_qflist[14]=11;        m_qfinvlist[14]=11915;    m_offset[14]=4;
    00186     m_qflist[15]=13;        m_qfinvlist[15]=10082;    m_offset[15]=5;
    00187     m_qflist[16]=16;        m_qfinvlist[16]=8192;        m_offset[16]=6;
    00188     m_qflist[17]=19;        m_qfinvlist[17]=6898;        m_offset[17]=7;
    00189     m_qflist[18]=22;        m_qfinvlist[18]=5957;        m_offset[18]=8;
    00190     m_qflist[19]=26;        m_qfinvlist[19]=5041;        m_offset[19]=10;
    00191     m_qflist[20]=32;        m_qfinvlist[20]=4096;        m_offset[20]=12;
    00192     m_qflist[21]=38;        m_qfinvlist[21]=3449;        m_offset[21]=14;
    00193     m_qflist[22]=45;        m_qfinvlist[22]=2912;        m_offset[22]=17;
    00194     m_qflist[23]=53;        m_qfinvlist[23]=2473;        m_offset[23]=20;
    00195     m_qflist[24]=64;        m_qfinvlist[24]=2048;        m_offset[24]=24;
    00196     m_qflist[25]=76;        m_qfinvlist[25]=1724;        m_offset[25]=29;
    00197     m_qflist[26]=90;        m_qfinvlist[26]=1456;        m_offset[26]=34;
    00198     m_qflist[27]=107;        m_qfinvlist[27]=1224;        m_offset[27]=40;
    00199     m_qflist[28]=128;        m_qfinvlist[28]=1024;        m_offset[28]=48;
    00200     m_qflist[29]=152;        m_qfinvlist[29]=862;        m_offset[29]=57;
    00201     m_qflist[30]=181;        m_qfinvlist[30]=724;        m_offset[30]=68;
    00202     m_qflist[31]=215;        m_qfinvlist[31]=609;        m_offset[31]=81;
    00203     m_qflist[32]=256;        m_qfinvlist[32]=512;        m_offset[32]=96;
    00204     m_qflist[33]=304;        m_qfinvlist[33]=431;        m_offset[33]=114;
    00205     m_qflist[34]=362;        m_qfinvlist[34]=362;        m_offset[34]=136;
    00206     m_qflist[35]=430;        m_qfinvlist[35]=304;        m_offset[35]=161;
    00207     m_qflist[36]=512;        m_qfinvlist[36]=256;        m_offset[36]=192;
    00208     m_qflist[37]=608;        m_qfinvlist[37]=215;        m_offset[37]=228;
    00209     m_qflist[38]=724;        m_qfinvlist[38]=181;        m_offset[38]=272;
    00210     m_qflist[39]=861;        m_qfinvlist[39]=152;        m_offset[39]=323;
    00211     m_qflist[40]=1024;    m_qfinvlist[40]=128;        m_offset[40]=384;
    00212     m_qflist[41]=1217;    m_qfinvlist[41]=107;        m_offset[41]=456;
    00213     m_qflist[42]=1448;    m_qfinvlist[42]=90;        m_offset[42]=543;
    00214     m_qflist[43]=1722;    m_qfinvlist[43]=76;        m_offset[43]=646;
    00215     m_qflist[44]=2048;    m_qfinvlist[44]=64;        m_offset[44]=768;
    00216     m_qflist[45]=2435;    m_qfinvlist[45]=53;        m_offset[45]=913;
    00217     m_qflist[46]=2896;    m_qfinvlist[46]=45;        m_offset[46]=1086;
    00218     m_qflist[47]=3444;    m_qfinvlist[47]=38;        m_offset[47]=1292;
    00219     m_qflist[48]=4096;    m_qfinvlist[48]=32;        m_offset[48]=1536;
    00220     m_qflist[49]=4870;    m_qfinvlist[49]=26;        m_offset[49]=1826;
    00221     m_qflist[50]=5792;    m_qfinvlist[50]=22;        m_offset[50]=2172;
    00222     m_qflist[51]=6888;    m_qfinvlist[51]=19;        m_offset[51]=2583;
    00223     m_qflist[52]=8192;    m_qfinvlist[52]=16;        m_offset[52]=3072;
    00224     m_qflist[53]=9741;    m_qfinvlist[53]=13;        m_offset[53]=3653;
    00225     m_qflist[54]=11585;    m_qfinvlist[54]=11;        m_offset[54]=4344;
    00226     m_qflist[55]=13777;    m_qfinvlist[55]=9;        m_offset[55]=5166;
    00227     m_qflist[56]=16384;    m_qfinvlist[56]=8;        m_offset[56]=6144;
    00228     m_qflist[57]=19483;    m_qfinvlist[57]=6;        m_offset[57]=7306;
    00229     m_qflist[58]=23170;    m_qfinvlist[58]=5;        m_offset[58]=8689;
    00230     m_qflist[59]=27554;    m_qfinvlist[59]=4;        m_offset[59]=10333;
    00231 }
    00232
    00233
    00234
    00235 void CompCompressor::SelectQuantisers( PicArray& pic_data , SubbandList& bands ,
    00236                                                  OneDArray<unsigned int>& est_counts )
    00237 {
    00238     // Select all the quantizers
    00239     for ( int b=bands.Length() ; b>=1 ; --b )
    00240         est_counts[b] = SelectQuant( pic_data , bands , b );
    00241 }
    00242
    00243 int CompCompressor::SelectQuant(PicArray& pic_data,SubbandList& bands,const int band_num)
    00244 {
    00245
    00246     Subband& node=bands(band_num);
    00247
    00248     const int qf_start_idx = 4;
    00249
    00250     if (band_num==bands.Length())
    00251         AddSubAverage(pic_data,node.Xl(),node.Yl(),SUBTRACT);
    00252
    00253     int min_idx;
    00254     double bandmax=PicAbsMax(pic_data,node.Xp(),node.Yp(),node.Xl(),node.Yl());
    00255
    00256     if (bandmax>=1)
    00257         node.SetMax(int(floor(log(float(bandmax))/log(2.0))));
    00258     else
    00259         node.SetMax(0);
    00260     int length=4*node.Max()+5;//this is the number of quantisers that are possible
    00261
    00262     OneDArray<int> count0(length);
    00263     int count1;
    00264     OneDArray<int> countPOS(length);
    00265     OneDArray<int> countNEG(length);
    00266     OneDArray<double> error_total(length);
    00267     OneDArray<CostType> costs(length);
    00268     int quant_val;
    00269     ValueType val,abs_val;
    00270     int error;
    00271     double p0,p1;
    00272     double sign_entropy;
    00273
    00274     int xp=node.Xp();
    00275     int yp=node.Yp();
    00276     int xl=node.Xl();
    00277     int yl=node.Yl();
    00278     double vol;
    00279
    00280     if (bandmax < 1.0 )
    00281     {
    00282         //coefficients are zero so the subband can be skipped
    00283         node.SetQf(0,-1);//indicates that the subband is skipped
    00284
    00285         if ( band_num == bands.Length() )
    00286             AddSubAverage(pic_data,node.Xl(),node.Yl(),ADD);
    00287
    00288         return 0;
    00289     }
    00290     else
    00291     {
    00292         for ( int q=0 ; q<costs.Length() ; q++)
    00293         {
    00294             error_total[q] = 0.0;
    00295             count0[q] = 0;
    00296             countPOS[q] = 0;
    00297             countNEG[q] = 0;
    00298         }
    00299
    00300         //first, find to nearest integral number of bits using 1/4 of the data
    00302         vol=double((yl/2)*(xl/2));//vol is only 1/4 of the coeffs
    00303         count1=int(vol);
    00304         for ( int j=yp+1 ; j<yp+yl ; j+=2 )
    00305         {
    00306             for ( int i=xp+((j-yp)%4)/2 ; i<xp+xl ; i+=2)
    00307             {
    00308
    00309                 val = pic_data[j][i];
    00310                 quant_val = abs(val);
    00311                 abs_val = quant_val;
    00312
    00313                 for ( int q=qf_start_idx ; q<costs.Length() ; q+=4)
    00314                 {
    00315                     quant_val >>= (q/4);
    00316
    00317                     if (quant_val)
    00318                     {
    00319                         count0[q]+=quant_val;
    00320                         quant_val <<= (q/4);
    00321                         if (val>0)
    00322                             countPOS[q]++;
    00323                         else
    00324                             countNEG[q]++;
    00325                     }
    00326
    00327                     error = abs_val-quant_val;
    00328
    00329                     if ( quant_val != 0)
    00330                         error -= m_offset[q];
    00331
    00332                     error_total[q] +=  pow4( static_cast<double>(error) );
    00333                 }// q
    00334             }// i
    00335         }// j
    00336
    00337          //do entropy calculation etc
    00338         for ( int q=qf_start_idx ; q<costs.Length() ; q+=4 )
    00339         {
    00340             costs[q].MSE = error_total[q]/( vol*node.Wt()*node.Wt() );
    00341 //
    00342             costs[q].MSE = std::sqrt( costs[q].MSE );
    00343 //
    00344              //calculate probabilities and entropy
    00345             p0 = double( count0[q] )/double( count0[q]+count1 );
    00346             p1 = 1.0-p0;
    00347
    00348             if ( p0 != 0.0 && p1 != 0.0)
    00349                 costs[q].ENTROPY =- (p0*log(p0)+p1*log(p1))/log(2.0);
    00350             else
    00351                 costs[q].ENTROPY = 0.0;
    00352
    00353             //we want the entropy *per symbol*, not per bit ...
    00354             costs[q].ENTROPY *= double(count0[q]+count1);
    00355             costs[q].ENTROPY /= vol;
    00356
    00357             //now add in the sign entropy
    00358             if ( countPOS[q]+countNEG[q] != 0 )
    00359             {
    00360                 p0 = float(countNEG[q])/float( countPOS[q]+countNEG[q] );
    00361                 p1 = 1.0-p0;
    00362                 if ( p0 != 0.0 && p1 != 0.0)
    00363                     sign_entropy = -( (p0*log(p0)+p1*log(p1) ) / log(2.0));
    00364                 else
    00365                     sign_entropy = 0.0;
    00366             }
    00367             else
    00368                 sign_entropy=0.0;
    00369
    00370               //we want the entropy *per symbol*, not per bit ...
    00371             sign_entropy *= double( countNEG[q]+countPOS[q] );
    00372             sign_entropy /= vol;
    00373
    00374             costs[q].ENTROPY += sign_entropy;
    00375
    00376             //sort out correction factors
    00377             costs[q].ENTROPY *= m_encparams.EntropyFactors().Factor(band_num,m_fsort,m_csort);
    00378             costs[q].TOTAL = costs[q].MSE+m_lambda*costs[q].ENTROPY;
    00379
    00380         }// q
    00381
    00382         //find the qf with the lowest cost
    00383         min_idx=qf_start_idx;
    00384         for ( int q=qf_start_idx ; q<costs.Length() ; q+=4 )
    00385         {
    00386             if ( costs[q].TOTAL < costs[min_idx].TOTAL )
    00387                 min_idx=q;
    00388         }
    00389
    00390         //now repeat to get to 1/2 bit accuracy
    00392         for ( int q=std::max(0,min_idx-2) ; q<=std::min(costs.Last(),min_idx+2) ; q+=2 )
    00393         {
    00394             if ( q != min_idx )
    00395             {
    00396                 error_total[q] = 0.0;
    00397                 count0[q] = 0;
    00398                 countPOS[q] = 0;
    00399                 countNEG[q] = 0;
    00400             }
    00401         }
    00402
    00403         vol = double( (yl/2) * (xl/2) );
    00404         count1 = int(vol);
    00405         int top_idx = std::min(costs.Last(),min_idx+2);
    00406         int bottom_idx = std::max(0,min_idx-2);
    00407
    00408         for (int j=yp+1 ; j<yp+yl ; j+=2 )
    00409         {
    00410             for (int i=xp+1 ; i<xp+xl ; i+=2 )
    00411             {
    00412                 val = pic_data[j][i];
    00413                 abs_val = abs(val);
    00414
    00415                 for ( int q=bottom_idx ; q<=top_idx ; q+=2 )
    00416                 {
    00417                     if ( q != min_idx )
    00418                     {
    00419                         quant_val = int(abs_val);
    00420                         quant_val *= m_qfinvlist[q];
    00421                         quant_val >>= 17;
    00422
    00423                         if ( quant_val )
    00424                         {
    00425                             count0[q] += quant_val;
    00426                             quant_val *= m_qflist[q];
    00427
    00428                             if (val>0.0)
    00429                                 countPOS[q]++;
    00430                             else
    00431                                 countNEG[q]++;
    00432                         }
    00433
    00434                         error = abs_val-quant_val;
    00435
    00436                         if ( quant_val != 0 )
    00437                             error -= m_offset[q];
    00438
    00439                         error_total[q] += pow4( static_cast<double>(error) );
    00440                     }//end of if
    00441                 }//q
    00442             }//J
    00443         }//I
    00444
    00445          //do entropy calculation
    00446         for ( int q=bottom_idx ; q<=top_idx ; q+=2 )
    00447         {
    00448             if ( q != min_idx )
    00449             {
    00450                 costs[q].MSE = error_total[q] / (vol*node.Wt()*node.Wt());
    00451 //
    00452                 costs[q].MSE = std::sqrt( costs[q].MSE );
    00453 //
    00454
    00455
    00456                   //calculate probabilities and entropy
    00457                 p0 = double(count0[q]) / double(count0[q]+count1);
    00458                 p1 = 1.0-p0;
    00459
    00460                 if (p0 != 0.0 && p1 != 0.0)
    00461                     costs[q].ENTROPY =- (p0*log(p0) + p1*log(p1)) / log(2.0);
    00462                 else
    00463                     costs[q].ENTROPY = 0.0;
    00464                 //we want the entropy *per symbol*, not per bit ...
    00465                 costs[q].ENTROPY *= count0[q]+count1;
    00466                 costs[q].ENTROPY /= vol;
    00467
    00468                   //now add in the sign entropy
    00469                 if (countPOS[q]+countNEG[q] != 0)
    00470                 {
    00471                     p0 = double( countNEG[q] )/double( countPOS[q] + countNEG[q] );
    00472                     p1 = 1.0-p0;
    00473                     if (p0 != 0.0 && p1 != 0.0)
    00474                         sign_entropy =- ( (p0*log(p0)+p1*log(p1)) / log(2.0) );
    00475                     else
    00476                         sign_entropy = 0.0;
    00477                 }
    00478                 else
    00479                     sign_entropy = 0.0;
    00480
    00481                   //we want the entropy *per symbol*, not per bit ...
    00482                 sign_entropy *= double(countNEG[q]+countPOS[q]);
    00483                 sign_entropy /= vol;
    00484
    00485                 costs[q].ENTROPY += sign_entropy;
    00486                 //sort out correction factors
    00487                 costs[q].ENTROPY *= m_encparams.EntropyFactors().Factor(band_num,m_fsort,m_csort);
    00488                 costs[q].TOTAL = costs[q].MSE+m_lambda*costs[q].ENTROPY;
    00489             }
    00490         }//q
    00491
    00492          //find the qf with the lowest cost
    00493         for ( int q=bottom_idx ; q<=top_idx ; q+=2 )
    00494         {
    00495             if ( costs[q].TOTAL < costs[min_idx].TOTAL )
    00496                 min_idx = q;
    00497         }
    00498
    00499          //finally use 1/2 the values to get 1/4 bit accuracy
    00501
    00502         bottom_idx=std::max(0,min_idx-1);
    00503         top_idx=std::min(costs.Length()-1,min_idx+1);
    00504
    00505         for ( int q=bottom_idx ; q<=top_idx ; q++ )
    00506         {
    00507             error_total[q] = 0.0;
    00508             count0[q] = 0;
    00509             countPOS[q] = 0;
    00510             countNEG[q] = 0;
    00511         }
    00512
    00513         vol = double( (yl/2) * xl );
    00514         count1 = int( vol );
    00515
    00516         for (int j=yp ; j<yp+yl ; ++j )
    00517         {
    00518             for (int i=xp+1 ; i<xp+xl ; i+=2)
    00519             {
    00520                 val = pic_data[j][i];
    00521                 abs_val = abs(val);
    00522
    00523                 for (int q=bottom_idx;q<=top_idx;q++)
    00524                 {
    00525                     quant_val = int(abs_val);
    00526                     quant_val *= m_qfinvlist[q];
    00527                     quant_val >>= 17;
    00528
    00529                     if (quant_val)
    00530                     {
    00531                         count0[q] += quant_val;
    00532                         quant_val *= m_qflist[q];
    00533
    00534                         if ( val > 0 )
    00535                             countPOS[q]++;
    00536                         else
    00537                             countNEG[q]++;
    00538                     }
    00539
    00540                     error = abs_val - quant_val;
    00541
    00542                     if ( quant_val != 0 )
    00543                         error -= m_offset[q];
    00544
    00545                     error_total[q] +=  pow4( static_cast<double>(error) );
    00546                 }//q
    00547             }//i
    00548         }//j
    00549
    00550          //do entropy calculation
    00551         for ( int q=bottom_idx ; q<=top_idx ; q++ )
    00552         {
    00553             costs[q].MSE = error_total[q]/(vol*node.Wt()*node.Wt());
    00554 //
    00555             costs[q].MSE = std::sqrt( costs[q].MSE );
    00556 //
    00557              //calculate probabilities and entropy
    00558             p0 = double( count0[q] )/ double( count0[q]+count1 );
    00559             p1 = 1.0 - p0;
    00560
    00561             if ( p0 != 0.0 && p1 != 0.0)
    00562                 costs[q].ENTROPY = -( p0*log(p0)+p1*log(p1) ) / log(2.0);
    00563             else
    00564                 costs[q].ENTROPY = 0.0;
    00565
    00566               //we want the entropy *per symbol*, not per bit ...
    00567             costs[q].ENTROPY *= double(count0[q]+count1);
    00568             costs[q].ENTROPY /= vol;
    00569
    00570              //now add in the sign entropy
    00571             if ( countPOS[q] + countNEG[q] != 0 )
    00572             {
    00573                 p0 = double( countNEG[q] )/double( countPOS[q]+countNEG[q] );
    00574                 p1 = 1.0-p0;
    00575                 if ( p0 != 0.0 && p1 != 0.0)
    00576                     sign_entropy = -( (p0*log(p0)+p1*log(p1) ) / log(2.0));
    00577                 else
    00578                     sign_entropy = 0.0;
    00579             }
    00580             else
    00581                 sign_entropy = 0.0;
    00582
    00583               //we want the entropy *per symbol*, not per bit ...
    00584             sign_entropy *= double(countNEG[q]+countPOS[q]);
    00585             sign_entropy /= vol;
    00586
    00587             costs[q].ENTROPY += sign_entropy;
    00588
    00589             //sort out correction factors
    00590             costs[q].ENTROPY *= m_encparams.EntropyFactors().Factor(band_num,m_fsort,m_csort);
    00591             costs[q].TOTAL = costs[q].MSE+m_lambda*costs[q].ENTROPY;
    00592
    00593         }//q
    00594
    00595          //find the qf with the lowest cost
    00596         for ( int q=bottom_idx ; q<=top_idx ; q++ )
    00597         {
    00598             if ( costs[q].TOTAL < costs[min_idx].TOTAL )
    00599                 min_idx=q;
    00600         }
    00601
    00602         if ( costs[min_idx].ENTROPY == 0.0 )//then can skip after all
    00603             node.SetQf(0,-1);
    00604         else
    00605             node.SetQf(0,min_idx);
    00606
    00607         if ( band_num == bands.Length())
    00608             AddSubAverage(pic_data,node.Xl(),node.Yl(),ADD);
    00609
    00610         return int(costs[min_idx].ENTROPY*double(xl*yl));
    00611     }
    00612
    00613 }
    00614
    00615 ValueType CompCompressor::PicAbsMax(const PicArray& pic_data) const
    00616 {
    00617     //finds the maximum absolute value of the picture array
    00618     return PicAbsMax(pic_data,pic_data.FirstX() , pic_data.FirstY(),
    00619                      pic_data.LengthX(),pic_data.LengthY());
    00620 }
    00621
    00622 ValueType CompCompressor::PicAbsMax(const PicArray& pic_data,int xp, int yp ,int xl ,int yl) const
    00623 {
    00624
    00625     int first_x=std::max(pic_data.FirstX(),xp);
    00626     int first_y=std::max(pic_data.FirstY(),yp);
    00627     int last_x=std::min(pic_data.LastX(),xp+xl-1);
    00628     int last_y=std::min(pic_data.LastY(),yp+yl-1);
    00629     ValueType val=0;
    00630
    00631     for (int j=first_y ; j<=last_y; ++j)
    00632     {
    00633         for (int i=first_x ; i<=last_x; ++i)
    00634         {
    00635             val = std::max( val , pic_data[j][i] );
    00636         }// i
    00637     }// j
    00638
    00639     return val;
    00640 }
    00641
    00642 void CompCompressor::SetToVal(PicArray& pic_data,const Subband& node,ValueType val){
    00643
    00644     for (int j=node.Yp() ; j<node.Yp() + node.Yl() ; ++j)
    00645     {
    00646         for (int i=node.Xp(); i<node.Xp() + node.Xl() ; ++i)
    00647         {
    00648             pic_data[j][i] = val;
    00649         }// i
    00650     }// j
    00651
    00652 }
    00653
    00654 void CompCompressor::AddSubAverage(PicArray& pic_data,int xl,int yl,AddOrSub dirn){
    00655
    00656     ValueType last_val=2692;//corresponds to mid-grey in this DC band with these filters
    00657                             //NB this is hard-wired for a level 4 transform
    00658     ValueType last_val2;
    00659
    00660     if ( dirn == SUBTRACT )
    00661     {
    00662         for ( int j=0 ; j<yl ; j++)
    00663             {
    00664             for ( int i=0 ; i<xl ; i++)
    00665                 {
    00666                 last_val2 = pic_data[j][i];
    00667                 pic_data[j][i] -= last_val;
    00668                 last_val = last_val2;
    00669             }// i
    00670         }// j
    00671     }
    00672     else
    00673     {
    00674         for ( int j=0 ; j<yl ; j++)
    00675         {
    00676             for ( int i=0 ; i<xl; i++ )
    00677             {
    00678                 pic_data[j][i] += last_val;
    00679                 last_val = pic_data[j][i];
    00680             }// i
    00681         }// j
    00682
    00683     }
    00684 }