gimp-montypak/wavelet.c

   1 /*
   2  *   Wavelet Denoise filter for GIMP - The GNU Image Manipulation Program
   3  *
   4  *   Copyright (C) 2008 Monty
   5  *   Code based on research by Crystal Wagner and Prof. Ivan Selesnik,
   6  *   Polytechnic University, Brooklyn, NY
   7  *   See: http://taco.poly.edu/selesi/DoubleSoftware/
   8  *
   9  *   This program is free software; you can redistribute it and/or modify
  10  *   it under the terms of the GNU General Public License as published by
  11  *   the Free Software Foundation; either version 2 of the License, or
  12  *   (at your option) any later version.
  13  *
  14  *   This program is distributed in the hope that it will be useful,
  15  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  *   GNU General Public License for more details.
  18  *
  19  *   You should have received a copy of the GNU General Public License
  20  *   along with this program; if not, write to the Free Software
  21  *   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  22  *
  23  */
  24
  25
  26 #include <string.h>
  27 #include <string.h>
  28 #include <stdlib.h>
  29
  30 #include <libgimp/gimp.h>
  31 #include <libgimp/gimpui.h>
  32
  33 /* convolution code below assumes FSZ is even */
  34 #define FSZ 14
  35
  36 static float analysis_fs[2][3][FSZ]={
  37   {
  38     {
  39       +0.00000000000000,
  40       +0.00069616789827,
  41       -0.02692519074183,
  42       -0.04145457368920,
  43       +0.19056483888763,
  44       +0.58422553883167,
  45       +0.58422553883167,
  46       +0.19056483888763,
  47       -0.04145457368920,
  48       -0.02692519074183,
  49       +0.00069616789827,
  50       +0.00000000000000,
  51       +0.00000000000000,
  52       +0.00000000000000,
  53     },
  54     {
  55       +0.00000000000000,
  56       -0.00014203017443,
  57       +0.00549320005590,
  58       +0.01098019299363,
  59       -0.13644909765612,
  60       -0.21696226276259,
  61       +0.33707999754362,
  62       +0.33707999754362,
  63       -0.21696226276259,
  64       -0.13644909765612,
  65       +0.01098019299363,
  66       +0.00549320005590,
  67       -0.00014203017443,
  68       +0.00000000000000,
  69     },
  70     {
  71       +0.00000000000000,
  72       +0.00014203017443,
  73       -0.00549320005590,
  74       -0.00927404236573,
  75       +0.07046152309968,
  76       +0.13542356651691,
  77       -0.64578354990472,
  78       +0.64578354990472,
  79       -0.13542356651691,
  80       -0.07046152309968,
  81       +0.00927404236573,
  82       +0.00549320005590,
  83       -0.00014203017443,
  84       +0.00000000000000,
  85     }
  86   },
  87   {
  88     {
  89       +0.00000000000000,
  90       +0.00000000000000,
  91       +0.00069616789827,
  92       -0.02692519074183,
  93       -0.04145457368920,
  94       +0.19056483888763,
  95       +0.58422553883167,
  96       +0.58422553883167,
  97       +0.19056483888763,
  98       -0.04145457368920,
  99       -0.02692519074183,
 100       +0.00069616789827,
 101       +0.00000000000000,
 102       +0.00000000000000,
 103     },
 104     {
 105       +0.00000000000000,
 106       +0.00000000000000,
 107       -0.00014203017443,
 108       +0.00549320005590,
 109       +0.01098019299363,
 110       -0.13644909765612,
 111       -0.21696226276259,
 112       +0.33707999754362,
 113       +0.33707999754362,
 114       -0.21696226276259,
 115       -0.13644909765612,
 116       +0.01098019299363,
 117       +0.00549320005590,
 118       -0.00014203017443,
 119     },
 120     {
 121       +0.00000000000000,
 122       +0.00000000000000,
 123       +0.00014203017443,
 124       -0.00549320005590,
 125       -0.00927404236573,
 126       +0.07046152309968,
 127       +0.13542356651691,
 128       -0.64578354990472,
 129       +0.64578354990472,
 130       -0.13542356651691,
 131       -0.07046152309968,
 132       +0.00927404236573,
 133       +0.00549320005590,
 134       -0.00014203017443,
 135     }
 136   }
 137 };
 138
 139 static float analysis[2][3][FSZ]={
 140   {
 141     {
 142       +0.00000000000000,
 143       +0.00000000000000,
 144       +0.00017870679071,
 145       -0.02488304194507,
 146       +0.00737700819766,
 147       +0.29533805776119,
 148       +0.59529279993637,
 149       +0.45630440337458,
 150       +0.11239376309619,
 151       -0.01971220693439,
 152       -0.00813549683439,
 153       +0.00005956893024,
 154       +0.00000000000000,
 155       +0.00000000000000,
 156     },
 157     {
 158       +0.00000000000000,
 159       +0.00000000000000,
 160       -0.00012344587034,
 161       +0.01718853971559,
 162       -0.00675291099550,
 163       +0.02671809818132,
 164       -0.64763513288874,
 165       +0.47089724990858,
 166       +0.16040017815754,
 167       -0.01484700537727,
 168       -0.00588868840296,
 169       +0.00004311757177,
 170       +0.00000000000000,
 171       +0.00000000000000,
 172     },
 173     {
 174       +0.00000000000000,
 175       +0.00000000000000,
 176       +0.00001437252392,
 177       -0.00200122286479,
 178       +0.00027261232228,
 179       +0.06840460220387,
 180       +0.01936710587994,
 181       -0.68031992557818,
 182       +0.42976785708978,
 183       +0.11428688385011,
 184       +0.05057805218407,
 185       -0.00037033761102,
 186       +0.00000000000000,
 187       +0.00000000000000,
 188     }
 189   },
 190   {
 191     {
 192       +0.00000000000000,
 193       +0.00000000000000,
 194       +0.00005956893024,
 195       -0.00813549683439,
 196       -0.01971220693439,
 197       +0.11239376309619,
 198       +0.45630440337458,
 199       +0.59529279993637,
 200       +0.29533805776119,
 201       +0.00737700819766,
 202       -0.02488304194507,
 203       +0.00017870679071,
 204       +0.00000000000000,
 205       +0.00000000000000,
 206     },
 207
 208     {
 209       +0.00000000000000,
 210       +0.00000000000000,
 211       -0.00037033761102,
 212       +0.05057805218407,
 213       +0.11428688385011,
 214       +0.42976785708978,
 215       -0.68031992557818,
 216       +0.01936710587994,
 217       +0.06840460220387,
 218       +0.00027261232228,
 219       -0.00200122286479,
 220       +0.00001437252392,
 221       +0.00000000000000,
 222       +0.00000000000000,
 223     },
 224     {
 225       +0.00000000000000,
 226       +0.00000000000000,
 227       +0.00004311757177,
 228       -0.00588868840296,
 229       -0.01484700537727,
 230       +0.16040017815754,
 231       +0.47089724990858,
 232       -0.64763513288874,
 233       +0.02671809818132,
 234       -0.00675291099550,
 235       +0.01718853971559,
 236       -0.00012344587034,
 237       +0.00000000000000,
 238       +0.00000000000000,
 239     }
 240   }
 241 };
 242
 243 static float synthesis_fs[2][3][FSZ]={
 244   {
 245     {
 246       +0.00000000000000,
 247       +0.00000000000000,
 248       +0.00000000000000,
 249       +0.00069616789827,
 250       -0.02692519074183,
 251       -0.04145457368920,
 252       +0.19056483888763,
 253       +0.58422553883167,
 254       +0.58422553883167,
 255       +0.19056483888763,
 256       -0.04145457368920,
 257       -0.02692519074183,
 258       +0.00069616789827,
 259       +0.00000000000000,
 260     },
 261     {
 262       +0.00000000000000,
 263       -0.00014203017443,
 264       +0.00549320005590,
 265       +0.01098019299363,
 266       -0.13644909765612,
 267       -0.21696226276259,
 268       +0.33707999754362,
 269       +0.33707999754362,
 270       -0.21696226276259,
 271       -0.13644909765612,
 272       +0.01098019299363,
 273       +0.00549320005590,
 274       -0.00014203017443,
 275       +0.00000000000000,
 276     },
 277     {
 278       +0.00000000000000,
 279       -0.00014203017443,
 280       +0.00549320005590,
 281       +0.00927404236573,
 282       -0.07046152309968,
 283       -0.13542356651691,
 284       +0.64578354990472,
 285       -0.64578354990472,
 286       +0.13542356651691,
 287       +0.07046152309968,
 288       -0.00927404236573,
 289       -0.00549320005590,
 290       +0.00014203017443,
 291       +0.00000000000000,
 292     }
 293   },
 294   {
 295     {
 296       +0.00000000000000,
 297       +0.00000000000000,
 298       +0.00069616789827,
 299       -0.02692519074183,
 300       -0.04145457368920,
 301       +0.19056483888763,
 302       +0.58422553883167,
 303       +0.58422553883167,
 304       +0.19056483888763,
 305       -0.04145457368920,
 306       -0.02692519074183,
 307       +0.00069616789827,
 308       +0.00000000000000,
 309       +0.00000000000000,
 310     },
 311     {
 312       -0.00014203017443,
 313       +0.00549320005590,
 314       +0.01098019299363,
 315       -0.13644909765612,
 316       -0.21696226276259,
 317       +0.33707999754362,
 318       +0.33707999754362,
 319       -0.21696226276259,
 320       -0.13644909765612,
 321       +0.01098019299363,
 322       +0.00549320005590,
 323       -0.00014203017443,
 324       +0.00000000000000,
 325       +0.00000000000000,
 326     },
 327     {
 328       -0.00014203017443,
 329       +0.00549320005590,
 330       +0.00927404236573,
 331       -0.07046152309968,
 332       -0.13542356651691,
 333       +0.64578354990472,
 334       -0.64578354990472,
 335       +0.13542356651691,
 336       +0.07046152309968,
 337       -0.00927404236573,
 338       -0.00549320005590,
 339       +0.00014203017443,
 340       +0.00000000000000,
 341       +0.00000000000000,
 342     }
 343   }
 344 };
 345
 346 static float synthesis[2][3][FSZ]={
 347   {
 348     {
 349       +0.00000000000000,
 350       +0.00000000000000,
 351       +0.00005956893024,
 352       -0.00813549683439,
 353       -0.01971220693439,
 354       +0.11239376309619,
 355       +0.45630440337458,
 356       +0.59529279993637,
 357       +0.29533805776119,
 358       +0.00737700819766,
 359       -0.02488304194507,
 360       +0.00017870679071,
 361       +0.00000000000000,
 362       +0.00000000000000,
 363     },
 364     {
 365       +0.00000000000000,
 366       +0.00000000000000,
 367       +0.00004311757177,
 368       -0.00588868840296,
 369       -0.01484700537727,
 370       +0.16040017815754,
 371       +0.47089724990858,
 372       -0.64763513288874,
 373       +0.02671809818132,
 374       -0.00675291099550,
 375       +0.01718853971559,
 376       -0.00012344587034,
 377       +0.00000000000000,
 378       +0.00000000000000,
 379     },
 380     {
 381       +0.00000000000000,
 382       +0.00000000000000,
 383       -0.00037033761102,
 384       +0.05057805218407,
 385       +0.11428688385011,
 386       +0.42976785708978,
 387       -0.68031992557818,
 388       +0.01936710587994,
 389       +0.06840460220387,
 390       +0.00027261232228,
 391       -0.00200122286479,
 392       +0.00001437252392,
 393       +0.00000000000000,
 394       +0.00000000000000,
 395     }
 396   },
 397   {
 398     {
 399       +0.00000000000000,
 400       +0.00000000000000,
 401       +0.00017870679071,
 402       -0.02488304194507,
 403       +0.00737700819766,
 404       +0.29533805776119,
 405       +0.59529279993637,
 406       +0.45630440337458,
 407       +0.11239376309619,
 408       -0.01971220693439,
 409       -0.00813549683439,
 410       +0.00005956893024,
 411       +0.00000000000000,
 412       +0.00000000000000,
 413     },
 414     {
 415       +0.00000000000000,
 416       +0.00000000000000,
 417       +0.00001437252392,
 418       -0.00200122286479,
 419       +0.00027261232228,
 420       +0.06840460220387,
 421       +0.01936710587994,
 422       -0.68031992557818,
 423       +0.42976785708978,
 424       +0.11428688385011,
 425       +0.05057805218407,
 426       -0.00037033761102,
 427       +0.00000000000000,
 428       +0.00000000000000,
 429     },
 430     {
 431       +0.00000000000000,
 432       +0.00000000000000,
 433       -0.00012344587034,
 434       +0.01718853971559,
 435       -0.00675291099550,
 436       +0.02671809818132,
 437       -0.64763513288874,
 438       +0.47089724990858,
 439       +0.16040017815754,
 440       -0.01484700537727,
 441       -0.00588868840296,
 442       +0.00004311757177,
 443       +0.00000000000000,
 444       +0.00000000000000,
 445     }
 446   }
 447 };
 448
 449 typedef struct {
 450   float *x;
 451   int rows;
 452   int cols;
 453 } m2D;
 454
 455 typedef struct {
 456   m2D w[8][2][2];
 457 } wrank;
 458
 459 static m2D alloc_m2D(int m, int n){
 460   m2D ret;
 461   ret.rows = m;
 462   ret.cols = n;
 463   ret.x = calloc(m*n,sizeof(*(ret.x)));
 464   return ret;
 465 }
 466
 467 static void free_m2D(m2D *c){
 468   if(c->x) free(c->x);
 469   c->x = NULL;
 470 }
 471
 472 static void complex_threshold(m2D set[4], float T, int soft, int pt, int *pc, int pi, int(*check)(void)){
 473   int i,j;
 474   int N = set[0].rows*set[0].cols;
 475
 476   if(T>.01){
 477     float TT = T*T;
 478
 479     if(soft){
 480       for(i=0;i<N;){
 481         for(j=0;j<set[0].cols;i++,j++){
 482           float v00 = (set[0].x[i] + set[3].x[i]) * .7071067812;
 483           float v11 = (set[0].x[i] - set[3].x[i]) * .7071067812;
 484           float v01 = (set[1].x[i] + set[2].x[i]) * .7071067812;
 485           float v10 = (set[1].x[i] - set[2].x[i]) * .7071067812;
 486
 487           if(v00*v00 + v10*v10 < TT){
 488             v00 = 0.;
 489             v10 = 0.;
 490           }else{
 491             float y = sqrt(v00*v00 + v10*v10);
 492             float scale = y/(y+T);
 493             v00 *= scale;
 494             v10 *= scale;
 495           }
 496           if(v01*v01 + v11*v11 < TT){
 497             v01 = 0.;
 498             v11 = 0.;
 499           }else{
 500             float y = sqrt(v01*v01 + v11*v11);
 501             float scale = y/(y+T);
 502             v01 *= scale;
 503             v11 *= scale;
 504           }
 505
 506           set[0].x[i] = (v00 + v11) * .7071067812;
 507           set[3].x[i] = (v00 - v11) * .7071067812;
 508           set[1].x[i] = (v01 + v10) * .7071067812;
 509           set[2].x[i] = (v01 - v10) * .7071067812;
 510         }
 511         if(check && check())return;
 512       }
 513     }else{
 514       for(i=0;i<N;){
 515         for(j=0;j<set[0].cols;i++,j++){
 516           float v00 = (set[0].x[i] + set[3].x[i]) * .7071067812;
 517           float v11 = (set[0].x[i] - set[3].x[i]) * .7071067812;
 518           float v01 = (set[1].x[i] + set[2].x[i]) * .7071067812;
 519           float v10 = (set[1].x[i] - set[2].x[i]) * .7071067812;
 520
 521           if(v00*v00 + v10*v10 < TT){
 522             v00 = 0.;
 523             v10 = 0.;
 524           }
 525           if(v01*v01 + v11*v11 < TT){
 526             v01 = 0.;
 527             v11 = 0.;
 528           }
 529
 530           set[0].x[i] = (v00 + v11) * .7071067812;
 531           set[3].x[i] = (v00 - v11) * .7071067812;
 532           set[1].x[i] = (v01 + v10) * .7071067812;
 533           set[2].x[i] = (v01 - v10) * .7071067812;
 534         }
 535         if(check && check())return;
 536       }
 537     }
 538   }
 539   if(pt)gimp_progress_update((gfloat)((*pc)+=pi)/pt);
 540 }
 541
 542 /* assume the input is padded, return output that's padded for next stage */
 543 /* FSZ*2-2 padding, +1 additional padding if vector odd */
 544
 545 static m2D convolve_padded(const m2D x, float *f, int pt, int *pc, int pi, int(*check)(void)){
 546   int i,M = x.rows;
 547   int j,in_N = x.cols;
 548   int k;
 549   int comp_N = (in_N - FSZ + 3) / 2;
 550   int out_N = (comp_N+1)/2*2+FSZ*2-2;
 551
 552   m2D y=alloc_m2D(M,out_N);
 553   if(check && check())return y;
 554
 555   for(i=0;i<M;i++){
 556     float *xrow = x.x+i*in_N;
 557     float *yrow = y.x+i*out_N+FSZ-1;
 558     for(j=0;j<comp_N;j++){
 559       float a = 0;
 560       for(k=0;k<FSZ;k++)
 561         a+=xrow[k]*f[FSZ-k-1];
 562       xrow+=2;
 563       yrow[j]=a;
 564     }
 565     if(check && check())return y;
 566
 567     /* extend output padding */
 568     for(j=0;j<FSZ-1;j++){
 569       yrow--;
 570       yrow[0]=yrow[1];
 571     }
 572     for(j=FSZ-1+comp_N;j<out_N;j++)
 573       yrow[j] = yrow[j-1];
 574
 575     if(check && check())return y;
 576
 577   }
 578   if(pt)gimp_progress_update((gfloat)((*pc)+=pi)/pt);
 579   return y;
 580 }
 581
 582 static m2D convolve_transpose_padded(const m2D x, float *f, int pt, int *pc, int pi, int(*check)(void)){
 583   int i,M = x.rows;
 584   int j,in_N = x.cols;
 585   int k;
 586   int comp_N = (in_N - FSZ + 3) / 2;
 587   int out_N = (comp_N+1)/2*2+FSZ*2-2;
 588
 589   m2D y=alloc_m2D(out_N,M);
 590   if(check && check())return y;
 591
 592   for(i=0;i<M;i++){
 593     float *xrow = x.x+i*in_N;
 594     float *ycol = y.x + i + M*(FSZ-1);
 595     for(j=0;j<comp_N;j++){
 596       float a = 0;
 597       for(k=0;k<FSZ;k++)
 598         a+=xrow[k]*f[FSZ-k-1];
 599       xrow+=2;
 600       ycol[j*M]=a;
 601     }
 602     if(check && check())return y;
 603
 604     /* extend output padding */
 605     for(j=0;j<FSZ-1;j++){
 606       ycol -= M;
 607       ycol[0] = ycol[M];
 608     }
 609     for(j=FSZ-1+comp_N;j<out_N;j++)
 610       ycol[j*M] = ycol[(j-1)*M];
 611
 612     if(check && check())return y;
 613
 614   }
 615   if(pt)gimp_progress_update((gfloat)((*pc)+=pi)/pt);
 616   return y;
 617 }
 618
 619 /* discards the padding added by the matching convolution */
 620
 621 /* y will often be smaller than a full x expansion due to preceeding
 622    rounds of padding out to even values; this padding is also
 623    discarded */
 624 static void deconvolve_padded(m2D y, m2D x, float *f, int pt, int *pc, int pi, int(*check)(void)){
 625   int i;
 626   int j,in_N = x.cols;
 627   int k;
 628   int out_N = y.cols;
 629   int out_M = y.rows;
 630
 631   for(i=0;i<out_M;i++){
 632     float *xrow = x.x+i*in_N+FSZ/2;
 633     float *yrow = y.x+i*out_N;
 634
 635     for(j=0;j<out_N;j+=2){
 636       float a = 0;
 637       for(k=1;k<FSZ;k+=2)
 638         a+=xrow[k>>1]*f[FSZ-k-1];
 639       yrow[j]+=a;
 640       a=0.;
 641       for(k=1;k<FSZ;k+=2)
 642         a+=xrow[k>>1]*f[FSZ-k];
 643       yrow[j+1]+=a;
 644       xrow++;
 645     }
 646     if(check && check()) return;
 647   }
 648
 649   if(pt)gimp_progress_update((gfloat)((*pc)+=pi)/pt);
 650 }
 651
 652 /* discards the padding added by the matching convolution */
 653
 654 /* y will often be smaller than a full x expansion due to preceeding
 655    rounds of padding out to even values; this padding is also
 656    discarded */
 657 static void deconvolve_transpose_padded(m2D y, m2D x, float *f, int pt, int *pc, int pi, int(*check)(void)){
 658   int i;
 659   int j,in_N = x.cols;
 660   int k;
 661   int out_M = y.rows;
 662   int out_N = y.cols;
 663
 664   for(i=0;i<out_N;i++){
 665     float *xrow = x.x+i*in_N+FSZ/2;
 666     float *ycol = y.x+i;
 667     for(j=0;j<out_M;j+=2){
 668       float a = 0;
 669       for(k=1;k<FSZ;k+=2)
 670         a+=xrow[k>>1]*f[FSZ-k-1];
 671       ycol[j*out_N]+=a;
 672       a=0.;
 673       for(k=1;k<FSZ;k+=2)
 674         a+=xrow[k>>1]*f[FSZ-k];
 675       ycol[(j+1)*out_N]+=a;
 676       xrow++;
 677     }
 678     if(check && check()) return;
 679   }
 680
 681   if(pt)gimp_progress_update((gfloat)((*pc)+=pi)/pt);
 682 }
 683
 684 /* consumes and replaces lo if free_input set, otherwise overwrites */
 685 static void forward_threshold(m2D lo[4], m2D y[4],
 686                               float af[2][3][FSZ], float sf[2][3][FSZ],
 687                               float T, int soft,
 688                               int free_input, int pt, int *pc, int pi,
 689                               int(*check)(void)){
 690   m2D x[4] = {{NULL,0,0},{NULL,0,0},{NULL,0,0},{NULL,0,0}};
 691   m2D temp[4] = {{NULL,0,0},{NULL,0,0},{NULL,0,0},{NULL,0,0}};
 692   m2D tempw[4] = {{NULL,0,0},{NULL,0,0},{NULL,0,0},{NULL,0,0}};
 693   m2D temp2[4] = {{NULL,0,0},{NULL,0,0},{NULL,0,0},{NULL,0,0}};
 694   int r,c,i;
 695
 696   for(i=0;i<4;i++){
 697     x[i] = lo[i];
 698     lo[i] = (m2D){NULL,0,0};
 699   }
 700
 701   for(i=0;i<4;i++){
 702     temp[i]  = convolve_transpose_padded(x[i], af[i>>1][2], pt, pc, pi, check);
 703     if(check && check())goto cleanup;
 704     tempw[i] = convolve_padded(temp[i], af[i&1][2], pt, pc, pi, check); /* w 7 */
 705     if(check && check())goto cleanup;
 706   }
 707
 708   r = tempw[0].rows;
 709   c = tempw[0].cols;
 710   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 711   if(check && check())goto cleanup;
 712
 713   for(i=0;i<4;i++){
 714     temp2[i]=alloc_m2D(c*2 - FSZ*3 + 2, r);
 715     y[i] = alloc_m2D(temp2[i].rows, temp2[i].cols*2 - FSZ*3 + 2);
 716     if(check && check())goto cleanup;
 717     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][2], pt, pc, pi, check);
 718     if(check && check())goto cleanup;
 719     free_m2D(tempw+i);
 720     tempw[i] = convolve_padded(temp[i], af[i&1][1], pt, pc, pi, check); /* w6 */
 721     if(check && check())goto cleanup;
 722   }
 723   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 724   if(check && check())goto cleanup;
 725
 726   for(i=0;i<4;i++){
 727     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][1], pt, pc, pi, check);
 728     if(check && check())goto cleanup;
 729     free_m2D(tempw+i);
 730     tempw[i] = convolve_padded(temp[i], af[i&1][0], pt, pc, pi, check); /* w5 */
 731     if(check && check())goto cleanup;
 732     free_m2D(temp+i);
 733   }
 734   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 735   if(check && check())goto cleanup;
 736
 737   for(i=0;i<4;i++){
 738     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][0], pt, pc, pi, check);
 739     if(check && check())goto cleanup;
 740     free_m2D(tempw+i);
 741     deconvolve_padded(y[i],temp2[i],sf[i>>1][2], pt, pc, pi, check);
 742     if(check && check())goto cleanup;
 743     free_m2D(temp2+i);
 744     temp[i]   = convolve_transpose_padded(x[i], af[i>>1][1], pt, pc, pi, check);
 745     if(check && check())goto cleanup;
 746     temp2[i]  = alloc_m2D(c*2 - FSZ*3 + 2, r);
 747     if(check && check())goto cleanup;
 748     tempw[i]  = convolve_padded(temp[i], af[i&1][2], pt, pc, pi, check); /* w4 */
 749     if(check && check())goto cleanup;
 750   }
 751   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 752   if(check && check())goto cleanup;
 753
 754   for(i=0;i<4;i++){
 755     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][2], pt, pc, pi, check);
 756     if(check && check())goto cleanup;
 757     free_m2D(tempw+i);
 758     tempw[i] = convolve_padded(temp[i], af[i&1][1], pt, pc, pi, check); /* w3 */
 759     if(check && check())goto cleanup;
 760   }
 761   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 762   if(check && check())goto cleanup;
 763
 764   for(i=0;i<4;i++){
 765     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][1], pt, pc, pi, check);
 766     if(check && check())goto cleanup;
 767     free_m2D(tempw+i);
 768     tempw[i] = convolve_padded(temp[i], af[i&1][0], pt, pc, pi, check); /* w2 */
 769     if(check && check())goto cleanup;
 770     free_m2D(temp+i);
 771   }
 772   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 773   if(check && check())goto cleanup;
 774
 775   for(i=0;i<4;i++){
 776     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][0], pt, pc, pi, check);
 777     if(check && check())goto cleanup;
 778     free_m2D(tempw+i);
 779     deconvolve_padded(y[i],temp2[i],sf[i>>1][1], pt, pc, pi, check);
 780     if(check && check())goto cleanup;
 781     free_m2D(temp2+i);
 782     temp[i]  = convolve_transpose_padded(x[i], af[i>>1][0], pt, pc, pi, check);
 783     if(check && check())goto cleanup;
 784     if(free_input) free_m2D(x+i);
 785     temp2[i] = alloc_m2D(c*2 - FSZ*3 + 2, r);
 786     if(check && check())goto cleanup;
 787     tempw[i] = convolve_padded(temp[i], af[i&1][2], pt, pc, pi, check); /* w1 */
 788     if(check && check())goto cleanup;
 789   }
 790   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 791   if(check && check())goto cleanup;
 792
 793   for(i=0;i<4;i++){
 794     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][2], pt, pc, pi, check);
 795     if(check && check())goto cleanup;
 796     free_m2D(tempw+i);
 797     tempw[i] = convolve_padded(temp[i], af[i&1][1], pt, pc, pi, check); /* w0 */
 798     if(check && check())goto cleanup;
 799   }
 800   complex_threshold(tempw, T, soft, pt, pc, pi, check);
 801   if(check && check())goto cleanup;
 802
 803   for(i=0;i<4;i++){
 804     deconvolve_transpose_padded(temp2[i],tempw[i],sf[i&1][1], pt, pc, pi, check);
 805     if(check && check())goto cleanup;
 806     deconvolve_padded(y[i],temp2[i],sf[i>>1][0], pt, pc, pi, check);
 807     if(check && check())goto cleanup;
 808     lo[i] = convolve_padded(temp[i], af[i&1][0], pt, pc, pi, check); /* lo */
 809     if(check && check())goto cleanup;
 810     free_m2D(temp+i);
 811   }
 812
 813  cleanup:
 814   for(i=0;i<4;i++){
 815     if(free_input)free_m2D(x+i);
 816     free_m2D(temp+i);
 817     free_m2D(tempw+i);
 818     free_m2D(temp2+i);
 819   }
 820 }
 821
 822 /* consumes/replaces lo */
 823 static void finish_backward(m2D lo[4], m2D y[4], float sf[2][3][FSZ], int pt, int *pc, int pi, int(*check)(void)){
 824   int r=lo[0].rows,c=lo[0].cols,i;
 825   m2D temp = {NULL,0,0};
 826
 827   for(i=0;i<4;i++){
 828     temp=alloc_m2D(c*2 - FSZ*3 + 2, r);
 829     deconvolve_transpose_padded(temp,lo[i],sf[i&1][0], pt, pc, pi, check);
 830     free_m2D(lo+i);
 831     if(check && check()){
 832       free_m2D(&temp);
 833       return;
 834     }
 835     deconvolve_padded(y[i],temp,sf[i>>1][0], pt, pc, pi, check);
 836     free_m2D(&temp);
 837     lo[i]=y[i];
 838     y[i]=(m2D){NULL,0,0};
 839     if(check && check()) return;
 840   }
 841 }
 842
 843 static m2D transform_threshold(m2D x, int J, float T[16], int soft, int pt, int *pc, int (*check)(void)){
 844   int i,j;
 845   m2D partial_y[J][4];
 846   m2D lo[4];
 847
 848   for(j=0;j<J;j++)
 849     for(i=0;i<4;i++)
 850       partial_y[j][i]=(m2D){NULL,0,0};
 851
 852   lo[0] = lo[1] = lo[2] = lo[3] = x;
 853
 854   forward_threshold(lo, partial_y[0], analysis_fs, synthesis_fs, T[0], soft, 0, pt, pc, (1<<18), check);
 855   if(check && check()) goto cleanup3;
 856
 857   for(j=1;j<J;j++){
 858     forward_threshold(lo, partial_y[j], analysis, synthesis, T[j], soft, 1, pt, pc, (1<<(9-j>0?(9-j)*2:1)), check);
 859     if(check && check()) goto cleanup3;
 860   }
 861
 862   for(j=J-1;j;j--){
 863     finish_backward(lo, partial_y[j], synthesis, pt, pc, (1<<(9-j>0?(9-j)*2:1)), check);
 864     if(check && check()) goto cleanup3;
 865   }
 866
 867   finish_backward(lo, partial_y[0], synthesis_fs, pt, pc, (1<<18), check);
 868   if(check && check()) goto cleanup3;
 869
 870   {
 871     int end = lo[0].rows*lo[0].cols;
 872     float *y = lo[0].x;
 873     float *p1 = lo[1].x;
 874     float *p2 = lo[2].x;
 875     float *p3 = lo[3].x;
 876
 877     for(j=0;j<end;j++)
 878       y[j]+=p1[j]+p2[j]+p3[j];
 879
 880     if(check && check()) goto cleanup3;
 881   }
 882
 883  cleanup3:
 884
 885   for(j=0;j<J;j++)
 886     for(i=0;i<4;i++)
 887       free_m2D(&partial_y[j][i]);
 888
 889   free_m2D(lo+1);
 890   free_m2D(lo+2);
 891   free_m2D(lo+3);
 892   return lo[0];
 893
 894 }
 895
 896 static int wavelet_filter(int width, int height, float *buffer,
 897                             int pr, float T[16], int soft, int (*check)(void)){
 898
 899   int J=4;
 900   int i,j;
 901   m2D xc={NULL,0,0};
 902   m2D yc={NULL,0,0};
 903   int pc=0;
 904   int pt=0;
 905   float *ptr = buffer;
 906
 907   if(check && check())goto abort;
 908
 909   /* we want J to be as 'deep' as the image to eliminate
 910      splotchiness with deep coarse settings */
 911
 912   while(1){
 913     int mult = 1 << (J+1);
 914     if(width/mult < 1) break;
 915     if(height/mult < 1) break;
 916     J++;
 917   }
 918
 919   if(J>15)J=15;
 920   if(pr)
 921     for(i=0;i<J;i++)
 922       pt+=108*(1<<(9-i>0?(9-i)*2:1));
 923
 924   /* Input matrix must be pre-padded for first stage convolutions;
 925      odd->even padding goes on the bottom/right */
 926
 927   xc = alloc_m2D((height+1)/2*2+FSZ*2-2,
 928                  (width+1)/2*2+FSZ*2-2),yc;
 929   if(check && check())goto abort;
 930
 931   /* populate and pad input matrix */
 932   for(i=0;i<height;i++){
 933     float *row=xc.x + (i+FSZ-1)*xc.cols;
 934
 935     /* X padding */
 936     for(j=0;j<FSZ-1;j++)
 937       row[j] = *ptr * .5f;
 938
 939     /* X filling */
 940     for(;j<width+FSZ-1;j++){
 941       row[j] = *ptr * .5f;
 942       ptr++;
 943     }
 944
 945     /* X padding */
 946     for(;j<xc.cols;j++)
 947       row[j] = row[j-1];
 948   }
 949
 950   /* Y padding */
 951   for(i=FSZ-2;i>=0;i--){
 952     float *pre=xc.x + (i+1)*xc.cols;
 953     float *row=xc.x + i*xc.cols;
 954     for(j=0;j<xc.cols;j++)
 955       row[j]=pre[j];
 956   }
 957   for(i=xc.rows-FSZ+1;i<xc.rows;i++){
 958     float *pre=xc.x + (i-1)*xc.cols;
 959     float *row=xc.x + i*xc.cols;
 960     for(j=0;j<xc.cols;j++)
 961       row[j]=pre[j];
 962   }
 963
 964   if(check && check())goto abort;
 965   yc=transform_threshold(xc,J,T,soft,pt,&pc,check);
 966   if(check && check())goto abort;
 967
 968   /* pull filtered values back out of padded matrix */
 969   ptr = buffer;
 970   for(i=0;i<height;i++){
 971     float *row = yc.x + (i+FSZ-1)*yc.cols + FSZ-1;
 972
 973     for(j=0;j<width;j++)
 974       *ptr++ = row[j]*.5f;
 975   }
 976
 977  abort:
 978   free_m2D(&yc);
 979   free_m2D(&xc);
 980   return (check && check());
 981 }
 982