libdirac_encoder/prefilter.cpp

   1 /* ***** BEGIN LICENSE BLOCK *****
   2 *
   3 * $Id: prefilter.cpp,v 1.3 2008/06/19 10:28:53 tjdwave Exp $ $Name:  $
   4 *
   5 * Version: MPL 1.1/GPL 2.0/LGPL 2.1
   6 *
   7 * The contents of this file are subject to the Mozilla Public License
   8 * Version 1.1 (the "License"); you may not use this file except in compliance
   9 * with the License. You may obtain a copy of the License at
  10 * http://www.mozilla.org/MPL/
  11 *
  12 * Software distributed under the License is distributed on an "AS IS" basis,
  13 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License for
  14 * the specific language governing rights and limitations under the License.
  15 *
  16 * The Original Code is BBC Research and Development code.
  17 *
  18 * The Initial Developer of the Original Code is the British Broadcasting
  19 * Corporation.
  20 * Portions created by the Initial Developer are Copyright (C) 2008.
  21 * All Rights Reserved.
  22 *
  23 * Contributor(s): Thomas Davies (Original Author),
  24 *
  25 * Alternatively, the contents of this file may be used under the terms of
  26 * the GNU General Public License Version 2 (the "GPL"), or the GNU Lesser
  27 * Public License Version 2.1 (the "LGPL"), in which case the provisions of
  28 * the GPL or the LGPL are applicable instead of those above. If you wish to
  29 * allow use of your version of this file only under the terms of the either
  30 * the GPL or LGPL and not to allow others to use your version of this file
  31 * under the MPL, indicate your decision by deleting the provisions above
  32 * and replace them with the notice and other provisions required by the GPL
  33 * or LGPL. If you do not delete the provisions above, a recipient may use
  34 * your version of this file under the terms of any one of the MPL, the GPL
  35 * or the LGPL.
  36 * ***** END LICENSE BLOCK ***** */
  37
  38 #include<libdirac_encoder/prefilter.h>
  39 #include<libdirac_common/arrays.h>
  40
  41 using namespace dirac;
  42
  43 void dirac::CWMFilter( Picture& picture, const int strength )
  44 {
  45     CWMFilterComponent( picture.Data(Y_COMP), strength );
  46     CWMFilterComponent( picture.Data(U_COMP), strength );
  47     CWMFilterComponent( picture.Data(V_COMP), strength );
  48 }
  49
  50 void dirac::CWMFilterComponent( PicArray& pic_data, const int strength )
  51 {
  52     // Do centre-weighted median denoising
  53
  54     PicArray pic_copy( pic_data );
  55
  56     const int width( 3 );
  57     const int offset( (width-1)/2 );
  58     const int centre_weight = std::max(1, (width*width+1)-strength );
  59     const int list_length = centre_weight+(width*width)-1;
  60
  61     ValueType* val_list = new ValueType[list_length];
  62
  63     for (int j=offset; j<pic_data.LengthY()-offset; ++j){
  64         for (int i=offset; i<pic_data.LastX()-offset; ++i){
  65             // Make the value list
  66             int pos=0;
  67             for (; pos<centre_weight-1; ++pos)
  68                 val_list[pos] = pic_copy[j][i];
  69
  70             for (int s=-offset; s<=offset; ++s){
  71                 for (int r=-offset; r<=offset; ++r){
  72                     val_list[pos]=pic_copy[j+s][i+r];
  73                     pos++;
  74                 }// r
  75             }// s
  76
  77             pic_data[j][i] = Median( val_list, list_length );
  78         }// i
  79     }// j
  80
  81     delete[] val_list;
  82 }
  83
  84 ValueType dirac::Median( const ValueType* val_list, const int length)
  85 {
  86
  87
  88     OneDArray<ValueType> ordered_vals( length );
  89
  90     // Place the values in order
  91     int pos=0;
  92     ordered_vals[0] = val_list[0];
  93     for (int i=1 ; i<length ; ++i )
  94     {
  95         for (int k=0 ; k<i ; ++k)
  96         {
  97             if (val_list[i]<ordered_vals[k])
  98             {
  99                 pos=k;
 100                 break;
 101             }
 102             else
 103                 pos=k+1;
 104         }// k
 105
 106         if ( pos==i)
 107             ordered_vals[i] = val_list[i];
 108         else
 109         {
 110             for (int k=i-1 ; k>=pos ; --k )
 111             {
 112                 ordered_vals[k+1] = ordered_vals[k];
 113             }// k
 114             ordered_vals[pos] = val_list[i];
 115         }
 116     }// i
 117
 118     // return the middle value
 119     if ( length%2!=0 )
 120         return ordered_vals[(length-1)/2];
 121     else
 122         return (ordered_vals[(length/2)-1]+ordered_vals[length/2]+1)>>1;
 123
 124 }
 125
 126 /*************************************************************/
 127
 128
 129 void VFilter( PicArray& pic_data, const OneDArray<int>& filter, const int bits );
 130 void HFilter( PicArray& pic_data, const OneDArray<int>& filter, const int bits );
 131
 132 double sinxoverx( const double val )
 133 {
 134     if ( 0.0f == val )
 135         return 1.0;
 136     else
 137         return sin(val)/val;
 138
 139 }
 140
 141 OneDArray<int> MakeLPRectFilter( const float bw, const int bits )
 142 {
 143     const int tl = 8;
 144     const float pi = 3.1415926535;
 145
 146     OneDArray<double> double_filter( Range( -tl, tl ) );
 147     OneDArray<int> int_filter( Range( -tl, tl) );
 148
 149     // Use the Hanning window
 150     for (int i=double_filter.First(); i<=double_filter.Last(); ++i)
 151     {
 152         double_filter[i] = cos( (pi*i)/
 153                                    (double_filter.Length()+1) );
 154     }
 155
 156     // Apply sinc function
 157     for (int i=double_filter.First(); i<=double_filter.Last(); ++i)
 158     {
 159         double_filter[i] *= sinxoverx( pi*1.0*bw*i );
 160     }
 161
 162     // Get DC gain = 1<<bits
 163     double sum = 0.0;
 164     for (int i=double_filter.First(); i<=double_filter.Last(); ++i)
 165         sum += double_filter[i];
 166
 167     for (int i=double_filter.First(); i<=double_filter.Last(); ++i)
 168     {
 169         double_filter[i] *= double(1<<(bits+4));
 170         double_filter[i] /= sum;
 171     }
 172
 173     // Turn the float filter into an integer filter
 174     for (int i=double_filter.First(); i<=double_filter.Last(); ++i)
 175     {
 176         int_filter[i] = double_filter[i]>0 ? int( double_filter[i]+0.5 ) : -int( -double_filter[i]+0.5 );
 177         int_filter[i] = (int_filter[i]+8)>>4;
 178     }
 179
 180     return int_filter;
 181 }
 182
 183 void dirac::LPFilter( PicArray& pic_data, const float qf, const int strength )
 184 {
 185     float bw = (std::min( std::max( qf+3.0f-float(strength), 1.0f ), 10.0f ))/10.0;
 186
 187     // filter with 14-bit accuracy
 188     OneDArray<int> filter=MakeLPRectFilter(bw, 14);
 189
 190     HFilter( pic_data, filter, 14 );
 191     VFilter( pic_data, filter, 14 );
 192
 193 }
 194
 195 void HFilter( PicArray& pic_data, const OneDArray<int>& filter, const int bits )
 196 {
 197     ValueType* line_data = new ValueType[pic_data.LengthX()];
 198     const int offset = (1<<(bits-1));
 199
 200     int sum;
 201
 202     for (int j=0; j<pic_data.LengthY(); ++j)
 203     {
 204         // Do the first bit
 205         for (int i=0; i<filter.Last(); ++i)
 206         {
 207             sum = offset;
 208             for (int k=filter.Last(); k>=filter.First(); --k)
 209                 sum += filter[k]*pic_data[j][std::max(i-k,0)];
 210             sum >>= bits;
 211             sum = std::min( 127, std::max( -128, sum) );
 212             line_data[i] = ValueType( sum );
 213         }// i
 214         // Do the middle bit
 215         for (int i=filter.Last(); i<=pic_data.LastX()+filter.First(); ++i)
 216         {
 217             sum = offset;
 218             for (int k=filter.Last(); k>=filter.First(); --k)
 219                 sum += filter[k]*pic_data[j][i-k];
 220             sum >>= bits;
 221             sum = std::min( 127, std::max( -128, sum) );
 222             line_data[i] = ValueType( sum );
 223         }// i
 224         // Do the last bit
 225         for (int i=pic_data.LastX()+filter.First()+1; i<pic_data.LengthX(); ++i)
 226         {
 227             sum = offset;
 228             for (int k=filter.Last(); k>=filter.First(); --k)
 229                 sum += filter[k]*pic_data[j][std::min(i-k,pic_data.LastX())];
 230             sum >>= bits;
 231             sum = std::min( 127, std::max( -128, sum) );
 232             line_data[i] = ValueType( sum );
 233         }// i
 234
 235         // Copy data back
 236
 237         for (int i=0; i<pic_data.LengthX(); ++i )
 238             pic_data[j][i] = line_data[i];
 239
 240     }// j
 241
 242     delete[] line_data;
 243 }
 244
 245 void VFilter( PicArray& pic_data, const OneDArray<int>& filter, const int bits )
 246 {
 247     PicArray tmp_data( pic_data );
 248     const int offset = (1<<(bits-1));
 249
 250     int sum;
 251
 252     // Do the first bit
 253     for (int j=0; j<filter.Last(); ++j)
 254     {
 255
 256         for (int i=0; i<pic_data.LengthX(); ++i)
 257         {
 258             sum = offset;
 259             for (int k=filter.Last(); k>=filter.First(); --k)
 260                 sum += filter[k]*pic_data[std::max(j-k,0)][i];
 261             sum >>= bits;
 262             sum = std::min( 127, std::max( -128, sum) );
 263             tmp_data[j][i] = ValueType( sum );
 264         }// i
 265
 266     }// j
 267
 268     // Do the middle bit
 269     for (int j=filter.Last(); j<=pic_data.LastY()+filter.First(); ++j)
 270     {
 271
 272         for (int i=0; i<pic_data.LengthX(); ++i)
 273         {
 274             sum = offset;
 275             for (int k=filter.Last(); k>=filter.First(); --k)
 276                 sum += filter[k]*pic_data[j-k][i];
 277             sum >>= bits;
 278             sum = std::min( 127, std::max( -128, sum) );
 279             tmp_data[j][i] = ValueType( sum );
 280         }// i
 281
 282     }// j
 283
 284     // Do the last bit
 285     for (int j=pic_data.LastY()+filter.First()+1; j<pic_data.LengthY(); ++j)
 286     {
 287
 288         for (int i=0; i<pic_data.LengthX(); ++i)
 289         {
 290             sum = offset;
 291             for (int k=filter.Last(); k>=filter.First(); --k)
 292                 sum += filter[k]*pic_data[std::min(j-k,pic_data.LastY())][i];
 293             sum >>= bits;
 294             sum = std::min( 127, std::max( -128, sum) );
 295             tmp_data[j][i] = ValueType( sum );
 296         }// i
 297
 298     }// j
 299
 300     // Copy data across
 301     pic_data = tmp_data;
 302
 303 }
 304
 305 /***************************************************************************/
 306
 307 ValueType DiagFilterBchkD( const PicArray& pic,
 308                            const int xpos, const int ypos,
 309                            const TwoDArray<int>& filter,
 310                            const int shift)
 311 {
 312     // Half the filter length
 313     const int len2 = 6;
 314
 315     const int height = pic.LengthY();
 316     const int width = pic.LengthX();
 317
 318     int uplus, uneg, vplus, vneg;
 319     int val = (1<<(shift-1));
 320
 321     // Do 0 position horizontally
 322     val += filter[0][0]*pic[ypos][xpos];
 323
 324     for (int i=1; i<=len2;++i){
 325             uplus = xpos + i;
 326             uplus = (uplus>=width ? width-1 : uplus);
 327             uneg = xpos - i;
 328             uneg = (uneg<0 ? 0 : uneg );
 329             val += filter[0][i]*(pic[ypos][uplus]+pic[ypos][uneg] );
 330     }
 331
 332     // Do other positions vertically//
 333     //////////////////////////////////
 334
 335     for (int j=1; j<=len2;++j){
 336         vplus = ypos + j;
 337         vplus = ( vplus>=height ? height-1 : vplus);
 338         vneg = ypos - j;
 339         vneg = (vneg<0 ? 0 : vneg );
 340
 341         // Do 0 position horizontally
 342         val += filter[j][0]*(pic[vneg][xpos]+pic[vplus][xpos]);
 343         for (int i=1; i<=len2;++i){
 344             uplus = xpos + i;
 345             uplus = (uplus>=width ? width-1 : uplus);
 346             uneg = xpos - i;
 347             uneg = (uneg<0 ? 0 : uneg );
 348             val += filter[j][i]*(pic[vneg][uplus]+pic[vneg][uneg]+
 349                                  pic[vplus][uplus]+pic[vplus][uneg] );
 350         }
 351     }
 352
 353     val >>= shift;
 354
 355     return ValueType(val);
 356 }
 357
 358 ValueType DiagFilterD( const PicArray& pic,
 359                        const int xpos, const int ypos,
 360                        const TwoDArray<int>& filter,
 361                        const int shift)
 362 {
 363     // Half the filter length
 364     const int len2 = 6;
 365
 366     int uplus, uneg, vplus, vneg;
 367     int val = (1<<(shift-1));
 368
 369     // Do 0 position horizontally
 370     val += filter[0][0]*pic[ypos][xpos];
 371
 372     for (int i=1; i<=len2;++i){
 373             uplus = xpos + i;
 374             uneg = xpos - i;
 375             val += filter[0][i]*(pic[ypos][uplus]+pic[ypos][uneg] );
 376     }
 377
 378     // Do other positions vertically//
 379     //////////////////////////////////
 380
 381     for (int j=1; j<=len2;++j){
 382         vplus = ypos + j;
 383         vneg = ypos - j;
 384
 385         // Do 0 position horizontally
 386         val += filter[j][0]*(pic[vneg][xpos]+pic[vplus][xpos]);
 387         for (int i=1; i<=len2;++i){
 388             uplus = xpos + i;
 389             uneg = xpos - i;
 390             val += filter[j][i]*(pic[vneg][uplus]+pic[vneg][uneg]+
 391                                  pic[vplus][uplus]+pic[vplus][uneg] );
 392         }
 393     }
 394
 395     val >>= shift;
 396
 397     return ValueType(val);
 398 }
 399
 400
 401 TwoDArray<int> GetDiagLPFilter( const float bw )
 402 {
 403     TwoDArray<int> f( 7, 7 );
 404
 405     // Bandwidth quantised to range 0.2-1
 406     int qbf = int( bw*10.0 + 0.5 );
 407     qbf = std::min( std::max( qbf, 2 ) , 10 );
 408
 409     switch (qbf){
 410
 411     case 1 :
 412         f[0][0]=1651; f[0][1]=1544;  f[0][2]=1259;  f[0][3]=887;  f[0][4]=530;  f[0][5]=260;  f[0][6]=99;
 413         f[1][0]=1544; f[1][1]=1442;  f[1][2]=1170;  f[1][3]=817;  f[1][4]=480;  f[1][5]=229;  f[1][6]=83;
 414         f[2][0]=1259; f[2][1]=1170;  f[2][2]=935;   f[2][3]=634;  f[2][4]=354;  f[2][5]=153;  f[2][6]=45;
 415         f[3][0]=887;  f[3][1]=817;   f[3][2]=634;   f[3][3]=405;  f[3][4]=202;  f[3][5]=70;   f[3][6]=11;
 416         f[4][0]=530;  f[4][1]=480;   f[4][2]=354;   f[4][3]=202;  f[4][4]=80;   f[4][5]=15;   f[4][6]=0;
 417         f[5][0]=260;  f[5][1]=229;   f[5][2]=153;   f[5][3]=70;   f[5][4]=15;   f[5][5]=0;    f[5][6]=0;
 418         f[6][0]=99;   f[6][1]=83;    f[6][2]=45;    f[6][3]=11;   f[6][4]=0;    f[6][5]=0;    f[6][6]=0;
 419
 420         break;
 421
 422     case 2:
 423
 424         f[0][0]=2855;  f[0][1]=2540;  f[0][2]=1775;  f[0][3]=947;  f[0][4]=364;  f[0][5]=89;  f[0][6]=10;
 425         f[1][0]=2540;  f[1][1]=2251;  f[1][2]=1551;  f[1][3]=804;  f[1][4]=290;  f[1][5]=59;  f[1][6]=1;
 426         f[2][0]=1775;  f[2][1]=1551;  f[2][2]=1020;  f[2][3]=475;  f[2][4]=130;  f[2][5]=3;   f[2][6]=-10;
 427         f[3][0]=947;   f[3][1]=804;   f[3][2]=475;   f[3][3]=165;  f[3][4]=5;    f[3][5]=-22; f[3][6]=-6;
 428         f[4][0]=364;   f[4][1]=290;   f[4][2]=130;   f[4][3]=5;    f[4][4]=-28;  f[4][5]=-10; f[4][6]=0;
 429         f[5][0]=89;    f[5][1]=59;    f[5][2]=3;     f[5][3]=-22;  f[5][4]=-10;  f[5][5]=0;   f[5][6]=0;
 430         f[6][0]=10;    f[6][1]=1;     f[6][2]=-10;   f[6][3]=-6;   f[6][4]=0;    f[6][5]=0;   f[6][6]=0;
 431
 432         break;
 433
 434     case 3:
 435
 436         f[0][0]=5767;  f[0][1]=4718;  f[0][2]=2498;  f[0][3]=745;  f[0][4]=72;   f[0][5]=5;   f[0][6]=23;
 437         f[1][0]=4718;  f[1][1]=3796;  f[1][2]=1875;  f[1][3]=423;  f[1][4]=-58;  f[1][5]=-41; f[1][6]=7;
 438         f[2][0]=2498;  f[2][1]=1875;  f[2][2]=643;   f[2][3]=-146; f[2][4]=-241; f[2][5]=-88; f[2][6]=-9;
 439         f[3][0]=745;   f[3][1]=423;   f[3][2]=-146;  f[3][3]=-367; f[3][4]=-220; f[3][5]=-51; f[3][6]=-2;
 440         f[4][0]=72;    f[4][1]=-58;   f[4][2]=-241;  f[4][3]=-220; f[4][4]=-78;  f[4][5]=-5;  f[4][6]=0;
 441         f[5][0]=5;     f[5][1]=-41;   f[5][2]=-88;   f[5][3]=-51;  f[5][4]=-5;   f[5][5]=0;   f[5][6]=0;
 442         f[6][0]=23;    f[6][1]=7;     f[6][2]=-9;    f[6][3]=-2;   f[6][4]=0;    f[6][5]=0;   f[6][6]=0;
 443
 444         break;
 445
 446     case 4:
 447
 448         f[0][0]=10534;  f[0][1]=7642;  f[0][2]=2603;  f[0][3]=194;  f[0][4]=56;   f[0][5]=120;  f[0][6]=28;
 449         f[1][0]=7642;   f[1][1]=5237;  f[1][2]=1218;  f[1][3]=-383; f[1][4]=-153; f[1][5]=40;   f[1][6]=2;
 450         f[2][0]=2603;   f[2][1]=1218;  f[2][2]=-771;  f[2][3]=-958; f[2][4]=-269; f[2][5]=-3;   f[2][6]=-7;
 451         f[3][0]=194;    f[3][1]=-383;  f[3][2]=-958;  f[3][3]=-541; f[3][4]=-18;  f[3][5]=48;   f[3][6]=4;
 452         f[4][0]=56;     f[4][1]=-153;  f[4][2]=-269;  f[4][3]=-18;  f[4][4]=96;   f[4][5]=22;   f[4][6]=0;
 453         f[5][0]=120;    f[5][1]=40;    f[5][2]=-3;    f[5][3]=48;   f[5][4]=22;   f[5][5]=0;    f[5][6]=0;
 454         f[6][0]=28;     f[6][1]=2;     f[6][2]=-7;    f[6][3]=4;    f[6][4]=0;    f[6][5]=0;    f[6][6]=0;
 455
 456         break;
 457
 458     case 5 :
 459
 460         f[0][0]=16421;  f[0][1]=10159;  f[0][2]=1716;  f[0][3]=33;   f[0][4]=325;  f[0][5]=57;  f[0][6]=6;
 461         f[1][0]=10159;  f[1][1]=5309;   f[1][2]=-580;  f[1][3]=-747; f[1][4]=44;   f[1][5]=-43; f[1][6]=-25;
 462         f[2][0]=1716;   f[2][1]=-580;   f[2][2]=-2310; f[2][3]=-763; f[2][4]=100;  f[2][5]=-19; f[2][6]=-12;
 463         f[3][0]=33;     f[3][1]=-747;   f[3][2]=-763;  f[3][3]=308;  f[3][4]=326;  f[3][5]=27;  f[3][6]=1;
 464         f[4][0]=325;    f[4][1]=44;     f[4][2]=100;   f[4][3]=326;  f[4][4]=84;   f[4][5]=-14; f[4][6]=0;
 465         f[5][0]=57;     f[5][1]=-43;    f[5][2]=-19;   f[5][3]=27;   f[5][4]=-14;  f[5][5]=0;   f[5][6]=0;
 466         f[6][0]=6;      f[6][1]=-25;    f[6][2]=-12;   f[6][3]=1;    f[6][4]=0;    f[6][5]=0;   f[6][6]=0;
 467
 468         break;
 469
 470     case 6 :
 471
 472         f[0][0]=23511;  f[0][1]=11883;  f[0][2]=566;   f[0][3]=524;  f[0][4]=231;  f[0][5]=18;  f[0][6]=41;
 473         f[1][0]= 11883; f[1][1]=3647;   f[1][2]=-2496; f[1][3]=-361; f[1][4]=-96;  f[1][5]=-97; f[1][6]=1;
 474         f[2][0]=566;    f[2][1]=-2496;  f[2][2]=-2329; f[2][3]=459;  f[2][4]=152;  f[2][5]=-7;  f[2][6]=18;
 475         f[3][0]=524;    f[3][1]=-361;   f[3][2]=459;   f[3][3]=979;  f[3][4]=33;   f[3][5]=-28; f[3][6]=3;
 476         f[4][0]=231;    f[4][1]=-96;    f[4][2]=152;   f[4][3]=33;   f[4][4]=-184; f[4][5]=-15; f[4][6]=0;
 477         f[5][0]=18;     f[5][1]=-97;    f[5][2]=-7;    f[5][3]=-28;  f[5][4]=-15;  f[5][5]=0;   f[5][6]=0;
 478         f[6][0]=41;     f[6][1]=1;      f[6][2]=18;    f[6][3]=3;    f[6][4]=0;    f[6][5]=0;   f[6][6]=0;
 479
 480         break;
 481
 482     case 7 :
 483
 484         f[0][0]=32188;  f[0][1]=12652;  f[0][2]=3;      f[0][3]=921;  f[0][4]=1;    f[0][5]=128;  f[0][6]=0;
 485         f[1][0]=12652;  f[1][1]=295;    f[1][2]=-3414;  f[1][3]=-2;   f[1][4]=-343; f[1][5]=-1;   f[1][6]=-37;
 486         f[2][0]=3;      f[2][1]=-3414;  f[2][2]=-212;   f[2][3]=1273; f[2][4]=1;    f[2][5]=98;   f[2][6]=0;
 487         f[3][0]=921;    f[3][1]=-2;     f[3][2]=1273;   f[3][3]=110;  f[3][4]=-363; f[3][5]=0;    f[3][6]=-8;
 488         f[4][0]=1;      f[4][1]=-343;   f[4][2]=1;      f[4][3]=-363; f[4][4]=-29;  f[4][5]=29;   f[4][6]=0;
 489         f[5][0]=128;    f[5][1]=-1;     f[5][2]=98;     f[5][3]=0;    f[5][4]=29;   f[5][5]=0;    f[5][6]=0;
 490         f[6][0]=0;      f[6][1]=-37;    f[6][2]=0;      f[6][3]=-8;   f[6][4]=0;    f[6][5]=0;    f[6][6]=0;
 491
 492         break;
 493
 494     case 8 :
 495
 496         f[0][0]=41902;  f[0][1]=12084;  f[0][2]=435;    f[0][3]=610;   f[0][4]=188;  f[0][5]=34;  f[0][6]=37;
 497         f[1][0]=12084;  f[1][1]=-4268;  f[1][2]=-2715;  f[1][3]=-286;  f[1][4]=-144; f[1][5]=-84; f[1][6]=-2;
 498         f[2][0]=435;    f[2][1]=-2715;  f[2][2]=2809;   f[2][3]=640;   f[2][4]=127;  f[2][5]=10;  f[2][6]=17;
 499         f[3][0]=610;    f[3][1]=-286;   f[3][2]=640;    f[3][3]=-1250; f[3][4]=-45;  f[3][5]=-26; f[3][6]=2;
 500         f[4][0]=188;    f[4][1]=-144;   f[4][2]=127;    f[4][3]=-45;   f[4][4]=259;  f[4][5]=-8;  f[4][6]=0;
 501         f[5][0]=34;     f[5][1]=-84;    f[5][2]=10;     f[5][3]=-26;   f[5][4]=-8;   f[5][5]=0;   f[5][6]=0;
 502         f[6][0]=37;     f[6][1]=-2;     f[6][2]=17;     f[6][3]=2;     f[6][4]=0;    f[6][5]=0;   f[6][6]=0;
 503
 504         break;
 505
 506     case 9 :
 507
 508         f[0][0]=53098;  f[0][1]=10449;  f[0][2]=1546;  f[0][3]=73;    f[0][4]=342;  f[0][5]=38;  f[0][6]=12;
 509         f[1][0]=10449;  f[1][1]=-9060;  f[1][2]=-873;  f[1][3]=-727;  f[1][4]=52;   f[1][5]=-65; f[1][6]=-20;
 510         f[2][0]=1546;   f[2][1]=-873;   f[2][2]=4261;  f[2][3]=-627;  f[2][4]=137;  f[2][5]=-27; f[2][6]=-7;
 511         f[3][0]=73;     f[3][1]=-727;   f[3][2]=-627;  f[3][3]=-804;  f[3][4]=328;  f[3][5]=14;  f[3][6]=2;
 512         f[4][0]=342;    f[4][1]=52;     f[4][2]=137;   f[4][3]=328;   f[4][4]=-83;  f[4][5]=-20; f[4][6]=0;
 513         f[5][0]=38;     f[5][1]=-65;    f[5][2]=-27;   f[5][3]=14;    f[5][4]=-20;  f[5][5]=0;   f[5][6]=0;
 514         f[6][0]=12;     f[6][1]=-20;    f[6][2]=-7;    f[6][3]=2;     f[6][4]=0;    f[6][5]=0;   f[6][6]=0;
 515
 516         break;
 517
 518     default:// case 10
 519
 520         for (int j=0; j<f.LengthY(); ++j ){
 521             for (int i=0; i<f.LengthX(); ++i ){
 522                 f[j][i] = 0;
 523             }
 524         }
 525         f[0][0] = 65536;
 526
 527     }
 528
 529     return f;
 530
 531 }
 532
 533
 534 // Does a diagnonal prefilter
 535 void dirac::DiagFilter( PicArray& pic_data, const float qf, const int strength ){
 536     // One quadrant of the filter taps
 537
 538
 539     float ffactor = (8.0+strength-4.0 - qf )/5.0;
 540     int factor = std::max(0, std::min( 256, int( ffactor*256.0 ) ) ) ;
 541
 542     float bw = (1.0-ffactor)*0.6+0.4;
 543
 544     //std::cout<<std::endl<<"Diagonal prefiltering with bandwidth = "<<bw;
 545
 546     if (bw>0.9)
 547         return;
 548
 549     TwoDArray<int> filter = GetDiagLPFilter( bw );
 550
 551     filter[0][0] = ( factor*filter[0][0] + ( (1<<8)-factor )*(1<<16) + (1<<7) ) >> 8;
 552
 553     for (int i=1;i<7; ++i )
 554         filter[0][i] = ( factor*filter[0][i] + (1<<7) ) >> 8;
 555
 556     for (int j=1;j<7; ++j )
 557         for (int i=0;i<7; ++i )
 558             filter[j][i] = ( factor*filter[j][i] + (1<<7) ) >> 8;
 559
 560
 561     PicArray tmp_data( pic_data.LengthY(), pic_data.LengthX(), pic_data.CSort() );
 562
 563     const int shift = 16;
 564
 565     for (int j=0; j<7;++j)
 566         for (int i=0; i<pic_data.LengthX();++i)
 567             tmp_data[j][i] = DiagFilterBchkD( pic_data, i, j, filter, shift);
 568
 569     for (int j=7; j<pic_data.LengthY()-7;++j){
 570         for (int i=0; i<7;++i)
 571             tmp_data[j][i] = DiagFilterBchkD( pic_data, i, j, filter, shift );
 572
 573         for (int i=7; i<pic_data.LengthX()-7;++i)
 574             tmp_data[j][i] = DiagFilterD( pic_data, i, j, filter, shift );
 575
 576         for (int i=pic_data.LengthX()-7; i<pic_data.LengthX();++i)
 577             tmp_data[j][i] = DiagFilterBchkD( pic_data, i, j, filter, shift );
 578
 579     }
 580
 581     for (int j=pic_data.LengthY()-7; j<pic_data.LengthY();++j)
 582         for (int i=0; i<pic_data.LengthX();++i)
 583             tmp_data[j][i] = DiagFilterBchkD( pic_data, i, j, filter, shift );
 584
 585     pic_data = tmp_data;
 586
 587 }
 588