smooth.cc

   1
   2 #define NB_MAX_PICS 128
   3
   4 #include <iostream.h>
   5 #include <fstream.h>
   6
   7 extern "C"
   8 {
   9 #include "math.h"
  10 #include "string.h"
  11 }
  12
  13 #include "pictures.h"
  14
  15 int main(int argc, char **argv)
  16 {
  17   Picture *pictures[NB_MAX_PICS], *result, *ref_pic;
  18   char result_name[64], ref_name[64];
  19   int nb_pics, x, y, d;
  20   float av, av0, te, te0, sum, sum_sq, v;
  21   float *array;
  22   int k, width, height;
  23
  24   strcpy(result_name, "./result.ppm");
  25   strcpy(ref_name, "");
  26   nb_pics = 0;
  27
  28   while(argc > 1)
  29     {
  30       if(strcmp(argv[1], "-p")==0)
  31         {
  32           argv += 1; argc -= 1;
  33           while(argc > 1)
  34             {
  35               if(nb_pics == NB_MAX_PICS)
  36                  {
  37                    cerr<<"Too many pics.\n";
  38                    exit(1);
  39                  }
  40               pictures[nb_pics] = new Picture;
  41               cerr<<"Loading "<<argv[1]<<" ... "; cerr.flush();
  42               pictures[nb_pics]->LoadPpm(argv[1]);
  43               cerr<<"done.\n";
  44               if(nb_pics==0)
  45                 {
  46                   width = pictures[nb_pics]->SizeX;
  47                   height = pictures[nb_pics]->SizeY;
  48                 }
  49               else
  50                 {
  51                   if((width != pictures[nb_pics]->SizeX) ||
  52                      (height != pictures[nb_pics]->SizeY))
  53                     {
  54                       cerr<<"Size mismatch.\n";
  55                       exit(1);
  56                     }
  57                 }
  58               nb_pics++;
  59               argv += 1; argc -= 1;
  60             }
  61         }
  62       else if(strcmp(argv[1], "-o")==0)
  63         {
  64           if(argc>2)
  65             {
  66               strncpy(result_name, argv[2], 64);
  67               argv += 2; argc -= 2;
  68               cout<<"Output sent to "<<result_name<<".\n";
  69             }
  70           else
  71             {
  72               cerr<<"No filename.\n";
  73               exit(1);
  74             }
  75         }
  76       else if(strcmp(argv[1], "-r")==0)
  77         {
  78           if(argc>2)
  79             {
  80               strncpy(ref_name, argv[2], 64);
  81               argv += 2; argc -= 2;
  82               cout<<"Using "<<ref_name<<" as reference picture.\n";
  83             }
  84           else
  85             {
  86               cerr<<"No filename.\n";
  87               exit(1);
  88             }
  89         }
  90       else
  91         {
  92           cerr<<"Unknown option "<<argv[1]<<".\n";
  93           exit(1);
  94         }
  95     }
  96
  97   cerr<<"Width="<<width<<" Height="<<height<<"\n";
  98   cerr<<"Smoothing ... "; cerr.flush();
  99
 100   array = new float[nb_pics];
 101
 102   result = new Picture(width, height);
 103
 104   for(y=0; y<height; y++) for(x=0; x<width; x++)
 105     {
 106       sum = 0;
 107       for(k=0; k<nb_pics; k++)
 108         sum += pictures[k]->GetComponent(x, y, PIC_RED);
 109       for(d=0; d<PIC_DEPTH; d++)
 110         result->SetComponent(x, y, d, sum/float(nb_pics));
 111     }
 112
 113   cerr<<"done.\n";
 114
 115   if(strlen(ref_name) > 0)
 116     {
 117       cerr<<"Lightning correction ... "; cerr.flush();
 118       ref_pic = new Picture;
 119       ref_pic->LoadPpm(ref_name);
 120
 121       sum = 0; sum_sq = 0;
 122       for(x=0; x<width; x++) for(y=0; y<height; y++)
 123         {
 124           v = ref_pic->GetComponent(x, y, PIC_RED);
 125           sum += v;
 126           sum_sq += v*v;
 127         }
 128       av0 = sum/float(width * height);
 129       te0 = sqrt(sum_sq/float(width * height) - av0*av0);
 130
 131       for(x=0; x<width; x++)
 132         {
 133           sum = 0; sum_sq = 0;
 134           for(y=0; y<height; y++)
 135             {
 136               v = ref_pic->GetComponent(x, y, PIC_RED);
 137               sum += v;
 138               sum_sq += v*v;
 139             }
 140
 141           av = sum/float(height);
 142           te = sqrt(sum_sq/float(height) - av*av);
 143
 144           for(y=0; y<height; y++)
 145             {
 146               v = result->GetComponent(x, y, PIC_RED);
 147               //if(te>0) v = av0 + (v - av)*(te0/te);
 148               //else
 149               v = av0 + (v - av);
 150               if(v < 0) v = 0; else if(v > 1) v = 1;
 151               for(d=0; d<PIC_DEPTH; d++) result->SetComponent(x, y, d, v);
 152             }
 153         }
 154       delete ref_pic;
 155       cerr<<"done.\n";
 156     }
 157   else cerr<<"No column-lightning correction.\n";
 158
 159   for(k=0; k<nb_pics; k++) delete pictures[k];
 160
 161   cerr<<"Saving ... "; cerr.flush();
 162   result->SavePpm(result_name);
 163   cerr<<"done.\n";
 164   delete result;
 165 }