im.hh

   1 #ifndef IM_HH
   2 #define IM_HH
   3
   4 #include <assert.h>
   5 #include <climits>
   6 #include <CImg.h>
   7 #include "util.hh"
   8
   9 using namespace cimg_library;
  10
  11 /** General image processing tools. */
  12 namespace im {
  13
  14   /** A struct for computing the sum of pixels in an image. */
  15   template <typename T>
  16   struct PixelSum {
  17
  18     /** Initialize the summer with a reference to an image. */
  19     PixelSum(const CImg<T> &img) : img(img), sum(0), count(0) { }
  20
  21     /** Process a pixel by summing. */
  22     void operator()(int x, int y)
  23     {
  24       sum += img(x, y);
  25       count++;
  26     }
  27
  28     const CImg<T> &img; //!< The image from which the pixel values are fetched
  29     T sum; //!< The sum of the pixels
  30     int count; //!< The number of pixels summed
  31   };
  32
  33
  34   /**
  35       Create a rectangular zero-sum peak filter.
  36
  37       The filter is first filled with the value of \c -sign, and the
  38       center value of the filter is set so that the sum is zero.
  39
  40       \param width = the width of the filter (must be odd)
  41       \param sign = the sign of the peak
  42       \return the filter image
  43   */
  44   template <typename T>
  45   CImg<T>
  46   peak_filter(int width, int sign = 1)
  47   {
  48     assert(width % 2 == 1);
  49     CImg<T> filter(width, width);
  50     filter.fill(-sign);
  51     filter(filter.width / 2, filter.height / 2) -= filter.sum();
  52
  53     return filter;
  54   }
  55
  56   /** Set all negative values of the image to zero.
  57       \param img = the image to be processed
  58   */
  59   template <typename T>
  60   void
  61   zero_negatives(CImg<T> &img)
  62   {
  63     for (int i = 0; i < (int)img.size(); i++)
  64       if (img[i] < 0)
  65         img[i] = 0;
  66   }
  67
  68   /** Weigh an image with a centered gaussian.
  69       \param img = the image to weight
  70       \param xc = center-x of the gaussian
  71       \param yc = center-y of the gaussian
  72       \param x_sigma2 = the variance in x-dimension
  73       \param y_sigma2 = the variance in y-dimension
  74   */
  75   template <typename T>
  76   void
  77   weight_gaussian(CImg<T> &img, float xc, float yc,
  78                   float x_sigma2, float y_sigma2)
  79   {
  80     cimg_mapXY(img, x, y) {
  81       float exponent = 0;
  82       if (x_sigma2 > 0)
  83         exponent += -0.5 * (xc - x) * (xc - x) / x_sigma2;
  84       if (y_sigma2 > 0)
  85         exponent += -0.5 * (yc - y) * (yc - y) / y_sigma2;
  86       img(x,y) *= exp(exponent);
  87     }
  88   }
  89
  90   /**
  91       Compute the Hough transform.
  92
  93       Each point in the resulting image correspond to a straight line
  94       in the original image.  The x-axis represents the angle (theta)
  95       of the normal of the line, and the y-axis represents the
  96       distance (rho) between the line and the center of the image.
  97       The width of the result will be \c num_thetas, and the height of
  98       the result will be (2 * \c max_rho + 1).  The center row of the
  99       Hough image corresponds to lines crossing the center of the
 100       original image.
 101
 102       \param src = the image to transform
 103       \param theta1 = the smallest value of theta
 104       \param theta2 = the largest value of theta
 105       \param num_thetas = the width of the resulting image
 106       \param max_rho = the maximum distance considered
 107       \return Hough transform of \c src
 108   */
 109   template<typename T>
 110   CImg<T> hough(const CImg<T> &src, float theta1, float theta2,
 111                 int num_thetas, int max_rho)
 112   {
 113     int num_rhos = max_rho * 2 + 1;
 114     int rho_center = num_rhos / 2;
 115
 116     CImg<T> result = CImg<T>(num_thetas, num_rhos);
 117     result.fill(0);
 118
 119     cimg_mapXY(src, x, y) {
 120
 121       // Skip non-positive pixels
 122       if (src(x,y) <= 0)
 123         continue;
 124
 125       float rho_x = (float)x - src.width/2;
 126       float rho_y = (float)y - src.height/2;
 127       float rho0 = std::sqrt(rho_x * rho_x + rho_y * rho_y);
 128       float theta0 = std::atan2(rho_y, rho_x) / M_PI * 180;
 129
 130       if (theta0 < 0) {
 131         theta0 += 180;
 132         rho0 = - rho0;
 133       }
 134
 135       // Iterate all thetas and compute corresponding rho
 136       float theta_delta = (theta2 - theta1) / (num_thetas - 1);
 137       for (int t = 0; t < num_thetas; t++) {
 138         float theta = theta1 + t * theta_delta;
 139         float r = rho0 * std::cos((theta0 - theta) / 180 * M_PI) + rho_center;
 140         if (r < 0 || r >= result.dimy())
 141           continue;
 142
 143         int ir = (int)floor(r);
 144         float w = r - ir;
 145
 146         result(t, ir) += src(x, y) * (1 - w);
 147         if (ir + 1 < (int)result.height)
 148           result(t, ir+1) += src(x, y) * w;
 149       }
 150     }
 151
 152     return result;
 153   }
 154
 155
 156   /** Median of the values in a one-dimensional row image.
 157    *
 158    * \note It is safe to specify ranges outside the image.  The range
 159    * will be clipped.
 160    *
 161    * \param img = the source image
 162    * \param x1 = the start of the range
 163    * \param x2 = the end of the range
 164    */
 165   template <typename T>
 166   T
 167   median(CImg<T> &img, int x1, int x2)
 168   {
 169     if (x1 < 0)
 170       x1 = 0;
 171     if (x2 >= (int)img.width)
 172       x2 = img.width - 1;
 173
 174     std::vector<T> values(x2 - x1 + 1);
 175     for (int x = x1; x <= x2; x++)
 176       values[x - x1] = img(x);
 177     return util::median(values);
 178   }
 179
 180
 181   /** Remove a peak from a one-dimensional image.
 182    *
 183    * \param img = the image to process
 184    * \param x0 = the location of the peak
 185    * \param win = the width of the median analysis window
 186    * \param value = the value to set the peak neighbourhood to
 187    */
 188   template <typename T>
 189   void
 190   median_peak_remove(CImg<T> &img, int x0, int win, T value = 0)
 191   {
 192     // Compute the threshold to stop the removing
 193     T threshold = im::median(img, x0 - win, x0 + win);
 194
 195     // Remove values in the right neighbourhood
 196     for (int i = 0; i < win; i++) {
 197       int x = x0 + i;
 198       if (x >= (int)img.width)
 199         break;
 200       if (img(x) < threshold)
 201         break;
 202       img(x) = value;
 203     }
 204
 205     // Remove values in the left neighbourhood
 206     for (int i = 1; i < win; i++) {
 207       int x = x0 - i;
 208       if (x < 0)
 209         break;
 210       if (img(x) < threshold)
 211         break;
 212       img(x) = value;
 213     }
 214   }
 215
 216   /** Find a maximum from a one-dimensional image.
 217    * \param img = the image to find maximum from
 218    * \param x1 = the start of the range (default: 0)
 219    * \param x2 = the end of the range (default: the last pixel of the image)
 220    * \return the position of the maximum
 221    */
 222   template <typename T>
 223   int
 224   find_max1(const CImg<T> &img, int x1 = 0, int x2 = -1)
 225   {
 226     if (x2 < 0)
 227       x2 = img.width - 1;
 228
 229     int best_x = x1;
 230     T best = img(best_x);
 231
 232     while (x1 < x2) {
 233       x1++;
 234       if (img(x1) > best) {
 235         best = img(x1);
 236         best_x = x1;
 237       }
 238     }
 239
 240     return best_x;
 241   }
 242
 243   /** Find maximum along a row of a two-dimensional image.
 244    * \param img = the image to search the maximum from
 245    * \param y = the row to analyse
 246    * \param x1 = the start of the range (default: 0)
 247    * \param x2 = the end of the range (default: the last pixel of the row)
 248    * \return the location of the maximum
 249    */
 250   template <typename T>
 251   int
 252   find_max_x(const CImg<T> &img, int y, int x1 = 0, int x2 = -1)
 253   {
 254     if (x2 < 0)
 255       x2 = img.width - 1;
 256
 257     int best_x = x1;
 258     T best = img(best_x, y);
 259
 260     while (x1 < x2) {
 261       x1++;
 262       if (img(x1, y) > best) {
 263         best = img(x1, y);
 264         best_x = x1;
 265       }
 266     }
 267
 268     return best_x;
 269   }
 270
 271   /** Paste an image to another image.
 272    * \param src = the image to be pasted
 273    * \param tgt = the target images
 274    * \param x0 = target location
 275    * \param y0 = target location
 276    * \return reference to the target image
 277    */
 278   template <typename T>
 279   T&
 280   paste_image(const T &src, T &tgt, int x0, int y0)
 281   {
 282     for (int x = 0; x < (int)src.width; x++) {
 283       for (int y = 0; y < (int)src.height; y++) {
 284         tgt(x0 + x, y0 + y) = src(x, y);
 285       }
 286     }
 287     return tgt;
 288   }
 289
 290   /** Compute the sum of each column of the image.
 291    * \param img = the source image
 292    * \return a one-dimensional image containing the sum of each column
 293    */
 294   template<typename T>
 295   CImg<T> sum_y(const CImg<T> &img)
 296   {
 297     CImg<T> result(img.width, 1);
 298     for (int x = 0; x < (int)img.width; x++) {
 299       T sum = 0;
 300       for (int y = 0; y < (int)img.height; y++)
 301         sum += img(x,y);
 302       result(x) = sum;
 303     }
 304     return result;
 305   }
 306
 307   /** Set range of values in an one-dimensional image.
 308    *
 309    * \note It is safe to specify ranges outside the image.  The range
 310    * will be clipped appropriately.
 311    *
 312    * \param img = the target image
 313    * \param x1 = the start of the range
 314    * \param x2 = the end of the range
 315    * \param value = the value to set
 316    * \return a reference to the target image
 317    */
 318   template <typename T>
 319   CImg<T>&
 320   set_range1(CImg<T> &img, int x1 = 0, int x2 = INT_MAX, T value = 0)
 321   {
 322     if (x1 < 0)
 323       x1 = 0;
 324     if (x2 > (int)img.width - 1)
 325       x2 = img.width - 1;
 326
 327     for (int x = x1; x <= x2; x++)
 328       img(x) = value;
 329
 330     return img;
 331   }
 332
 333   /** Compute the maximum of each row.
 334    *
 335    * \warning The validity of the range is not checked.
 336    *
 337    * \param img = the source image
 338    * \param x1 = the start of the range
 339    * \param x2 = the end of the range (negative: use width of image)
 340    * \return a one dimensional row image containing the maximums
 341    */
 342   template <typename T>
 343   CImg<T>
 344   max_x(const CImg<T> &img, int x1 = 0, int x2 = -1)
 345   {
 346     CImg<T> result(img.height);
 347
 348     if (x2 < 0)
 349       x2 = img.width - 1;
 350
 351     for (int y = 0; y < (int)img.height; y++) {
 352       T max = img(x1, y);
 353       for (int x = x1 + 1; x <= x2; x++) {
 354         if (img(x, y) > max)
 355           max = img(x, y);
 356       }
 357       result(y) = max;
 358     }
 359
 360     return result;
 361   }
 362
 363   /** Compute the sum of the pixels along a line
 364    *
 365    * \bug The float coordinates should be handled more elegantly so
 366    * that the line is processed along the longer axis and the other
 367    * coordinate is computed for each middle location.
 368    *
 369    * \warning Does not check the image boundaries.
 370    * \param img = the source image
 371    * \param line = the line to sum
 372    * \return The sum of the pixel values along the line.
 373    */
 374   template <typename T>
 375   T
 376   line_sum(CImg<T> &img, const geom::Line &line)
 377   {
 378     // Check single point
 379     if ((int)(line.a.x + 0.5) == (int)(line.b.x + 0.5) &&
 380         (int)(line.a.y + 0.5) == (int)(line.b.y + 0.5))
 381       return img((int)(line.a.x + 0.5), (int)(line.b.y + 0.5));
 382
 383     // Compute line
 384     float dx = line.b.x - line.a.x;
 385     float dy = line.b.y - line.a.y;
 386     float div = cimg::max(cimg::abs(dx), cimg::abs(dy));
 387     dx /= div;
 388     dy /= div;
 389
 390     // Add value to points along the line
 391     T result = 0;
 392     for (int i = 0; i <= div; i++) {
 393       int x = (int)(line.a.x + i * dx + 0.5);
 394       int y = (int)(line.a.y + i * dy + 0.5);
 395       result += img(x, y);
 396     }
 397
 398     return result;
 399   }
 400
 401 };
 402
 403 #endif /* IM_HH */