Image Interpolation with Contour Stencils
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
imcoarsen.c
Go to the documentation of this file.
1 
20 #include <math.h>
21 #include <string.h>
22 #include <ctype.h>
23 
24 #include "imageio.h"
25 
26 #define VERBOSE 0
27 
29 #define NUMSTDS 4
30 
31 
33 typedef struct
34 {
36  uint32_t *Data;
38  int Width;
40  int Height;
41 } image;
42 
43 typedef enum
44 {
50 
52 typedef struct
53 {
55  char *InputFile;
57  char *OutputFile;
59  int JpegQuality;
65  float ScaleFactor;
67  float PsfSigma;
69 
70 
71 int ParseParams(programparams *Param, int argc, char *argv[]);
72 int Coarsen(image v, image u, programparams Param);
73 
76 {
77  puts("Image coarsening utility, P. Getreuer 2010-2011\n");
78  puts("Usage: imcoarsen [options] <input file> <output file>\n\n"
79  "Only " READIMAGE_FORMATS_SUPPORTED " images are supported.\n");
80  puts("Options:");
81  puts(" -x <number> the coarsening factor (>= 1.0, may be non-integer)");
82  puts(" -p <number> sigma_h, blur size of the Gaussian point spread");
83  puts(" function in units of output pixels.");
84  puts(" -b <ext> boundary handling extension, choices for <ext> are");
85  puts(" const constant extension");
86  puts(" hsym half-sample symmetric (default)");
87  puts(" wsym whole-sample symmetric");
88  puts(" per periodic");
89  puts(" -g <grid> grid to use for resampling, choices for <grid> are");
90  puts(" centered grid with centered alignment (default)");
91  puts(" topleft the top-left anchored grid\n");
92 #ifdef USE_LIBJPEG
93  puts(" -q <number> quality for saving JPEG images (0 to 100)\n");
94 #endif
95  puts("Example: coarsen by factor 2.5");
96  puts(" imcoarsen -x 2.5 -p 0.35 frog.bmp coarse.bmp");
97 }
98 
99 
100 int main(int argc, char *argv[])
101 {
102  programparams Param;
103  image u = {NULL, 0, 0}, v = {NULL, 0, 0};
104  int Status = 1;
105 
106 
107  if(!ParseParams(&Param, argc, argv))
108  return 0;
109 
110  /* Read the input image */
111  if(!(u.Data = (uint32_t *)ReadImage(&u.Width, &u.Height, Param.InputFile,
113  goto Catch;
114 
115  if(Param.ScaleFactor >= u.Width || Param.ScaleFactor >= u.Height)
116  {
117  ErrorMessage("Image is too small for scale factor.\n");
118  goto Catch;
119  }
120 
121  /* Allocate the output image */
122  v.Width = (int)ceil(u.Width / Param.ScaleFactor);
123  v.Height = (int)ceil(u.Height / Param.ScaleFactor);
124 #if VERBOSE > 0
125  printf("%dx%d input -> %dx%d output\n", u.Width, u.Height, v.Width, v.Height);
126 #endif
127 
128  if(!(v.Data = (uint32_t *)Malloc(sizeof(uint32_t)*
129  ((long int)v.Width)*((long int)v.Height))))
130  goto Catch;
131 
132  /* Convolution followed by downsampling */
133  if(!Coarsen(v, u, Param))
134  goto Catch;
135 
136  /* Write the output image */
137  if(!WriteImage(v.Data, v.Width, v.Height, Param.OutputFile,
139  goto Catch;
140 #if VERBOSE > 0
141  else
142  printf("Output written to \"%s\".\n", Param.OutputFile);
143 #endif
144 
145  Status = 0; /* Finished successfully, set exit status to zero. */
146 
147 Catch:
148  Free(v.Data);
149  Free(u.Data);
150  return Status;
151 }
152 
153 
154 float Sqr(float x)
155 {
156  return x*x;
157 }
158 
165 static int ConstExtension(int N, int i)
166 {
167  if(i < 0)
168  return 0;
169  else if(i >= N)
170  return N - 1;
171  else
172  return i;
173 }
174 
175 
182 static int HSymExtension(int N, int i)
183 {
184  while(1)
185  {
186  if(i < 0)
187  i = -1 - i;
188  else if(i >= N)
189  i = (2*N - 1) - i;
190  else
191  return i;
192  }
193 }
194 
195 
202 static int WSymExtension(int N, int i)
203 {
204  while(1)
205  {
206  if(i < 0)
207  i = -i;
208  else if(i >= N)
209  i = (2*N - 2) - i;
210  else
211  return i;
212  }
213 }
214 
215 
222 static int PerExtension(int N, int i)
223 {
224  while(1)
225  {
226  if(i < 0)
227  i += N;
228  else if(i >= N)
229  i -= N;
230  else
231  return i;
232  }
233 }
234 
235 
236 int (*ExtensionMethod[4])(int, int) =
237  {ConstExtension, HSymExtension, WSymExtension, PerExtension};
238 
239 
241 {
242  int (*Extension)(int, int) = ExtensionMethod[Param.Boundary];
243  const float PsfRadius = NUMSTDS*Param.PsfSigma*Param.ScaleFactor;
244  const int PsfWidth = 1 + ((PsfRadius == 0) ? 1 : (int)ceil(2*PsfRadius));
245  float *Temp = NULL, *PsfBuf = NULL;
246  float ExpDenom, Weight, Sum[4], DenomSum;
247  float XStart, YStart, X, Y;
248  uint32_t Pixel;
249  int IndexX0, IndexY0, SrcOffset, DestOffset;
250  int x, y, n, c, Success = 0;
251 
252 
253  if(!(Temp = (float *)Malloc(sizeof(float)*4*v.Width*u.Height))
254  || !(PsfBuf = (float *)Malloc(sizeof(float)*PsfWidth)))
255  goto Catch;
256 
257  ExpDenom = 2 * Sqr(Param.PsfSigma*Param.ScaleFactor);
258 
259  if(Param.CenteredGrid)
260  {
261  XStart = (1/Param.ScaleFactor - 1)/2;
262  YStart = (1/Param.ScaleFactor - 1)/2;
263  }
264  else
265  XStart = YStart = 0;
266 
267  for(x = 0; x < v.Width; x++)
268  {
269  X = (-XStart + x)*Param.ScaleFactor;
270  IndexX0 = (int)ceil(X - PsfRadius - 0.5f);
271  DenomSum = 0;
272 
273  /* Evaluate the PSF */
274  for(n = 0; n < PsfWidth; n++)
275  {
276  PsfBuf[n] = Sqr(X - (IndexX0 + n));
277 
278  if(!n || PsfBuf[n] < DenomSum)
279  DenomSum = PsfBuf[n];
280  }
281 
282  if(ExpDenom > 0)
283  for(n = 0; n < PsfWidth; n++)
284  PsfBuf[n] = (float)exp((DenomSum - PsfBuf[n]) / ExpDenom);
285  else if(IndexX0 == (int)floor(X - PsfRadius + 0.5f))
286  { /* If PsfSigma = 0, sample the nearest neighbor */
287  PsfBuf[0] = 1;
288  PsfBuf[1] = 0;
289  }
290  else /* At a half integer, average the two nearest neighbors */
291  PsfBuf[0] = PsfBuf[1] = 0.5f;
292 
293  DenomSum = 0;
294 
295  for(n = 0; n < PsfWidth; n++)
296  DenomSum += PsfBuf[n];
297 
298  for(y = 0, SrcOffset = 0, DestOffset = x; y < u.Height;
299  y++, SrcOffset += u.Width, DestOffset += v.Width)
300  {
301  Sum[0] = Sum[1] = Sum[2] = Sum[3] = 0;
302 
303  for(n = 0; n < PsfWidth; n++)
304  {
305  Weight = PsfBuf[n];
306  Pixel = u.Data[Extension(u.Width, IndexX0 + n) + SrcOffset];
307 
308  for(c = 0; c < 4; c++)
309  Sum[c] += (float)((uint8_t *)&Pixel)[c] * Weight;
310  }
311 
312  for(c = 0; c < 4; c++)
313  Temp[4*DestOffset + c] = Sum[c] / DenomSum;
314  }
315  }
316 
317  for(y = 0; y < v.Height; y++, v.Data += v.Width)
318  {
319  Y = (-YStart + y)*Param.ScaleFactor;
320  IndexY0 = (int)ceil(Y - PsfRadius - 0.5f);
321  DenomSum = 0;
322 
323  /* Evaluate the PSF */
324  for(n = 0; n < PsfWidth; n++)
325  {
326  PsfBuf[n] = Sqr(Y - (IndexY0 + n));
327 
328  if(!n || PsfBuf[n] < DenomSum)
329  DenomSum = PsfBuf[n];
330  }
331 
332  if(ExpDenom > 0)
333  for(n = 0; n < PsfWidth; n++)
334  PsfBuf[n] = (float)exp((DenomSum - PsfBuf[n]) / ExpDenom);
335  else if(IndexY0 == (int)floor(Y - PsfRadius + 0.5f))
336  { /* If PsfSigma = 0, sample the nearest neighbor */
337  PsfBuf[0] = 1;
338  PsfBuf[1] = 0;
339  }
340  else /* At a half integer, average the two nearest neighbors */
341  PsfBuf[0] = PsfBuf[1] = 0.5f;
342 
343  DenomSum = 0;
344 
345  for(n = 0; n < PsfWidth; n++)
346  DenomSum += PsfBuf[n];
347 
348  for(x = 0; x < v.Width; x++)
349  {
350  Sum[0] = Sum[1] = Sum[2] = Sum[3] = 0;
351 
352  for(n = 0; n < PsfWidth; n++)
353  {
354  SrcOffset = x + v.Width*Extension(u.Height, IndexY0 + n);
355  Weight = PsfBuf[n];
356 
357  for(c = 0; c < 4; c++)
358  Sum[c] += Temp[4*SrcOffset + c] * Weight;
359  }
360 
361  for(c = 0; c < 4; c++)
362  ((uint8_t *)&Pixel)[c] = (int)(Sum[c] / DenomSum + 0.5f);
363 
364  v.Data[x] = Pixel;
365  }
366  }
367 
368  Success = 1;
369 Catch:
370  Free(PsfBuf);
371  Free(Temp);
372  return Success;
373 }
374 
375 
376 int ParseParams(programparams *Param, int argc, char *argv[])
377 {
378  static char *DefaultOutputFile = (char *)"out.bmp";
379  char *OptionString;
380  char OptionChar;
381  int i;
382 
383 
384  if(argc < 2)
385  {
387  return 0;
388  }
389 
390  /* Set parameter defaults */
391  Param->InputFile = 0;
392  Param->OutputFile = DefaultOutputFile;
393  Param->JpegQuality = 99;
394  Param->ScaleFactor = 1;
395  Param->PsfSigma = 0.35f;
396  Param->CenteredGrid = 1;
397  Param->Boundary = BOUNDARY_HSYMMETRIC;
398 
399  for(i = 1; i < argc;)
400  {
401  if(argv[i] && argv[i][0] == '-')
402  {
403  if((OptionChar = argv[i][1]) == 0)
404  {
405  ErrorMessage("Invalid parameter format.\n");
406  return 0;
407  }
408 
409  if(argv[i][2])
410  OptionString = &argv[i][2];
411  else if(++i < argc)
412  OptionString = argv[i];
413  else
414  {
415  ErrorMessage("Invalid parameter format.\n");
416  return 0;
417  }
418 
419  switch(OptionChar)
420  {
421  case 'x':
422  Param->ScaleFactor = (float)atof(OptionString);
423 
424  if(Param->ScaleFactor < 1)
425  {
426  ErrorMessage("Invalid scale factor.\n");
427  return 0;
428  }
429  break;
430  case 'p':
431  Param->PsfSigma = (float)atof(OptionString);
432 
433  if(Param->PsfSigma < 0.0)
434  {
435  ErrorMessage("Point spread blur size must be nonnegative.\n");
436  return 0;
437  }
438  break;
439  case 'b':
440  if(!strcmp(OptionString, "const"))
441  Param->Boundary = BOUNDARY_CONSTANT;
442  else if(!strcmp(OptionString, "hsym"))
443  Param->Boundary = BOUNDARY_HSYMMETRIC;
444  else if(!strcmp(OptionString, "wsym"))
445  Param->Boundary = BOUNDARY_WSYMMETRIC;
446  else if(!strcmp(OptionString, "per"))
447  Param->Boundary = BOUNDARY_PERIODIC;
448  else
449  {
450  ErrorMessage("Boundary extension must be either \"const\", \"hsym\", or \"wsym\".\n");
451  return 0;
452  }
453  break;
454  case 'g':
455  if(!strcmp(OptionString, "centered")
456  || !strcmp(OptionString, "center"))
457  Param->CenteredGrid = 1;
458  else if(!strcmp(OptionString, "topleft")
459  || !strcmp(OptionString, "top-left"))
460  Param->CenteredGrid = 0;
461  else
462  {
463  ErrorMessage("Grid must be either \"centered\" or \"topleft\".\n");
464  return 0;
465  }
466  break;
467 
468 #ifdef USE_LIBJPEG
469  case 'q':
470  Param->JpegQuality = atoi(OptionString);
471 
472  if(Param->JpegQuality <= 0 || Param->JpegQuality > 100)
473  {
474  ErrorMessage("JPEG quality must be between 0 and 100.\n");
475  return 0;
476  }
477  break;
478 #endif
479  case '-':
481  return 0;
482  default:
483  if(isprint(OptionChar))
484  ErrorMessage("Unknown option \"-%c\".\n", OptionChar);
485  else
486  ErrorMessage("Unknown option.\n");
487 
488  return 0;
489  }
490 
491  i++;
492  }
493  else
494  {
495  if(!Param->InputFile)
496  Param->InputFile = argv[i];
497  else
498  Param->OutputFile = argv[i];
499 
500  i++;
501  }
502  }
503 
504  if(!Param->InputFile)
505  {
507  return 0;
508  }
509 
510  return 1;
511 }