DGtalTools  1.5.beta
at-u0-v1.cpp
1 
32 #include <iostream>
33 
34 #include <sstream>
35 #include <string>
36 #include <functional>
37 #include <boost/format.hpp>
38 
39 #include "CLI11.hpp"
40 
41 #include "DGtal/base/Common.h"
42 #include "DGtal/helpers/StdDefs.h"
43 #include "DGtal/io/readers/GenericReader.h"
44 #include "DGtal/io/writers/GenericWriter.h"
45 
46 #include "ATu0v1.h"
47 
187 using namespace std;
188 using namespace DGtal;
189 
190 int main( int argc, char* argv[] )
191 {
192  using namespace Z2i;
193 
194 
195  // parse command line using CLI ----------------------------------------------
196  CLI::App app;
197  string f1;
198  string f2 {"AT"};
199  string inpainting_mask;
200  double l;
201  double l1 {0.3125};
202  double l2 {0.0005};
203  double lr {sqrt(2)};
204  double a {1.0};
205  double epsilon;
206 
207  double e1 {2.0};
208  double e2 {0.25};
209  double er {2.0};
210  int verb {0};
211  int nbiter = {10};
212  int pix_sz = {1};
213  string scv {"0xff0000"};
214  string isnr;
215 
216 
217  stringstream ssDescr;
218  ssDescr << "Computes a piecewise smooth approximation of a grey-level or color image, by optimizing the Ambrosio-Tortorelli functional (with u a 0-form and v a 1-form).";
219  ssDescr << "Usage: " << argv[0] << " -i toto.pgm\n"
220  << "Computes the Ambrosio-Tortorelli reconstruction/segmentation of an input image."
221  << "It outputs 2 or 3 images (of basename given by option --output) giving the"
222  << " reconstructed image u, and other images superposing u and the discontinuities v."
223  << endl << endl
224  << " / "
225  << endl
226  << " | a.(u-g)^2 + v^2 |grad u|^2 + le.|grad v|^2 + (l/4e).(1-v)^2 "
227  << endl
228  << " / "
229  << endl
230  << "Discretized as (u 0-form, v 1-form, A vertex-edge bdry, B edge-face bdy)" << endl
231  << "E(u,v) = a(u-g)^t (u-g) + u^t A^t diag(v)^2 A^t u + l e v^t (A A^t + B^t B) v + l/(4e) (1-v)^t (1-v)" << endl
232  << endl
233  << "Example: ./at-u0-v1 -i ../Images/cerclesTriangle64b02.pgm -o tmp -a 0.05 -e 1 --lambda-1 0.1 --lambda-2 0.00001";
234  app.description(ssDescr.str());
235 
236 
237 
238  app.add_option("-i,--input,1", f1, "the input image PPM filename." )
239  ->required()
240  ->check(CLI::ExistingFile);
241  app.add_option("--inpainting-mask,-m", inpainting_mask, "the input inpainting mask filename." );
242  app.add_option("--output,-o", f2, "the output image basename.", true);
243  auto lambdaOpt = app.add_option("--lambda,-l",l, "the parameter lambda.");
244  app.add_option("--lambda-1,-1",l1, "the initial parameter lambda (l1).", true);
245  app.add_option("--lambda-2,-2",l2, "the final parameter lambda (l2).", true );
246  app.add_option("--lambda-ratio,-q",lr, "the division ratio for lambda from l1 to l2.", true);
247  app.add_option("--alpha,-a",a, "the parameter alpha.", true);
248  auto epsOpt = app.add_option("--epsilon,-e", "the initial and final parameter epsilon of AT functional at the same time.");
249 
250  app.add_option("--epsilon-1",e1, "the initial parameter epsilon.", true);
251  app.add_option("--epsilon-2",e2, "the final parameter epsilon.", true);
252  app.add_option("--epsilon-r",er, "sets the ratio between two consecutive epsilon values of AT functional.", true);
253 
254  app.add_option("--nbiter,-n",nbiter, "the maximum number of iterations.", true );
255  auto snrOpt = app.add_option("--image-snr", isnr, "the input image without deterioration if you wish to compute the SNR.");
256  app.add_option("--pixel-size,-p", pix_sz, "the pixel size for outputing images (useful when one wants to see the discontinuities v on top of u).", true);
257  app.add_option("--color-v,-c",scv, "the color chosen for displaying the singularities v (e.g. red is 0xff0000).", true );
258  app.add_option("--verbose,-v", verb, "the verbose level (0: silent, 1: less silent, etc).", true );
259 
260  app.get_formatter()->column_width(40);
261  CLI11_PARSE(app, argc, argv);
262  // END parse command line using CLI ----------------------------------------------
263 
264 
265 
266  Color color_v( (unsigned int) std::stoul( scv, nullptr, 16 ), 255 );
267  if ( lambdaOpt->count()) l1 = l2 = l;
268  if ( l2 > l1 ) l2 = l1;
269  if ( lr <= 1.0 ) lr = sqrt(2);
270  if ( epsOpt->count() > 0 ){
271  e1 = e2 = epsilon;
272 
273  }
274  bool snr = snrOpt->count() > 0;
275 
276 
277  bool color_image = f1.size() > 4 && f1.compare( f1.size() - 4, 4, ".ppm" ) == 0;
278  bool grey_image = f1.size() > 4 && f1.compare( f1.size() - 4, 4, ".pgm" ) == 0;
279  if ( ! color_image && ! grey_image )
280  {
281  trace.error() << "Input image file must be either a PGM (grey-level) or a PPM (color) image with these extensions."
282  << endl;
283  return 2;
284  }
285 
286  KSpace K;
287  ATu0v1< KSpace > AT( verb );
288  Domain domain;
289 
290  typedef ATu0v1<KSpace>::Calculus Calculus;
291  typedef ImageContainerBySTLVector<Domain, Color> ColorImage;
292  typedef ImageContainerBySTLVector<Domain, unsigned char> GreyLevelImage;
293  //---------------------------------------------------------------------------
294  if ( color_image )
295  {
296  trace.beginBlock("Reading PPM image");
297  ColorImage image = PPMReader<ColorImage>::importPPM( f1 );
298  trace.endBlock();
299  trace.beginBlock("Building AT");
300  domain = image.domain();
301  K.init( domain.lowerBound(), domain.upperBound() - Point::diagonal( 1 ), true );
302  AT.init( K );
303  AT.addInput( image, [] ( Color c ) -> double { return ((double) c.red()) / 255.0; } );
304  AT.addInput( image, [] ( Color c ) -> double { return ((double) c.green()) / 255.0; } );
305  AT.addInput( image, [] ( Color c ) -> double { return ((double) c.blue()) / 255.0; } );
306  trace.endBlock();
307  }
308  else if ( grey_image )
309  {
310  trace.beginBlock("Reading PGM image");
311  GreyLevelImage image = PGMReader<GreyLevelImage>::importPGM( f1 );
312  trace.endBlock();
313  trace.beginBlock("Building AT");
314  domain = image.domain();
315  K.init( domain.lowerBound(), domain.upperBound() - Point::diagonal( 1 ), true );
316  AT.init( K );
317  AT.addInput( image, [] (unsigned char c ) { return ((double) c) / 255.0; } );
318  trace.endBlock();
319  }
320 
321  //---------------------------------------------------------------------------
322  if ( snr && color_image )
323  {
324  trace.beginBlock("Reading ideal PPM image");
325  ColorImage image = PPMReader<ColorImage>::importPPM( isnr );
326  trace.endBlock();
327  AT.addInput( image, [] ( Color c ) -> double { return ((double) c.red()) / 255.0; }, true );
328  AT.addInput( image, [] ( Color c ) -> double { return ((double) c.green()) / 255.0; }, true );
329  AT.addInput( image, [] ( Color c ) -> double { return ((double) c.blue()) / 255.0; }, true );
330  }
331  else if ( snr && grey_image )
332  {
333  trace.beginBlock("Reading ideal PGM image");
334  GreyLevelImage image = PGMReader<GreyLevelImage>::importPGM( isnr );
335  trace.endBlock();
336  AT.addInput( image, [] (unsigned char c ) { return ((double) c) / 255.0; }, true );
337  }
338 
339  //---------------------------------------------------------------------------
340  // Prepare zoomed output domain
341  Domain out_domain( pix_sz * domain.lowerBound(),
342  pix_sz * domain.upperBound() + Point::diagonal( pix_sz - 1) );
343  //---------------------------------------------------------------------------
344  AT.setUFromInput();
345  double g_snr = snr ? AT.computeSNR() : 0.0;
346 
347  if ( inpainting_mask.size() > 0 )
348  {
349  string fm = inpainting_mask;
350  trace.beginBlock("Reading inpainting mask");
351  GreyLevelImage mask = GenericReader<GreyLevelImage>::import( fm );
352  trace.endBlock();
353  Calculus::PrimalForm0 m( AT.calculus );
354  for ( Calculus::Index index = 0; index < m.myContainer.rows(); index++)
355  {
356  auto cell = m.getSCell( index );
357  double col = ((double) mask( K.sCoords( cell ) )) / 255.0;
358  m.myContainer( index ) = col > 0.0 ? 1.0 : 0.0;
359  }
360  AT.setAlpha( a, m );
361  if ( grey_image )
362  {
363  ostringstream ossGM;
364  ossGM << boost::format("%s-g-mask.pgm") %f2;
365  GreyLevelImage image_mg( domain );
366  Calculus::DualForm2 mg = AT.primal_h0 * functions::dec::diagonal( m ) * AT.getG( 0 );
367  functions::dec::form2ToGreyLevelImage
368  ( AT.calculus, mg, image_mg, 0.0, 1.0, 1 );
369  PGMWriter<GreyLevelImage>::exportPGM( ossGM.str(), image_mg );
370  }
371  else if ( color_image )
372  {
373  ostringstream ossGM;
374  ossGM << boost::format("%s-g-mask.ppm") %f2;
375  ColorImage image_mg( domain );
376  Calculus::DualForm2 mg0 = AT.primal_h0 * functions::dec::diagonal( m ) * AT.getG( 0 );
377  Calculus::DualForm2 mg1 = AT.primal_h0 * functions::dec::diagonal( m ) * AT.getG( 1 );
378  Calculus::DualForm2 mg2 = AT.primal_h0 * functions::dec::diagonal( m ) * AT.getG( 2 );
379  functions::dec::threeForms2ToRGBColorImage
380  ( AT.calculus, mg0, mg1, mg2, image_mg, 0.0, 1.0, 1 );
381  PPMWriter<ColorImage, functors::Identity >::exportPPM( ossGM.str(), image_mg );
382  }
383  }
384  else
385  AT.setAlpha( a );
386 
387  trace.info() << AT << std::endl;
388  double n_v = 0.0;
389  double eps = 0.0;
390  while ( l1 >= l2 )
391  {
392  trace.info() << "************ lambda = " << l1 << " **************" << endl;
393  AT.setLambda( l1 );
394  for ( eps = e1; eps >= e2; eps /= er )
395  {
396  trace.info() << " ======= epsilon = " << eps << " ========" << endl;
397  AT.setEpsilon( eps );
398  int n = 0;
399  do {
400  trace.progressBar( n, nbiter );
401  AT.solveU();
402  AT.solveV();
403  AT.checkV();
404  n_v = AT.computeVariation();
405  } while ( ( n_v > 0.0001 ) && ( ++n < nbiter ) );
406  trace.progressBar( n, nbiter );
407  trace.info() << "[#### last variation = " << n_v << " " << endl;
408  }
409  if ( grey_image )
410  {
411  if ( verb > 0 ) trace.beginBlock("Writing u[0] as PGM image");
412  ostringstream ossU, ossV, ossW;
413  ossU << boost::format("%s-a%.5f-l%.7f-u.pgm") % f2 % a % l1;
414  ossV << boost::format("%s-a%.5f-l%.7f-u-v.pgm") % f2 % a % l1;
415  ossW << boost::format("%s-a%.5f-l%.7f-u-v.ppm") % f2 % a % l1;
416  Calculus::DualForm2 u = AT.primal_h0 * AT.getU( 0 );
417  Calculus::DualForm1 v = AT.primal_h1 * AT.getV();
418  // Restored image
419  GreyLevelImage image_u( domain );
420  functions::dec::form2ToGreyLevelImage
421  ( AT.calculus, u, image_u, 0.0, 1.0, 1 );
422  PGMWriter<GreyLevelImage>::exportPGM( ossU.str(), image_u );
423  // Zoomed restored image with discontinuities (in black).
424  GreyLevelImage image_uv( out_domain );
425  functions::dec::form2ToGreyLevelImage
426  ( AT.calculus, u, image_uv, 0.0, 1.0, pix_sz );
427  functions::dec::dualForm1ToGreyLevelImage
428  ( AT.calculus, v, image_uv, 0.0, 1.0, pix_sz );
429  PGMWriter<GreyLevelImage>::exportPGM( ossV.str(), image_uv );
430  // Zoomed restored image with discontinuities (in specified color).
431  ColorImage cimage( out_domain );
432  functions::dec::threeForms2ToRGBColorImage
433  ( AT.calculus, u, u, u, cimage, 0.0, 1.0, pix_sz );
434  functions::dec::dualForm1ToRGBColorImage
435  ( AT.calculus, v, cimage, color_v, 0.0, 1.0, pix_sz );
436  PPMWriter<ColorImage, functors::Identity >::exportPPM( ossW.str(), cimage );
437  if ( verb > 0 ) trace.endBlock();
438  }
439  else if ( color_image )
440  {
441  if ( verb > 0 ) trace.beginBlock("Writing u[0,1,2] as PGM image");
442  ostringstream ossU, ossV;
443  ossU << boost::format("%s-a%.5f-l%.7f-u.ppm") % f2 % a % l1;
444  ossV << boost::format("%s-a%.5f-l%.7f-u-v.ppm") % f2 % a % l1;
445  Calculus::DualForm2 u0 = AT.primal_h0 * AT.getU( 0 );
446  Calculus::DualForm2 u1 = AT.primal_h0 * AT.getU( 1 );
447  Calculus::DualForm2 u2 = AT.primal_h0 * AT.getU( 2 );
448  Calculus::DualForm1 v = AT.primal_h1 * AT.getV();
449  // Restored image
450  ColorImage image_u( domain );
451  functions::dec::threeForms2ToRGBColorImage
452  ( AT.calculus, u0, u1, u2, image_u, 0.0, 1.0, 1 );
453  PPMWriter<ColorImage, functors::Identity >::exportPPM( ossU.str(), image_u );
454  ColorImage image_uv( out_domain );
455  functions::dec::threeForms2ToRGBColorImage
456  ( AT.calculus, u0, u1, u2, image_uv, 0.0, 1.0, pix_sz );
457  functions::dec::dualForm1ToRGBColorImage
458  ( AT.calculus, v, image_uv, color_v, 0.0, 1.0, pix_sz );
459  PPMWriter<ColorImage, functors::Identity >::exportPPM( ossV.str(), image_uv );
460  if ( verb > 0 ) trace.endBlock();
461  }
462  // Compute SNR if possible
463  if ( snr )
464  {
465  double u_snr = AT.computeSNR();
466  trace.info() << "- SNR of u = " << u_snr << " SNR of g = " << g_snr << endl;
467  }
468  l1 /= lr;
469  }
470  return 0;
471 }
main
int main(int argc, char **argv)
Definition: 3dHeightMapViewer.cpp:206
DGtal::Color::green
void green(const unsigned char aGreenValue)
DGtal::Color::red
void red(const unsigned char aRedValue)
DGtal::Color::blue
void blue(const unsigned char aBlueValue)
DGtal::DiscreteExteriorCalculus::Index
LinearAlgebraBackend::DenseVector::Index Index
Domain
typename Self::Domain Domain
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
DGtal::KhalimskySpaceND::sCoords
Point sCoords(const SCell &c) const
DGtal::Trace::error
std::ostream & error()
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::progressBar
void progressBar(const double currentValue, const double maximalValue)
DGtal::Trace::endBlock
double endBlock()
DGtal::trace
Trace trace(traceWriterTerm)
DGtal::ATu0v1
Aim: This class solves Ambrosio-Tortorelli functional in a plane for u a (vector of) 0-form(s) and v ...
Definition: ATu0v1.h:75
Companion project associated to DGtal ; Generated on Mon Aug 26 2024 21:31:30 for DGtalTools by   doxygen 1.9.1