recover.c File Reference

00071 {
00072         int i,j,limi=pi+di,limj=pj+dj;
00073 
00074         for(i=pi;i<limi;i++) {
00075 
00076                 for(j=pj;j<limj;j++) {
00077 
00078                         if(im[i][j]>max)
00079                                 im[i][j]=max;
00080                         else if(im[i][j]<min)
00081                                 im[i][j]=min;
00082                 }
00083         }
00084 }

00124 {
00125         int layer,i,j;
00126         float my_mse=0,my_thres;
00127 
00128         // must be even for the implemented layering.
00129         const int di=MISSING_BL_SIZE,dj=MISSING_BL_SIZE,num_layers=NUM_LAYERS;
00130 
00131         float noise_thres=5;
00132         const float dt=(float).1;
00133         int main_iter,c,C=INNER_ITER_COUNT;
00134         float min=0,max=255;
00135 
00136         // get the max and min pixel value from the image (ignore missing regions
00137         // in min, max computation) to help clipping below.
00138         // This is unnecessarily paranoid, you can comment this out, i.e., set min=0,max=255,
00139         // and the results shouldn't change much.
00140         find_max_min_avoid(recovered,Ni,Nj,pi,pj,di,dj,&max,&min);
00141 
00142         // if threshold is below the noise threshold then no iterations will be done unless
00143         // MIN_ITER_COUNT is set. If set we will force iterations when needed.
00144         if(MIN_ITER_COUNT>(int)((1.0/dt)*(threshold-noise_thres))) {
00145 
00146                 threshold=noise_thres+MIN_ITER_COUNT*dt;
00147 
00148                 printf("Threshold too low, setting iteration count to MIN_ITER_COUNT\n");
00149         }
00150         
00151         my_thres=threshold;
00152 
00153         main_iter=(int)((1.0/dt)*(threshold-noise_thres)+.5);
00154         main_iter=(main_iter<0)?0:main_iter;
00155 
00156         if(main_iter>MAX_ITER_COUNT) {
00157 
00158                 main_iter=MAX_ITER_COUNT;
00159                 printf("\nMAX_ITER_COUNT hit, will do reduced number of iterations. Please adjust\n");
00160                 printf("parameters if you want more iterations\n\n");
00161         }
00162 
00163         for(j=0;j<main_iter;j++) { // progressive thresholds loop.
00164 
00165                 for(layer=0;layer<num_layers;layer++) { //layer loop
00166 
00167                         // recover layer, keeping everything else unchanged.
00168                         if(C==1) {
00169 
00170                                 // This is more efficient if C==1.
00171                                 // For most cases C=1 provides good solutions.
00172                                 process_layer(recovered,Ni,Nj,recovered,pi+layer,pj+layer,di-2*layer,dj-2*layer,my_thres);
00173                         }
00174                         else {
00175 
00176                                 for(c=0;c<C;c++) {
00177 
00178                                         // This one does C iterations with a given insignificant set.
00179                                         // If C is set to high value this becomes roughly equivalent to
00180                                         // solving the linear equations shown in the paper.
00181                                         // 
00182                                         // Solving the equations tends to take the sparsity constraints
00183                                         // (established through thresholding) more seriously than doing a single iteration.
00184                                         process_layer_w_signif_sets(recovered,Ni,Nj,recovered,pi+layer,pj+layer,di-2*layer,dj-2*layer,my_thres);
00185                                 }
00186                         }
00187 
00188                         // clip layer pixels to min, max calculated above.
00189                         clip_to_min_max_region(recovered,pi+layer,pj+layer,di-2*layer,dj-2*layer,min,max);
00190                 }
00191 
00192 
00193                 // print out progress.
00194                 if(PRINT_ITERS) {
00195 
00196                         my_mse=mse_region(original,recovered,pi,pj,di,dj);
00197 
00198                         printf("iteration=%3d, C=%3d, thres=%8.3f, PSNR %8.3f dB\n",j,C,my_thres,psnr(my_mse));
00199                 }
00200 
00201                 // Reduce threshold.
00202                 my_thres-=dt;
00203         }
00204 
00205         printf("\nPerformed %5d iterations\n\n",j);
00206 
00207         // recovered is clipped already, so it should be 0-255.
00208         // just round to integers.
00209         for(i=pi;i<pi+di;i++) {
00210 
00211                 for(j=pj;j<pj+dj;j++) {
00212 
00213                         recovered[i][j]=(float)((int)(recovered[i][j]+.5));
00214                 }
00215         }
00216 
00217         // calculate mse.
00218         my_mse=mse_region(original,recovered,pi,pj,di,dj);
00219 
00220         return(my_mse);
00221 }

00090 {
00091         int i,j,limi=pi+di,limj=pj+dj;
00092         float my_min=MAX,my_max=MIN;
00093 
00094         for(i=0;i<Ni;i++) {
00095 
00096                 for(j=0;j<Nj;j++) {
00097 
00098                         if((i>=pi)&&(i<limi)&&(j>=pj)&&(j<limj))
00099                                 continue;
00100                         if(im[i][j]<my_min)
00101                                 my_min=im[i][j];
00102                         if(im[i][j]>my_max)
00103                                 my_max=im[i][j];
00104                 }
00105         }
00106         *max=my_max;
00107         *min=my_min;
00108 }

00230 {
00231         int i,j;
00232         float **original,**recovered;
00233         int di=MISSING_BL_SIZE,dj=MISSING_BL_SIZE;
00234         int pi,pj,Ni,Nj;
00235 
00236         float mean,var,threshold,my_mse;
00237 
00238 
00239         // argument processing with rudimentary checking.
00240         if(argc!=8) {
00241 
00242                 printf(USAGE);
00243                 exit(1);
00244         }
00245 
00246         Ni=atoi(argv[3]);
00247         Nj=atoi(argv[4]);
00248         if((Ni<0)||(Nj<0)) {
00249 
00250                 printf("Negative image dimensions.\n");
00251                 exit(1);
00252         }
00253 
00254         pi=atoi(argv[5]);
00255         pj=atoi(argv[6]);
00256 
00257         if((pi<0)||(pi+di>=Ni)||(pj<0)||(pj+dj>=Nj)) {
00258 
00259                 printf("Missing block is out of bounds.\n");
00260                 exit(1);
00261         }
00262 
00263         DCT_SIZE=atoi(argv[7]);
00264 
00265         if(DCT_SIZE<0) {
00266 
00267                 printf("DCT_SIZE is negative.\n");
00268                 exit(1);
00269         }
00270 
00271         printf("------------------------------------------------------------------------------\n\n");   
00272         printf("Welcome to the image recovery proof of concept demo by Onur G. Guleryuz.\n");
00273         printf("This code is intended for research purposes.\n");
00274         printf("It is not implemented efficiently and \n");
00275         printf("it is not intended to be optimal or bulletproof.\n");
00276         printf("Please contact me at my academic email if you discover problems.\n\n");
00277         printf("Please experiment with it by changing the variables that control its operation\n");
00278         printf("on different types of image regions. The main algorithm is explained in:\n");
00279         printf("Onur G. Guleryuz, ``Nonlinear Approximation Based Image Recovery Using \n");
00280         printf("Adaptive Sparse Reconstructions,'' Proc. IEEE Int'l Conf. on Image Proc. \n");
00281         printf("(ICIP2003), Barcelona, Spain, Sept. 2003. (Please check my website for a \n");
00282         printf("journal paper that explains the algorithm in greater detail.) \n\n");
00283         printf("------------------------------------------------------------------------------\n");     
00284 
00285 
00286         if((pi<DCT_SIZE)||(pi+MISSING_BL_SIZE>Ni-DCT_SIZE)||\
00287                 (pj<DCT_SIZE)||(pj+MISSING_BL_SIZE>Nj-DCT_SIZE)) {
00288 
00289                 printf("WARNING: It is best to have a DCT_SIZE-1 pixel buffer from the boundary \n");
00290                 printf("         for good prediction results with the current layering strategy.\n\n");
00291         }
00292         
00293         // load original        
00294         original=load_image(argv[1],Ni,Nj);
00295         recovered=allocate_2d_float(Ni,Nj,0);
00296 
00297         for(i=0;i<Ni;i++) {
00298 
00299                 for(j=0;j<Nj;j++) {
00300 
00301                         recovered[i][j]=original[i][j];
00302                 }
00303         }
00304 
00305         // calculate surrounding mean and var to initialize missing region
00306         // and T_0
00307         mean=find_surrounding_mean(original,Ni,Nj,pi,pj,di,dj,1,1);
00308 
00309         var=find_surrounding_variance(original,Ni,Nj,pi,pj,di,dj,1,1,mean);
00310 
00311         // set T_0.
00312         threshold=(float)(sqrt(var));
00313 
00314         // We are doing very simplistic T_0 estimation. If T_0 is set too low
00315         // then algorithm will not perform well. One can get rid of the T_0 estimation
00316         // and set T_0 to 255 at the expense of some iterations (estimation hopefully
00317         // results in a T_0 that is just high enough).
00318         printf("Threshold = %8.3f, if too low please set manually or \n",threshold);
00319         printf("           increase the neighborhood size for variance computation.\n\n");
00320 
00321 
00322         // round mean, assuming a positive number.
00323         mean=(int)(mean+.5);
00324 
00325         // fill missing region with estimated mean value
00326         for(i=pi;i<pi+di;i++) {
00327 
00328                 for(j=pj;j<pj+dj;j++) {
00329 
00330                         recovered[i][j]=mean;
00331                 }
00332         }
00333 
00334         printf("Region at row=%5d column=%5d, size %5d x %5d\n\n",pi,pj,di,dj);
00335 
00336         // PSNR of mean filled image.
00337         my_mse=mse_region(original,recovered,pi,pj,di,dj);
00338         printf("Mean filled region, PSNR %8.3f dB\n",psnr(my_mse));
00339 
00340         save_image("mean_filled.raw",recovered,Ni,Nj);
00341         printf("(saved in image mean_filled.raw)\n\n");
00342 
00343         // fill missing region.
00344         my_mse=fill_layers(original,recovered,Ni,Nj,pi,pj,threshold);
00345         printf("Recovered region, PSNR %8.3f dB\n\n",psnr(my_mse));
00346 
00347         // save recovered image.
00348         save_image(argv[2],recovered,Ni,Nj);
00349 
00350         return(0);
00351 }

00048 {
00049         int i,j,limi=pi+di,limj=pj+dj,norm;
00050         float tmp1=0,tmp2;
00051 
00052         for(i=pi;i<limi;i++) {
00053 
00054                 for(j=pj;j<limj;j++) {
00055 
00056                         tmp2=im1[i][j]-im2[i][j];
00057                         tmp1+=tmp2*tmp2;
00058                 }
00059         }
00060 
00061         norm=di*dj;
00062         if(norm)
00063                 tmp1/=norm;
00064 
00065         return(tmp1);
00066 }

Defines
#define	psnr(mse) (10*log10(65025.0/(mse)))
#define	MIN ((float)-1e30)
#define	MAX ((float)1e30)
#define	MISSING_BL_SIZE 16
#define	NUM_LAYERS 8
#define	INNER_ITER_COUNT 1
#define	MIN_ITER_COUNT 10
#define	MAX_ITER_COUNT 10000
#define	PRINT_ITERS 1
#define	USAGE "recover_demo <raw_gray_scale_image_name> <raw_gray_scale_filled_image_name> <num_rows> <num_columns> <lost_block_row> <lost_block_col> <dct_size>\n"
Functions
float	mse_region (float im1, float im2, int pi, int pj, int di, int dj)
void	clip_to_min_max_region (float **im, int pi, int pj, int di, int dj, float min, float max)
void	find_max_min_avoid (float *im, int Ni, int Nj, int pi, int pj, int di, int dj, float max, float *min)
float	fill_layers (float original, float recovered, int Ni, int Nj, int pi, int pj, float threshold)
int	main (int argc, char **argv)
Variables
int	DCT_SIZE = 16

recover.c File Reference

Defines

Functions

Variables

Define Documentation

Function Documentation

Variable Documentation