Image Interpolation with Contour Stencils

imdiff.c

Go to the documentation of this file.
00001 
00032 #include <math.h>
00033 #include <string.h>
00034 #include <ctype.h>
00035 
00036 #include "imageio.h"
00037 #include "conv.h"
00038 
00039 
00041 #define DISPLAY_SCALING     255
00042 
00043 #define MSSIM_K1            0.01
00044 #define MSSIM_K2            0.03
00045 
00046 #define MSSIM_C1            (MSSIM_K1*MSSIM_K1)
00047 #define MSSIM_C2            (MSSIM_K2*MSSIM_K2)
00048 
00049 
00051 typedef enum {DEFAULT_METRICS, MAX_METRIC, MSE_METRIC, RMSE_METRIC,
00052     PSNR_METRIC, MSSIM_METRIC} metric;
00053 
00055 typedef struct
00056 {
00058     char *FileA;
00060     char *FileB;    
00062     int JpegQuality;
00064     metric Metric;
00066     int SeparateChannels;
00068     int Pad;
00069         
00071     char *DifferenceFile;
00073     float D;
00074 } programparams;    
00075     
00076 
00077 static int ParseParams(programparams *Param, int argc, char *argv[]);
00078 void MakeDifferenceImage(float *A, const float *B, 
00079     int Width, int Height, int NumChannels, float D);
00080 void BasicMetrics(float *Max, float *Mse, const float *A, const float *B, 
00081     int Width, int Height, int NumChannels, int Pad);
00082 float ComputeMssim(const float *A, const float *B, 
00083     int Width, int Height, int NumChannels, int Pad);
00084     
00085 
00087 void PrintHelpMessage()
00088 {
00089     printf("Image difference calculator, P. Getreuer 2010-2011\n\n");
00090     printf("Usage: imdiff [options] <exact file> <distorted file>\n\n"
00091         "Only " READIMAGE_FORMATS_SUPPORTED " images are supported.\n\n");
00092     printf("Options:\n");
00093     printf("   -m <metric>  metric to use for comparison, choices are\n");
00094     printf("        max     maximum absolute difference, max_n |A_n - B_n|\n");
00095     printf("        mse     mean squared error, 1/N sum |A_n - B_n|^2\n");
00096     printf("        rmse    root mean squared error, (MSE)^1/2\n");
00097     printf("        psnr    peak signal-to-noise ratio, -10 log10(MSE/255^2)\n");
00098     printf("        mssim   mean structural similarity\n\n");
00099     printf("   -s           compute metric separately for each channel\n");
00100     printf("   -p <pad>     remove a margin of <pad> pixels before comparison\n");
00101     printf("   -D <number>  D parameter for difference image (explained below)\n\n");    
00102 #ifdef LIBJPEG_SUPPORT
00103     printf("   -q <number>  quality for saving JPEG images (0 to 100)\n\n");
00104 #endif    
00105     printf("Alternatively, a difference image is generated by the syntax\n"
00106            "   imdiff [-D <number>] <exact file> <distorted file> <output file>\n\n");
00107     printf("The difference image is computed as\n"
00108            "   D_n = 255/D (A_n - B_n) + 255/2.\n"
00109            "Values outside of the range [0,255] are saturated.\n\n");
00110     printf("Example:\n"
00111 #ifdef LIBPNG_SUPPORT
00112         "   imdiff -mpsnr frog-exact.png frog-4x.bmp\n");
00113 #else
00114         "   imdiff -mpsnr frog-exact.bmp frog-4x.bmp\n");
00115 #endif
00116 }   
00117 
00118 
00119 int main(int argc, char *argv[])
00120 {
00121     struct
00122     {
00123         float *Data;
00124         int Width;
00125         int Height;
00126     } A = {NULL, 0, 0}, B = {NULL, 0, 0};
00127 
00128     programparams Param;
00129     float Max, MaxC[3], Mse, MseC[3], Mssim;
00130     int Channel, Status = 1;
00131     
00132            
00133     if(!ParseParams(&Param, argc, argv))
00134         return 0;
00135     
00136     /* Read the exact image */
00137     if(!(A.Data = (float *)ReadImage(&A.Width, &A.Height, Param.FileA, 
00138         IMAGEIO_FLOAT | IMAGEIO_RGB | IMAGEIO_PLANAR)))
00139         goto Catch;
00140     
00141     /* Read the distorted image */
00142     if(!(B.Data = (float *)ReadImage(&B.Width, &B.Height, Param.FileB, 
00143         IMAGEIO_FLOAT | IMAGEIO_RGB | IMAGEIO_PLANAR)))
00144         goto Catch;
00145     
00146     if(A.Width != B.Width || A.Height != B.Height)
00147     {
00148         ErrorMessage("Image sizes don't match, %dx%d vs. %dx%d.\n", 
00149             A.Width, A.Height, B.Width, B.Height);
00150         goto Catch;
00151     }
00152     else if(A.Width <= 2*Param.Pad || A.Height <= 2*Param.Pad)
00153     {
00154         ErrorMessage(
00155             "Removal of %d-pixel padding removes entire %dx%d image.\n",
00156             Param.Pad, A.Width, A.Height);
00157         goto Catch;
00158     }
00159     
00160     if(Param.DifferenceFile)
00161     {
00162         MakeDifferenceImage(A.Data, B.Data, A.Width, A.Height, 3, Param.D);
00163         
00164         if(!(WriteImage(A.Data, A.Width, A.Height, Param.DifferenceFile, 
00165             IMAGEIO_FLOAT | IMAGEIO_RGB | IMAGEIO_PLANAR, Param.JpegQuality)))
00166             goto Catch;
00167     }
00168     else
00169     {
00170         Max = 0;
00171         Mse = 0;
00172         
00173         for(Channel = 0; Channel < 3; Channel++)
00174         {
00175             BasicMetrics(&MaxC[Channel], &MseC[Channel], 
00176                     A.Data + Channel*A.Width*A.Height, 
00177                     B.Data + Channel*B.Width*B.Height, 
00178                     A.Width, A.Height, 1, Param.Pad);
00179            
00180             if(MaxC[Channel] > Max)
00181                 Max = MaxC[Channel];
00182             
00183             Mse += MseC[Channel];
00184         }
00185         
00186         Mse /= 3;
00187         
00188         switch(Param.Metric)
00189         {
00190             case DEFAULT_METRICS:
00191                 if(!Param.SeparateChannels)
00192                 {
00193                     printf("Maximum absolute difference:  %g\n", DISPLAY_SCALING*Max);
00194                     printf("Peak signal-to-noise ratio:   %.4f\n", -10*log10(Mse));
00195                 }
00196                 else
00197                 {
00198                     printf("Maximum absolute difference:  %g %g %g\n", 
00199                         DISPLAY_SCALING*MaxC[0], DISPLAY_SCALING*MaxC[1],
00200                         DISPLAY_SCALING*MaxC[2]);
00201                     printf("Peak signal-to-noise ratio:   %.4f %.4f %.4f\n", 
00202                         -10*log10(MseC[0]), -10*log10(MseC[1]), -10*log10(MseC[2]));
00203                 }
00204                 
00205                 if(A.Width <= 2*(5 + Param.Pad) 
00206                     || A.Height <= 2*(5 + Param.Pad))
00207                     printf("Image size is too small to compute MSSIM.\n");
00208                 else
00209                 {
00210                     
00211                     Mssim = (Max == 0) ? 1 : ComputeMssim(A.Data, B.Data, 
00212                             A.Width, A.Height, 3, Param.Pad);
00213                 
00214                     if(Mssim != -1)
00215                         printf("Mean structural similarity:   %.4f\n", Mssim);
00216                 }
00217                 break;
00218             case MAX_METRIC:
00219                 if(!Param.SeparateChannels)
00220                     printf("%g\n", DISPLAY_SCALING*Max);
00221                 else
00222                     printf("%g %g %g\n", DISPLAY_SCALING*MaxC[0], 
00223                         DISPLAY_SCALING*MaxC[1], DISPLAY_SCALING*MaxC[2]);
00224                 break;
00225             case MSE_METRIC:
00226                 if(!Param.SeparateChannels)
00227                     printf("%.4f\n", DISPLAY_SCALING*DISPLAY_SCALING*Mse);
00228                 else
00229                     printf("%.4f %.4f %.4f\n", 
00230                         DISPLAY_SCALING*DISPLAY_SCALING*MseC[0],
00231                         DISPLAY_SCALING*DISPLAY_SCALING*MseC[1],
00232                         DISPLAY_SCALING*DISPLAY_SCALING*MseC[2]);
00233                 break;
00234             case RMSE_METRIC:
00235                 if(!Param.SeparateChannels)
00236                     printf("%.4f\n", DISPLAY_SCALING*sqrt(Mse));
00237                 else
00238                     printf("%.4f %.4f %.4f\n", 
00239                         DISPLAY_SCALING*sqrt(MseC[0]),
00240                         DISPLAY_SCALING*sqrt(MseC[1]),
00241                         DISPLAY_SCALING*sqrt(MseC[2]));
00242                 break;
00243             case PSNR_METRIC:
00244                 if(!Param.SeparateChannels)
00245                     printf("%.4f\n", -10*log10(Mse));
00246                 else
00247                 printf("%.4f %.4f %.4f\n", 
00248                     -10*log10(MseC[0]), -10*log10(MseC[1]), -10*log10(MseC[2]));
00249                 break;            
00250             case MSSIM_METRIC:
00251                 if(A.Width <= 2*(5 + Param.Pad) 
00252                     || A.Height <= 2*(5 + Param.Pad))
00253                     printf("Image size is too small to compute MSSIM.\n");
00254                 else
00255                 {
00256                     Mssim = (Max == 0) ? 1 : ComputeMssim(A.Data, B.Data, 
00257                             A.Width, A.Height, 3, Param.Pad);
00258                 
00259                     if(Mssim == -1)
00260                         ErrorMessage("Memory allocation failed.\n");
00261                     else                
00262                         printf("%.4f\n", Mssim);
00263                 }
00264                 break;
00265         }
00266     }
00267     
00268     Status = 0; /* Finished successfully */
00269 Catch:
00270     Free(B.Data);
00271     Free(A.Data);
00272     return Status;
00273 }
00274 
00275 
00277 void MakeDifferenceImage(float *A, const float *B, 
00278     int Width, int Height, int NumChannels, float D)
00279 {
00280     const int NumEl = NumChannels*Width*Height;
00281     int n;
00282     
00283     D /= 255;
00284     
00285     for(n = 0; n < NumEl; n++)
00286         A[n] = (A[n] - B[n])/D + (float)0.5;
00287 }
00288 
00289 
00291 void BasicMetrics(float *Max, float *Mse, const float *A, const float *B, 
00292     int Width, int Height, int NumChannels, int Pad)
00293 {
00294     float Diff, CurMax = 0;
00295     double AccumMse = 0;
00296     int x, y, Channel, n;
00297     
00298     
00299     for(Channel = 0; Channel < NumChannels; Channel++)
00300         for(y = Pad; y < Height - Pad; y++)
00301             for(x = Pad; x < Width - Pad; x++)
00302             {
00303                 n = x + Width*(y + Height*Channel);                
00304                 Diff = (float)fabs(A[n] - B[n]);
00305                     
00306                 if(CurMax < Diff)
00307                     CurMax = Diff;
00308                 
00309                 AccumMse += Diff*Diff;
00310             }
00311     
00312     *Max = CurMax;
00313     *Mse = (float)(AccumMse / (NumChannels*(Width - 2*Pad)*(Height - 2*Pad)));
00314 }
00315 
00316 
00318 float ComputeMssim(const float *A, const float *B, 
00319     int Width, int Height, int NumChannels, int Pad)
00320 {   
00321     /* 11-tap Gaussian filter with standard deviation 1.5 */
00322     const int R = 5;
00323     filter Window = GaussianFilter(1.5, R);    
00324     /* Boundary does not matter, convolution is used only in the interior */
00325     boundaryext Boundary = GetBoundaryExt("zpd");
00326     
00327     const int NumPixels = Width*Height;
00328     const int NumEl = NumChannels*NumPixels;
00329     float *Buffer = NULL, *MuA = NULL, *MuB = NULL, 
00330         *MuAA = NULL, *MuBB = NULL, *MuAB = NULL;
00331     double MuASqr, MuBSqr, MuAMuB, 
00332         SigmaASqr, SigmaBSqr, SigmaAB, Mssim = -1;
00333     int n, x, y, c; 
00334     
00335     
00336     if(IsNullFilter(Window)
00337         || !(Buffer = (float *)Malloc(sizeof(float)*NumPixels))
00338         || !(MuA = (float *)Malloc(sizeof(float)*NumEl))
00339         || !(MuB = (float *)Malloc(sizeof(float)*NumEl))
00340         || !(MuAA = (float *)Malloc(sizeof(float)*NumEl))
00341         || !(MuBB = (float *)Malloc(sizeof(float)*NumEl))
00342         || !(MuAB = (float *)Malloc(sizeof(float)*NumEl)))
00343         goto Catch;
00344     
00345     SeparableConv2D(MuA, Buffer, A, Window, Window, 
00346         Boundary, Width, Height, NumChannels);        
00347     SeparableConv2D(MuB, Buffer, B, Window, Window, 
00348         Boundary, Width, Height, NumChannels);
00349     
00350     for(n = 0; n < NumEl; n++)
00351     {        
00352         MuAA[n] = A[n]*A[n];
00353         MuBB[n] = B[n]*B[n];
00354         MuAB[n] = A[n]*B[n];
00355     }
00356     
00357     SeparableConv2D(MuAA, Buffer, MuAA, Window, Window, 
00358         Boundary, Width, Height, NumChannels);
00359     SeparableConv2D(MuBB, Buffer, MuBB, Window, Window, 
00360         Boundary, Width, Height, NumChannels);
00361     SeparableConv2D(MuAB, Buffer, MuAB, Window, Window, 
00362         Boundary, Width, Height, NumChannels);
00363     Mssim = 0;
00364     
00365     Pad += R;
00366     
00367     for(c = 0; c < NumChannels; c++)
00368         for(y = Pad; y < Height - Pad; y++)
00369             for(x = Pad; x < Width - Pad; x++)            
00370             {
00371                 n = x + Width*(y + Height*c);
00372                 MuASqr = MuA[n]*MuA[n];
00373                 MuBSqr = MuB[n]*MuB[n];
00374                 MuAMuB = MuA[n]*MuB[n];
00375                 SigmaASqr = MuAA[n] - MuASqr;
00376                 SigmaBSqr = MuBB[n] - MuBSqr;
00377                 SigmaAB = MuAB[n] - MuAMuB;
00378                 
00379                 Mssim += ((2*MuAMuB + MSSIM_C1)*(2*SigmaAB + MSSIM_C2))
00380                     / ((MuASqr + MuBSqr + MSSIM_C1)
00381                         *(SigmaASqr + SigmaBSqr + MSSIM_C2));
00382             }
00383     
00384     Mssim /= NumChannels*(Width - 2*Pad)*(Height - 2*Pad);
00385     
00386 Catch:
00387     FreeFilter(Window);
00388     Free(MuAB);
00389     Free(MuBB);
00390     Free(MuAA);
00391     Free(MuB);
00392     Free(MuA);
00393     Free(Buffer);
00394     return (float)Mssim;
00395 }
00396 
00397 
00398 
00399 static int ParseParams(programparams *Param, int argc, char *argv[])
00400 {
00401     char *OptionString;
00402     char OptionChar;
00403     int i;
00404 
00405     
00406     if(argc < 2)
00407     {
00408         PrintHelpMessage();
00409         return 0;
00410     }
00411     
00412     /* Set parameter defaults */
00413     Param->FileA = NULL;
00414     Param->FileB = NULL;    
00415     Param->Metric = DEFAULT_METRICS;    
00416     Param->SeparateChannels = 0;
00417     
00418     Param->Pad = 0;    
00419     Param->DifferenceFile = NULL;
00420     Param->JpegQuality = 95;
00421     Param->D = 20;
00422         
00423     for(i = 1; i < argc;)
00424     {
00425         if(argv[i] && argv[i][0] == '-')
00426         {
00427             if((OptionChar = argv[i][1]) == 0)
00428             {
00429                 ErrorMessage("Invalid parameter format.\n");
00430                 return 0;
00431             }
00432 
00433             if(argv[i][2])
00434                 OptionString = &argv[i][2];
00435             else if(++i < argc)
00436                 OptionString = argv[i];
00437             else
00438             {
00439                 ErrorMessage("Invalid parameter format.\n");
00440                 return 0;
00441             }
00442             
00443             switch(OptionChar)
00444             {
00445             case 'p':
00446                 Param->Pad = atoi(OptionString);
00447                 
00448                 if(Param->Pad < 0)
00449                 {
00450                     ErrorMessage("Pad must be nonnegative.\n");
00451                     return 0;
00452                 }
00453                 break;
00454             case 's':
00455                 Param->SeparateChannels = 1;
00456                 i--;
00457                 break;
00458             case 'D':
00459                 Param->D = (float)atof(OptionString);
00460 
00461                 if(Param->D <= 0)
00462                 {
00463                     ErrorMessage("D must be positive.\n");
00464                     return 0;
00465                 }
00466                 break;
00467             case 'm':
00468                 if(!strcmp(OptionString, "max"))
00469                     Param->Metric = MAX_METRIC;
00470                 else if(!strcmp(OptionString, "mse"))
00471                     Param->Metric = MSE_METRIC;
00472                 else if(!strcmp(OptionString, "rmse"))
00473                     Param->Metric = RMSE_METRIC;
00474                 else if(!strcmp(OptionString, "psnr"))
00475                     Param->Metric = PSNR_METRIC;
00476                 else if(!strcmp(OptionString, "mssim"))
00477                     Param->Metric = MSSIM_METRIC;
00478                 else
00479                     ErrorMessage("Unknown metric.\n");
00480                 break;
00481                 
00482 #ifdef LIBJPEG_SUPPORT
00483             case 'q':
00484                 Param->JpegQuality = atoi(OptionString);
00485 
00486                 if(Param->JpegQuality <= 0 || Param->JpegQuality > 100)
00487                 {
00488                     ErrorMessage("JPEG quality must be between 0 and 100.\n");
00489                     return 0;
00490                 }
00491                 break;
00492 #endif
00493             case '-':
00494                 PrintHelpMessage();
00495                 return 0;
00496             default:
00497                 if(isprint(OptionChar))
00498                     ErrorMessage("Unknown option \"-%c\".\n", OptionChar);
00499                 else
00500                     ErrorMessage("Unknown option.\n");
00501 
00502                 return 0;
00503             }
00504 
00505             i++;
00506         }
00507         else
00508         {
00509             if(!Param->FileA)
00510                 Param->FileA = argv[i];
00511             else if(!Param->FileB)
00512                 Param->FileB = argv[i];
00513             else
00514                 Param->DifferenceFile = argv[i];
00515 
00516             i++;
00517         }
00518     }
00519     
00520     if(!Param->FileA || !Param->FileB)
00521     {
00522         PrintHelpMessage();
00523         return 0;
00524     }
00525     
00526     return 1;
00527 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines