Chan-Vese Segmentation
chanvesecli.c
Go to the documentation of this file.
1 
20 #include <math.h>
21 #include <stdlib.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include "cliio.h"
25 #include "chanvese.h"
26 #include "gifwrite.h"
27 #include "rgb2ind.h"
28 
29 #define ROUNDCLAMP(x) ((x < 0) ? 0 : \
30  ((x > 1) ? 255 : (uint8_t)floor(255.0*(x) + 0.5)))
31 
32 #ifdef __GNUC__
33 
34  #define ATTRIBUTE_UNUSED __attribute__((unused))
35 #else
36  #define ATTRIBUTE_UNUSED
37 #endif
38 
40 typedef struct
41 {
43  const char *InputFile;
45  const char *OutputFile;
47  const char *OutputFile2;
49  int JpegQuality;
50 
52  image Phi;
54  chanveseopt *Opt;
55 
56  int IterPerFrame;
57 } programparams;
58 
60 typedef struct
61 {
62  const num *Image;
63  unsigned char *Plot;
64  int *Delays;
65  int IterPerFrame;
66  int NumFrames;
67 } plotparam;
68 
69 
70 static void PrintHelpMessage()
71 {
72  puts("chanvese, P. Getreuer 2011-2012\n"
73  "Chan-Vese segmentation IPOL demo\n\n"
74  "Usage: chanvese [param:value ...] input animation final \n\n"
75  "where \"input\" and \"final\" are "
76  READIMAGE_FORMATS_SUPPORTED " files\n"
77  "and \"animation\" is a GIF file.\n");
78  puts("Parameters\n");
79  puts(" mu:<number> length penalty (default 0.25)");
80  puts(" nu:<number> area penalty (default 0.0)");
81  puts(" lambda1:<number> fit weight inside the cuve (default 1.0)");
82  puts(" lambda2:<number> fit weight outside the curve (default 1.0)");
83  puts(" phi0:<file> read initial level set from an image or text file");
84  puts(" tol:<number> convergence tolerance (default 1e-3)");
85  puts(" maxiter:<number> maximum number of iterations (default 500)");
86  puts(" dt:<number> time step (default 0.5)\n");
87  puts(" iterperframe:<number> iterations per frame (default 10)\n");
88 #ifdef LIBJPEG_SUPPORT
89  puts(" jpegquality:<number> Quality for saving JPEG images (0 to 100)\n");
90 #endif
91  puts("Example:\n"
92 #ifdef LIBPNG_SUPPORT
93  " chanvese tol:1e-5 mu:0.5 input.png animation.gif final.png\n");
94 #else
95  " chanvese tol:1e-5 mu:0.5 input.bmp animation.gif final.bmp\n");
96 #endif
97 }
98 
99 
100 static int PlotFun(int State, int Iter, num Delta,
101  const num *c1, const num *c2, const num *Phi,
102  int Width, int Height, int NumChannels, void *ParamPtr);
103 static int ParseParam(programparams *Param, int argc, const char *argv[]);
104 static int PhiRescale(image *Phi);
105 
106 
107 int WriteBinary(image Phi, const char *File)
108 {
109  unsigned char *Temp = NULL;
110  const int NumPixels = Phi.Width*Phi.Height;
111  int i, Success;
112 
113  if(!(Temp = (unsigned char *)malloc(Phi.Width*Phi.Height)))
114  return 0;
115 
116  for(i = 0; i < NumPixels; i++)
117  Temp[i] = (Phi.Data[i] >= 0) ? 255 : 0;
118 
119  Success = WriteImage(Temp, Phi.Width, Phi.Height, File,
120  IMAGEIO_U8 | IMAGEIO_GRAYSCALE, 0);
121 
122  free(Temp);
123  return Success;
124 }
125 
126 
127 int WriteAnimation(plotparam *PlotParam, int Width, int Height,
128  const char *OutputFile)
129 {
130  const int NumPixels = Width*Height;
131  unsigned char *PlotInd = NULL;
132  unsigned char *Palette = NULL;
133  unsigned char **PlotIndFrames = NULL;
134  int Frame, Success = 0;
135 
136  if(!(PlotInd = (unsigned char *)malloc(Width*Height*PlotParam->NumFrames))
137  || !(Palette = (unsigned char *)calloc(3*256, 1))
138  || !(PlotIndFrames = (unsigned char **)
139  malloc(sizeof(unsigned char *)*PlotParam->NumFrames)))
140  goto Catch;
141 
142  for(Frame = 0; Frame < PlotParam->NumFrames; Frame++)
143  PlotIndFrames[Frame] = PlotInd + NumPixels*Frame;
144 
145  /* Quantize colors for GIF */
146  if(!Rgb2Ind(PlotInd, Palette, 255, PlotParam->Plot,
147  Width*Height*PlotParam->NumFrames))
148  goto Catch;
149 
150  /* Optimize animation */
151  FrameDifference(PlotIndFrames, Width, Height, PlotParam->NumFrames, 255);
152 
153  /* Write the output animation */
154  if(!GifWrite(PlotIndFrames, Width, Height, PlotParam->NumFrames,
155  Palette, 256, 255, PlotParam->Delays, OutputFile))
156  {
157  fprintf(stderr, "Error writing \"%s\".\n", OutputFile);
158  goto Catch;
159  }
160  else
161  printf("Output written to \"%s\".\n", OutputFile);
162 
163  Success = 1;
164 Catch:
165  if(PlotIndFrames)
166  free(PlotIndFrames);
167  if(Palette)
168  free(Palette);
169  if(PlotInd)
170  free(PlotInd);
171  return Success;
172 }
173 
174 
175 int main(int argc, char *argv[])
176 {
177  programparams Param;
178  plotparam PlotParam;
179  image f = NullImage;
180  num c1[3], c2[3];
181  int Status = 1;
182 
183  PlotParam.Plot = NULL;
184  PlotParam.Delays = NULL;
185 
186  if(!ParseParam(&Param, argc, (const char **)argv))
187  goto Catch;
188 
189  /* Read the input image */
190  if(!ReadImageObj(&f, Param.InputFile))
191  goto Catch;
192 
193  if(Param.Phi.Data &&
194  (f.Width != Param.Phi.Width || f.Height != Param.Phi.Height))
195  {
196  fprintf(stderr, "Size mismatch: "
197  "phi0 (%dx%d) does not match image size (%dx%d).\n",
198  Param.Phi.Width, Param.Phi.Height, f.Width, f.Height);
199  goto Catch;
200  }
201 
202  PlotParam.Image = f.Data;
203  PlotParam.IterPerFrame = Param.IterPerFrame;
204  PlotParam.NumFrames = 0;
205 
206  ChanVeseSetPlotFun(Param.Opt, PlotFun, (void *)&PlotParam);
207 
208  printf("Segmentation parameters\n");
209  printf("f : [%d x %d %s]\n",
210  f.Width, f.Height, (f.NumChannels == 1) ? "grayscale" : "RGB");
211  printf("phi0 : %s\n", (Param.Phi.Data) ? "custom" : "default");
212  ChanVesePrintOpt(Param.Opt);
213 #ifdef NUM_SINGLE
214  printf("datatype : single precision float\n");
215 #else
216  printf("datatype : double precision float\n");
217 #endif
218  printf("\n");
219 
220  if(!Param.Phi.Data)
221  {
222  if(!AllocImageObj(&Param.Phi, f.Width, f.Height, 1))
223  {
224  fprintf(stderr, "Out of memory.");
225  goto Catch;
226  }
227 
228  ChanVeseInitPhi(Param.Phi.Data, Param.Phi.Width, Param.Phi.Height);
229  }
230 
231  /* Perform the segmentation */
232  if(!ChanVese(Param.Phi.Data, f.Data,
233  f.Width, f.Height, f.NumChannels, Param.Opt))
234  {
235  fprintf(stderr, "Error in ChanVese.");
236  goto Catch;
237  }
238 
239  /* Compute the final region averages */
240  RegionAverages(c1, c2, Param.Phi.Data, f.Data,
241  f.Width, f.Height, f.NumChannels);
242 
243  printf("\nRegion averages\n");
244 
245  if(f.NumChannels == 1)
246  printf("c1 : %.4f\nc2 : %.4f\n\n", c1[0], c2[0]);
247  else if(f.NumChannels == 3)
248  printf("c1 : (%.4f, %.4f, %.4f)\nc2 : (%.4f, %.4f, %.4f)\n\n",
249  c1[0], c1[1], c1[2], c2[0], c2[1], c2[2]);
250 
251  if(Param.OutputFile2 && !WriteBinary(Param.Phi, Param.OutputFile2))
252  goto Catch;
253 
254  if(!WriteAnimation(&PlotParam, f.Width, f.Height, Param.OutputFile))
255  goto Catch;
256 
257  Status = 0;
258 Catch:
259  if(PlotParam.Plot)
260  free(PlotParam.Plot);
261  if(PlotParam.Delays)
262  free(PlotParam.Delays);
263  FreeImageObj(Param.Phi);
264  FreeImageObj(f);
265  ChanVeseFreeOpt(Param.Opt);
266  return Status;
267 }
268 
269 
270 /* Plot callback function */
271 static int PlotFun(int State, int Iter, ATTRIBUTE_UNUSED num Delta,
272  ATTRIBUTE_UNUSED const num *c1, ATTRIBUTE_UNUSED const num *c2,
273  const num *Phi, int Width, int Height,
274  ATTRIBUTE_UNUSED int NumChannels, void *ParamPtr)
275 {
276  const int NumPixels = Width*Height;
277  plotparam *PlotParam = (plotparam *)ParamPtr;
278  unsigned char *Plot, *Temp = NULL;
279  int *Delays = NULL;
280  num Red, Green, Blue, Alpha;
281  long i;
282  int x, y, Edge, NumFrames = PlotParam->NumFrames, Success = 0;
283  int il, ir, iu, id;
284 
285  switch(State)
286  {
287  case 0:
288  /* We print to stderr so that messages are displayed on the console
289  immediately, during the TvRestore computation. If we use stdout,
290  messages might be buffered and not displayed until after TvRestore
291  completes, which would defeat the point of having this real-time
292  plot callback. */
293  if(NumChannels == 1)
294  fprintf(stderr, " Iteration %4d Delta %7.4f c1 = %6.4f c2 = %6.4f\r",
295  Iter, Delta, *c1, *c2);
296  else
297  fprintf(stderr, " Iteration %4d Delta %7.4f\r", Iter, Delta);
298  break;
299  case 1: /* Converged successfully */
300  fprintf(stderr, "Converged in %d iterations. \n",
301  Iter);
302  break;
303  case 2: /* Maximum iterations exceeded */
304  fprintf(stderr, "Maximum number of iterations exceeded. \n");
305  break;
306  }
307 
308  if(State == 0 && (Iter % PlotParam->IterPerFrame) > 0)
309  return 1;
310 
311  if(!(Plot = (unsigned char *)realloc(PlotParam->Plot,
312  3*Width*Height*(PlotParam->NumFrames + 1)))
313  || !(Delays = (int *)realloc(PlotParam->Delays,
314  sizeof(int)*(PlotParam->NumFrames + 1)))
315  || !(Temp = (unsigned char *)malloc(Width*Height)))
316  {
317  fprintf(stderr, "Out of memory.\n");
318  goto Catch;
319  }
320 
321  PlotParam->Plot = Plot;
322  PlotParam->Delays = Delays;
323  Plot += 3*NumPixels*PlotParam->NumFrames;
324 
325  for(y = 0, i = 0; y < Height; y++)
326  for(x = 0; x < Width; x++, i++)
327  {
328  if(Phi[i] >= 0 &&
329  ((x > 0 && Phi[i - 1] < 0)
330  || (x + 1 < Width && Phi[i + 1] < 0)
331  || (y > 0 && Phi[i - Width] < 0)
332  || (y + 1 < Height && Phi[i + Width] < 0)))
333  Edge = 1; /* Inside the curve, on the edge */
334  else
335  Edge = 0;
336 
337  Temp[i] = Edge;
338  }
339 
340  for(y = 0, i = 0; y < Height; y++)
341  {
342  iu = (y == 0) ? 0 : -Width;
343  id = (y == Height - 1) ? 0 : Width;
344 
345  for(x = 0; x < Width; x++, i++)
346  {
347  il = (x == 0) ? 0 : -1;
348  ir = (x == Width - 1) ? 0 : 1;
349 
350  Red = PlotParam->Image[i];
351  Green = PlotParam->Image[i + NumPixels];
352  Blue = PlotParam->Image[i + 2*NumPixels];
353 
354  Red = 0.95f*Red;
355  Green = 0.95f*Green;
356  Blue = 0.95f*Blue;
357 
358  i = x + Width*y;
359  Alpha = (4*Temp[i]
360  + Temp[i + ir] + Temp[i + il]
361  + Temp[i + id] + Temp[i + iu])/4.0f;
362 
363  if(Alpha > 1)
364  Alpha = 1;
365 
366  Red = (1 - Alpha)*Red;
367  Green = (1 - Alpha)*Green;
368  Blue = (1 - Alpha)*Blue + Alpha;
369 
370  Plot[3*i + 0] = ROUNDCLAMP(Red);
371  Plot[3*i + 1] = ROUNDCLAMP(Green);
372  Plot[3*i + 2] = ROUNDCLAMP(Blue);
373  }
374  }
375 
376  PlotParam->Delays[NumFrames] = (State == 0) ? 12 : 120;
377  PlotParam->NumFrames++;
378 
379  Success = 1;
380 Catch:
381  if(Temp)
382  free(Temp);
383  return Success;
384 }
385 
386 
387 static int ParseParam(programparams *Param, int argc, const char *argv[])
388 {
389  const char *Option, *Value;
390  num NumValue;
391  char TokenBuf[256];
392  int k, kread, Skip;
393 
394 
395  /* Set parameter defaults */
396  Param->InputFile = NULL;
397  Param->OutputFile = NULL;
398  Param->JpegQuality = 85;
399  Param->Phi = NullImage;
400  Param->Opt = NULL;
401  Param->IterPerFrame = 10;
402 
403  if(!(Param->Opt = ChanVeseNewOpt()))
404  {
405  fprintf(stderr, "Out of memory.\n");
406  return 0;
407  }
408 
409  if(argc < 2)
410  {
411  PrintHelpMessage();
412  return 0;
413  }
414 
415  k = 1;
416 
417  while(k < argc)
418  {
419  Skip = (argv[k][0] == '-') ? 1 : 0;
420  kread = CliParseArglist(&Option, &Value, TokenBuf, sizeof(TokenBuf),
421  k, &argv[k][Skip], argc, argv, ":");
422 
423  if(!Option)
424  {
425  if(!Param->InputFile)
426  Option = (char *)"f";
427  else
428  Option = (char *)"u";
429  }
430 
431  if(Option[0] == '-') /* Argument begins with two dashes "--" */
432  {
433  PrintHelpMessage();
434  return 0;
435  }
436 
437  if(!strcmp(Option, "f") || !strcmp(Option, "input"))
438  {
439  if(!Value)
440  {
441  fprintf(stderr, "Expected a value for option %s.\n", Option);
442  return 0;
443  }
444  Param->InputFile = Value;
445  }
446  else if(!strcmp(Option, "u") || !strcmp(Option, "output"))
447  {
448  if(!Value)
449  {
450  fprintf(stderr, "Expected a value for option %s.\n", Option);
451  return 0;
452  }
453  Param->OutputFile = Value;
454  }
455  else if(!strcmp(Option, "tol"))
456  {
457  if(CliGetNum(&NumValue, Value, Option))
458  ChanVeseSetTol(Param->Opt, NumValue);
459  else
460  return 0;
461  }
462  else if(!strcmp(Option, "mu"))
463  {
464  if(CliGetNum(&NumValue, Value, Option))
465  ChanVeseSetMu(Param->Opt, NumValue);
466  else
467  return 0;
468  }
469  else if(!strcmp(Option, "nu"))
470  {
471  if(CliGetNum(&NumValue, Value, Option))
472  ChanVeseSetNu(Param->Opt, NumValue);
473  else
474  return 0;
475  }
476  else if(!strcmp(Option, "lambda1"))
477  {
478  if(CliGetNum(&NumValue, Value, Option))
479  ChanVeseSetLambda1(Param->Opt, NumValue);
480  else
481  return 0;
482  }
483  else if(!strcmp(Option, "lambda2"))
484  {
485  if(CliGetNum(&NumValue, Value, Option))
486  ChanVeseSetLambda2(Param->Opt, NumValue);
487  else
488  return 0;
489  }
490  else if(!strcmp(Option, "dt"))
491  {
492  if(CliGetNum(&NumValue, Value, Option))
493  ChanVeseSetDt(Param->Opt, NumValue);
494  else
495  return 0;
496  }
497  else if(!strcmp(Option, "maxiter"))
498  {
499  if(CliGetNum(&NumValue, Value, Option))
500  ChanVeseSetMaxIter(Param->Opt, (int)NumValue);
501  else
502  return 0;
503  }
504  else if(!strcmp(Option, "phi0"))
505  {
506  if(!Value)
507  {
508  fprintf(stderr, "Expected a value for option %s.\n", Option);
509  return 0;
510  }
511 
512  if(Param->Phi.Data)
513  FreeImageObj(Param->Phi);
514 
515  if(!(ReadMatrixFromFile(&Param->Phi, Value, PhiRescale)))
516  return 0;
517  }
518  else if(!strcmp(Option, "jpegquality"))
519  {
520  if(!CliGetNum(&NumValue, Value, Option))
521  return 0;
522  else if(NumValue < 0 || 100 < NumValue)
523  {
524  fprintf(stderr, "JPEG quality must be between 0 and 100.\n");
525  return 0;
526  }
527  else
528  Param->JpegQuality = (int)NumValue;
529  }
530  else if(!strcmp(Option, "iterperframe"))
531  {
532  if(!CliGetNum(&NumValue, Value, Option))
533  return 0;
534  else if(NumValue <= 0)
535  {
536  fprintf(stderr, "Iterations per frame must be postive.\n");
537  return 0;
538  }
539  else
540  Param->IterPerFrame = (int)NumValue;
541  }
542  else if(Skip)
543  {
544  fprintf(stderr, "Unknown option \"%s\".\n", Option);
545  return 0;
546  }
547  else
548  {
549  if(!Param->InputFile)
550  Param->InputFile = argv[k];
551  else if(!Param->OutputFile)
552  Param->OutputFile = argv[k];
553  else
554  Param->OutputFile2 = argv[k];
555 
556  kread = k;
557  }
558 
559  k = kread + 1;
560  }
561 
562  if(!Param->InputFile || !Param->OutputFile)
563  {
564  PrintHelpMessage();
565  return 0;
566  }
567 
568  return 1;
569 }
570 
571 
572 /* If phi is read from an image file, this function is called to rescale
573  it from the range [0,1] to [-10,10]. */
574 static int PhiRescale(image *Phi)
575 {
576  const long NumEl = ((long)Phi->Width) * ((long)Phi->Height);
577  num *Data = Phi->Data;
578  long n;
579 
580  for(n = 0; n < NumEl; n++)
581  Data[n] = 4*(2*Data[n] - 1);
582 
583  return 1;
584 }
585 
Generated on Fri Jul 5 2013 18:57:19 for Chan-Vese Segmentation by   doxygen 1.8.3.1