libcrn  3.9.5
A document image processing library
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRNImage.hpp
Go to the documentation of this file.
1 /* Copyright 2007-2015 Yann LEYDIER, INSA-Lyon, CoReNum, Université Paris Descartes
2  *
3  * This file is part of libCRN.
4  *
5  * libCRN is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU Lesser General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * libCRN is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public License
16  * along with libCRN. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * file: CRNImage.hpp
19  * \author Yann LEYDIER
20  */
21 
22 #ifndef CRNIMAGEHPP
23 #define CRNIMAGEHPP
24 
25 #include <iostream>
26 
27 #include <CRNGeometry/CRNRect.h>
28 #include <CRNMath/CRNMatrixInt.h>
32 #include <set>
33 
34 #ifdef _MSC_VER // remove annoying warnings in Microsoft Visual C++
35 # pragma warning(disable:4800)
36 #endif
37 
38 namespace crn
39 {
40  /**************************************************************************************
41  * Construction and copy
42  *************************************************************************************/
43 
44  /****************************************************************************/
55  template<typename T> Image<T>::Image(size_t w, size_t h, pixel_type val):
56  ImageBase{w, h},
57  pixels(w * h, val)
58  { }
59 
60  /****************************************************************************/
71  template<typename T> Image<T>::Image(size_t w, size_t h, const pixel_type *data):
72  ImageBase{w, h},
73  pixels(data, data + w * h)
74  { }
75 
79  template<typename T> template<typename Y> Image<T>::Image(const Image<Y> &img):
80  ImageBase(img.GetWidth(), img.GetHeight()),
81  pixels(img.begin(), img.end())
82  { }
83 
87  template<typename T> Image<T>::Image(const Image<typename BoolNotBool<T>::type> &img):
88  ImageBase(img.GetWidth(), img.GetHeight()),
89  pixels(img.GetWidth() * img.GetHeight())
90  {
91  for (auto tmp : Range(img))
92  pixels[tmp] = img.At(tmp) ? std::numeric_limits<T>::max() : T(0);
93  }
94 
99  template<typename T> template<typename Y> Image<T>::Image(const Image<Y> &img, const Rect &bbox):
100  ImageBase(bbox.GetWidth(), bbox.GetHeight()),
101  pixels(bbox.GetWidth() * bbox.GetHeight())
102  {
103  for (int y = 0; y < height; ++y)
104  std::copy_n(img.begin() + bbox.GetLeft() + (y + bbox.GetTop()) * img.GetWidth(), width, pixels.begin() + y * width);
105  }
106 
111  template<typename T> Image<T>::Image(const Image<typename BoolNotBool<T>::type> &img, const Rect &bbox):
112  ImageBase(bbox.GetWidth(), bbox.GetHeight()),
113  pixels(bbox.GetWidth() * bbox.GetHeight())
114  {
115  FOREACHPIXEL(x, y, *this)
116  At(x, y) = img.At(x + bbox.GetLeft(), y + bbox.GetTop()) ? std::numeric_limits<T>::max() : T(0);
117  }
118 
123  template<typename T> template<typename Y> Image<T>& Image<T>::operator=(const Image<Y> &img)
124  {
125  width = img.GetWidth();
126  height = img.GetHeight();
127  pixels.resize(width * height);
128  std::copy(img.begin(), img.end(), pixels.begin());
129  }
130 
135  template<typename T> Image<T>& Image<T>::operator=(const Image<typename BoolNotBool<T>::type> &img)
136  {
137  width = img.GetWidth();
138  height = img.GetHeight();
139  pixels.resize(width * height);
140  for (auto tmp : Range(img))
141  pixels[tmp] = img.At(tmp) ? std::numeric_limits<T>::max() : T(0);
142  }
143 
147  template<typename T> template<typename Y> void Image<T>::Assign(const Image<Y> &img)
148  {
149  width = img.GetWidth();
150  height = img.GetHeight();
151  pixels.resize(width * height);
152  for (auto tmp : Range(img))
153  pixels[tmp] = T(img.At(tmp));
154  }
155 
160  template<typename T> void Image<T>::SavePNG(const Path &fname) const
161  {
162  impl::SavePNG(*this, fname);
163  }
169  template<typename T> void Image<T>::SaveJPEG(const Path &fname, unsigned int qual) const
170  {
171  impl::SaveJPEG(*this, fname, Min(qual, 100u));
172  }
173 
178  template<typename T> bool Image<T>::operator==(const Image &other) const
179  {
180  if ((GetWidth() != other.GetWidth()) || (GetHeight() != other.GetHeight()))
181  return false;
182  return std::equal(begin(), end(), other.begin());
183  }
184 
185  /**************************************************************************************
186  * Arithmetic
187  *************************************************************************************/
188 
194  template<typename T> Image<T>& Image<T>::operator+=(const Image<T> &img)
195  {
196  if ((GetWidth() != img.GetWidth()) || (GetHeight() != img.GetHeight()))
197  throw ExceptionDimension("Image+=(Image): images have different sizes.");
198  for (auto i : Range(img))
199  At(i) = T(At(i) + img.At(i));
200  return *this;
201  }
202 
208  template<typename T> Image<T>& Image<T>::operator-=(const Image<T> &img)
209  {
210  if ((GetWidth() != img.GetWidth()) || (GetHeight() != img.GetHeight()))
211  throw ExceptionDimension("Image-=(Image): images have different sizes.");
212  for (auto i : Range(img))
213  At(i) = T(At(i) - img.At(i));
214  return *this;
215  }
216 
221  template<typename T> Image<T>& Image<T>::operator*=(double f)
222  {
223  for (auto &v : pixels)
224  v = T(DecimalType<T>(v) * f);
225  return *this;
226  }
227 
233  template<typename T> Image<T>& Image<T>::operator*=(const Image<T> &img)
234  {
235  if ((GetWidth() != img.GetWidth()) || (GetHeight() != img.GetHeight()))
236  throw ExceptionDimension("Image*=(Image): images have different sizes.");
237  for (auto i : Range(img))
238  At(i) = T(At(i) * img.At(i));
239  return *this;
240  }
241 
247  template<typename T> Image<T>& Image<T>::operator/=(const Image<T> &img)
248  {
249  if ((GetWidth() != img.GetWidth()) || (GetHeight() != img.GetHeight()))
250  throw ExceptionDimension("Image/=(Image): images have different sizes.");
251  for (auto i : Range(img))
252  At(i) = T(At(i) / img.At(i));
253  return *this;
254  }
255 
256  /**************************************************************************************
257  * Edition
258  *************************************************************************************/
259 
261  template<typename T> void Image<T>::Negative()
262  {
263  for (size_t tmp = 0; tmp < width * height; ++tmp)
264  pixels[tmp] = -pixels[tmp];
265  }
266 
270  template<typename T> void Image<T>::Complement(T maxval)
271  {
272  for (size_t tmp = 0; tmp < width * height; ++tmp)
273  pixels[tmp] = T(maxval - pixels[tmp]);
274  }
275 
282  template<typename T> template<typename Y> void Image<T>::Blit(const Image<Y> &src, const Rect &srczone, size_t dx, size_t dy)
283  {
284  if ((dx >= width) || (dy >= height)) // nothing to be done
285  return;
286  // source clipping
287  auto bbox = srczone & src.GetBBox();
288  // destination clipping
289  if (dx + bbox.GetWidth() > width)
290  bbox.SetWidth(int(width) - int(dx));
291  if (dy + bbox.GetHeight() > height)
292  bbox.SetWidth(int(height) - int(dy));
293  // copy
294  for (int y = 0; y < bbox.GetHeight(); ++y)
295  std::copy_n(src.begin() + bbox.GetLeft() + (y + bbox.GetTop()) * src.GetWidth(), bbox.GetWidth(), pixels.begin() + dx + (dy + y) * width);
296  }
297 
298  /****************************************************************************/
311  template<typename T> void Image<T>::FloodFill(size_t x, size_t y, const pixel_type &val, crn::DistanceType dist)
312  {
313  if ((x >= width) || (y >= height))
314  throw crn::ExceptionDomain("Image::FloodFill(): Coordinates out of bounds.");
315  size_t offset = x + y * width;
316  pixel_type oldval(pixels[offset]);
317  if (oldval == val)
318  return; // do not fill twice
319  pixels[offset] = val;
320  if ((dist == crn::DistanceType::D4) || (dist == crn::DistanceType::D8))
321  {
322  if (x > 0)
323  if (pixels[offset - 1] == oldval)
324  FloodFill(x - 1, y, val, dist);
325  if (x < width - 1)
326  if (pixels[offset + 1] == oldval)
327  FloodFill(x + 1, y, val, dist);
328  if (y > 0)
329  if (pixels[offset - width] == oldval)
330  FloodFill(x, y - 1, val, dist);
331  if (y < height - 1)
332  if (pixels[offset + width] == oldval)
333  FloodFill(x, y + 1, val, dist);
334  }
335  if (dist == crn::DistanceType::D8)
336  {
337  if ((x > 0) && (y > 0))
338  if (pixels[offset - 1 - width] == oldval)
339  FloodFill(x - 1, y - 1, val, dist);
340  if ((x > 0) && (y < height - 1))
341  if (pixels[offset - 1 + width] == oldval)
342  FloodFill(x - 1, y + 1, val, dist);
343  if ((x < width - 1) && (y > 0))
344  if (pixels[offset + 1 - width] == oldval)
345  FloodFill(x + 1, y - 1, val, dist);
346  if ((x < width - 1) && (y < height - 1))
347  if (pixels[offset + 1 + width] == oldval)
348  FloodFill(x + 1, y + 1, val, dist);
349  }
350  }
351 
352  /****************************************************************************/
363  template<typename T> void Image<T>::ScanFill(size_t x, size_t y, const pixel_type &val, crn::DistanceType dist)
364  {
365  if ((x >= width) || (y >= height))
366  throw crn::ExceptionDomain("Image::ScanFill(): Coordinates out of bounds.");
367  pixel_type oldval = pixels[x + y * width];
368  if (oldval == val)
369  return; // do not fill twice
370  auto todo = std::vector<std::pair<int, std::pair<int, int>>>{};
371  int bx;
372  for (bx = int(x); bx >= 0; --bx)
373  if (pixels[bx + y * width] != oldval)
374  break;
375  bx += 1;
376  int ex;
377  for (ex = int(x); ex < int(width); ++ex)
378  if (pixels[ex + y * width] != oldval)
379  break;
380  ex -= 1;
381  todo.emplace_back(int(y), std::make_pair(bx, ex));
382  while (!todo.empty())
383  {
384  auto pos(todo.back());
385  todo.pop_back();
386  auto checklimitup = pos.second.first - 1; // these two variables are used to avoid
387  auto checklimitdown = pos.second.first - 1; // adding the same stream more than once
388  for (int tx = pos.second.first; tx <= pos.second.second; ++tx)
389  {
390  // fill
391  pixels[tx + pos.first * width] = val;
392  // look up
393  if (pos.first > 0)
394  {
395  if ((tx > checklimitup) && (pixels[tx + (pos.first - 1) * width] == oldval))
396  {
397  // add to queue
398  for (bx = tx; bx >= 0; --bx)
399  if (pixels[bx + (pos.first - 1) * width] != oldval)
400  break;
401  bx += 1;
402  for (ex = tx; ex < width; ++ex)
403  if (pixels[ex + (pos.first - 1) * width] != oldval)
404  break;
405  ex -= 1;
406  todo.emplace_back(pos.first - 1, std::make_pair(bx, ex));
407  checklimitup = ex;
408  }
409  if (dist == crn::DistanceType::D8)
410  {
411  // look up-left (only for the first pixel of the stream)
412  if ((tx == pos.second.first) && (tx > 0) && (tx - 1 >= checklimitup))
413  {
414  if (pixels[tx - 1 + (pos.first - 1) * width] == oldval)
415  {
416  // add to queue
417  for (bx = tx - 1; bx >= 0; --bx)
418  if (pixels[bx + (pos.first - 1) * width] != oldval)
419  break;
420  bx += 1;
421  for (ex = tx - 1; ex < width; ++ex)
422  if (pixels[ex + (pos.first - 1) * width] != oldval)
423  break;
424  ex -= 1;
425  todo.emplace_back(pos.first - 1, std::make_pair(bx, ex));
426  checklimitup = ex;
427  }
428  }
429  // look up-right (only for the last pixel of the stream)
430  if ((tx == pos.second.second) && (tx < width - 1) && (tx + 1 > checklimitup))
431  {
432  if (pixels[tx + 1 + (pos.first - 1) * width] == oldval)
433  {
434  // add to queue
435  for (bx = tx + 1; bx >= 0; --bx)
436  if (pixels[bx + (pos.first - 1) * width] != oldval)
437  break;
438  bx += 1;
439  for (ex = tx + 1; ex < width; ++ex)
440  if (pixels[ex + (pos.first - 1) * width] != oldval)
441  break;
442  ex -= 1;
443  todo.emplace_back(pos.first - 1, std::make_pair(bx, ex));
444  checklimitup = ex;
445  }
446  }
447  }
448  } // look up
449  // look down
450  if ((pos.first < height - 1) && (tx > checklimitdown))
451  {
452  if (pixels[tx + (pos.first + 1) * width] == oldval)
453  {
454  // add to queue
455  for (bx = tx; bx >= 0; --bx)
456  if (pixels[bx + (pos.first + 1) * width] != oldval)
457  break;
458  bx += 1;
459  for (ex = tx; ex < width; ++ex)
460  if (pixels[ex + (pos.first + 1) * width] != oldval)
461  break;
462  ex -= 1;
463  todo.emplace_back(pos.first + 1, std::make_pair(bx, ex));
464  checklimitdown = ex;
465  }
466  if (dist == crn::DistanceType::D8)
467  {
468  // look up-left (only for the first pixel of the stream)
469  if ((tx == pos.second.first) && (tx > 0) && (tx - 1 >= checklimitdown))
470  {
471  if (pixels[tx - 1 + (pos.first + 1) * width] == oldval)
472  {
473  // add to queue
474  for (bx = tx - 1; bx >= 0; --bx)
475  if (pixels[bx + (pos.first + 1) * width] != oldval)
476  break;
477  bx += 1;
478  for (ex = tx - 1; ex < width; ++ex)
479  if (pixels[ex + (pos.first + 1) * width] != oldval)
480  break;
481  ex -= 1;
482  todo.emplace_back(pos.first + 1, std::make_pair(bx, ex));
483  checklimitdown = ex;
484  }
485  }
486  // look up-right (only for the last pixel of the stream)
487  if ((tx == pos.second.second) && (tx < width - 1) && (tx + 1 > checklimitdown))
488  {
489  if (pixels[tx + 1 + (pos.first + 1) * width] == oldval)
490  {
491  // add to queue
492  for (bx = tx + 1; bx >= 0; --bx)
493  if (pixels[bx + (pos.first + 1) * width] != oldval)
494  break;
495  bx += 1;
496  for (ex = tx + 1; ex < width; ++ex)
497  if (pixels[ex + (pos.first + 1) * width] != oldval)
498  break;
499  ex -= 1;
500  todo.emplace_back(pos.first + 1, std::make_pair(bx, ex));
501  checklimitdown = ex;
502  }
503  }
504  }
505  } // look down
506  }
507  } // while there are streams left
508  }
509 
515  template<typename T> void Image<T>::DrawRect(const Rect &r, pixel_type color, bool filled)
516  {
517  // clipping
518  const auto clip = r & GetBBox();
519  if (!clip.IsValid())
520  return;
521 
522  // top line
523  std::fill_n(pixels.begin() + clip.GetLeft() + clip.GetTop() * width, clip.GetWidth(), color);
524  // bottom line
525  std::fill_n(pixels.begin() + clip.GetLeft() + clip.GetBottom() * width, clip.GetWidth(), color);
526 
527  for (int y = clip.GetTop() + 1; y < clip.GetBottom(); ++y)
528  {
529  if (filled) // fill center
530  std::fill_n(pixels.begin() + clip.GetLeft() + y * width, clip.GetWidth(), color);
531  else // vertical lines
532  At(clip.GetLeft(), y) = At(clip.GetRight(), y) = color;
533  }
534  }
535 
543  template<typename T> void Image<T>::DrawLine(size_t x1, size_t y1, size_t x2, size_t y2, pixel_type color)
544  {
545  // compute increments
546  int lg_delta, sh_delta, cycle, lg_step, sh_step;
547  lg_delta = (int)x2 - (int)x1;
548  sh_delta = (int)y2 - (int)y1;
549  if (lg_delta < 0)
550  {
551  lg_delta = -lg_delta;
552  lg_step = -1;
553  }
554  else
555  lg_step = 1;
556 
557  if (sh_delta < 0)
558  {
559  sh_delta = -sh_delta;
560  sh_step = -1;
561  }
562  else
563  sh_step = 1;
564 
565  if (sh_delta < lg_delta)
566  {
567  cycle = lg_delta >> 1;
568  while (x1 != x2)
569  {
570  if ((x1 >= 0) && (x1 < width) && (y1 >= 0) && (y1 < height))
571  At(x1, y1) = color;
572 
573  x1 += lg_step;
574  cycle = cycle + sh_delta;
575  if (cycle > lg_delta)
576  {
577  y1 += sh_step;
578  cycle = cycle - lg_delta;
579  }
580  }
581  }
582  else
583  {
584  cycle = sh_delta >> 1;
585  int tmp = lg_delta;
586  lg_delta = sh_delta;
587  sh_delta = tmp;
588  tmp = lg_step;
589  lg_step = sh_step;
590  sh_step = tmp;
591  while (y1 != y2)
592  {
593  if ((x1 >= 0) && (x1 < width) && (y1 >= 0) && (y1 < height))
594  At(x1, y1) = color;
595 
596  y1 += lg_step;
597  cycle = cycle + sh_delta;
598  if (cycle > lg_delta)
599  {
600  x1 += sh_step;
601  cycle = cycle - lg_delta;
602  }
603  }
604  }
605  }
606 
607  /*****************************************************************************/
614  template<typename T> void Image<T>::ScaleToSize(size_t w, size_t h)
615  {
616  const auto defaultvalue = pixels.front();
617  const auto nullvalue = defaultvalue - defaultvalue;
618  // width scaling
619  auto wscale = std::vector<pixel_type>();
620  auto step = (double)width / (double)w;
621 
622  if (step == 1)
623  { // just copy (assuming the values are base types)
624  wscale.swap(pixels);
625  }
626  else
627  {
628  wscale = std::vector<pixel_type>(w * height, defaultvalue);
629  for (auto x = size_t(0); x < w; ++x)
630  {
631  auto min = int(double(x) * step);
632  auto coeffmin = 1 - (double(x) * step - double(min));
633 
634  auto max = int(double(x + 1) * step);
635  auto coeffmax = double(x + 1) * step - double(max);
636 
637  if (step < 1.0)
638  {
639  max = min + 1;
640  coeffmax = 1 - coeffmin;
641  }
642  if (max >= int(width))
643  coeffmax = 0;
644  for (size_t y = 0; y < height; ++y)
645  {
646  auto acc = DecimalType<pixel_type>(nullvalue);
647  auto coeff = 0.0;
648  auto yoff = y * width;
649 
650  acc += coeffmin * pixels[min + yoff];
651  coeff += coeffmin;
652  if (coeffmax)
653  {
654  acc += coeffmax * pixels[max + yoff];
655  coeff += coeffmax;
656  }
657  for (int k = min + 1; k < max; ++k)
658  {
659  acc += pixels[k + yoff];
660  coeff += 1;
661  }
662  // mean
663  wscale[x + y * w] = T(acc / coeff);
664  }
665  }
666  }
667 
668  // height scaling
669  step = (double)height / (double)h;
670 
671  if (step == 1)
672  { // just copy (assuming the values are base types)
673  pixels.swap(wscale);
674  }
675  else
676  {
677  pixels = std::vector<pixel_type>(w * h, defaultvalue);
678 
679  for (size_t y = 0; y < h; ++y)
680  {
681  auto min = int(double(y) * step);
682  auto coeffmin = 1 - (double(y) * step - double(min));
683 
684  auto max = int(double(y + 1) * step);
685  auto coeffmax = double(y + 1) * step - double(max);
686 
687  if (step < 1.0)
688  {
689  max = min + 1;
690  coeffmax = 1 - coeffmin;
691  }
692  if (max >= int(height))
693  coeffmax = 0;
694  for (size_t x = 0; x < w; ++x)
695  {
696  auto acc = DecimalType<pixel_type>(nullvalue);
697  auto coeff = 0.0;
698 
699  acc += coeffmin * wscale[x + min * w];
700  coeff += coeffmin;
701  if (coeffmax)
702  {
703  acc += coeffmax * wscale[x + max * w];
704  coeff += coeffmax;
705  }
706  for (int k = min + 1; k < max; ++k)
707  {
708  acc += wscale[x + k * w];
709  coeff += 1;
710  }
711 
712  // mean
713  pixels[x + y * w] = T(acc / coeff);
714  }
715  }
716  }
717 
718  width = w;
719  height = h;
720  }
721 
722  /*****************************************************************************/
729  template<typename T> void Image<T>::Flip(const Orientation &ori)
730  {
731  if (ori == Orientation::HORIZONTAL)
732  {
733  for (size_t y = 0; y < height; ++y)
734  {
735  std::reverse(pixels.begin() + y * width, pixels.begin() + (y + 1) * width);
736  }
737  }
738  else if (ori == Orientation::VERTICAL)
739  {
740  auto tmp = std::vector<T>(width);
741  for (size_t y = 0; y < height / 2; ++y)
742  {
743  std::copy_n(pixels.begin() + y * width, width, tmp.begin());
744  std::copy_n(pixels.begin() + (height - 1 - y) * width, width, pixels.begin() + y * width);
745  std::copy_n(tmp.begin(), width, pixels.begin() + (height - 1 - y) * width);
746  }
747  }
748  else
749  throw ExceptionInvalidArgument("void Image::Flip(const Orientation &ori): invalid orientation.");
750  }
751 
752  /****************************************************************************/
761  template<typename T> template<typename CMP> void Image<T>::Dilate(const MatrixInt &strel, CMP cmp)
762  {
763  if (!(strel.GetRows() & 1) || !(strel.GetCols() & 1))
764  throw ExceptionDimension("void Image::Dilate(const MatrixInt &strel): even matrix dimensions.");
765 
766  auto halfw = int(strel.GetCols()) / 2;
767  auto halfh = int(strel.GetRows()) / 2;
768 
769  auto newpix = std::vector<pixel_type>(width * height, pixels.front());
770  for (size_t y = 0; y < height; ++y)
771  {
772  for (size_t x = 0; x < width; ++x)
773  {
774  auto basey = int(y) - halfh;
775  auto basex = int(x) - halfw;
776  //auto pix = std::numeric_limits<pixel_type>::max();
777  auto pix = [&strel, basex, basey, this]()
778  {
779  for (size_t cy = 0; cy < strel.GetRows(); cy++)
780  for (size_t cx = 0; cx < strel.GetCols(); cx++)
781  {
782  if (!strel[cy][cx])
783  continue;
784  auto tx = basex + int(cx);
785  if (tx < 0) continue;
786  if (tx >= int(GetWidth())) continue;
787  auto ty = basey + int(cy);
788  if (ty < 0) continue;
789  if (ty >= int(GetHeight())) continue;
790 
791  return pixels[tx + ty * width];
792  }
793  return pixels[basex + basey * width];
794  }();
795  for (size_t cy = 0; cy < strel.GetRows(); ++cy)
796  for (size_t cx = 0; cx < strel.GetCols(); ++cx)
797  {
798  if (!strel[cy][cx])
799  continue;
800  auto tx = basex + int(cx);
801  if (tx < 0) continue;
802  if (tx >= int(GetWidth())) continue;
803  auto ty = basey + int(cy);
804  if (ty < 0) continue;
805  if (ty >= int(GetHeight())) continue;
806 
807  pixel_type tmppix = pixels[tx + ty * width];
808  if (cmp(tmppix, pix))
809  pix = tmppix;
810  }
811  newpix[x + y * width] = pix;
812  }
813  }
814  pixels.swap(newpix);
815  }
816 
817  /****************************************************************************/
827  template<typename T> template<typename CMP> void Image<T>::Erode(const MatrixInt &strel, CMP cmp)
828  {
829  if (!(strel.GetRows() & 1) || !(strel.GetCols() & 1))
830  throw ExceptionDimension("void Image::Erode(MatrixInt &strel): even matrix dimensions.");
831 
832  auto halfw = int(strel.GetCols()) / 2;
833  auto halfh = int(strel.GetRows()) / 2;
834 
835  auto newpix = std::vector<pixel_type>(width * height, pixels.front());
836  for (size_t y = 0; y < height; y++)
837  for (size_t x = 0; x < width; x++)
838  {
839  const auto basey = int(y) - halfh;
840  const auto basex = int(x) - halfw;
841  //auto pix = std::numeric_limits<pixel_type>::min();
842  auto pix = [&strel, basex, basey, this]()
843  {
844  for (size_t cy = 0; cy < strel.GetRows(); cy++)
845  for (size_t cx = 0; cx < strel.GetCols(); cx++)
846  {
847  if (!strel[cy][cx])
848  continue;
849  auto tx = basex + int(cx);
850  if (tx < 0) continue;
851  if (tx >= int(GetWidth())) continue;
852  auto ty = basey + int(cy);
853  if (ty < 0) continue;
854  if (ty >= int(GetHeight())) continue;
855 
856  return pixels[tx + ty * width];
857  }
858  return pixels[basex + basey * width];
859  }();
860  for (size_t cy = 0; cy < strel.GetRows(); cy++)
861  for (size_t cx = 0; cx < strel.GetCols(); cx++)
862  {
863  if (!strel[cy][cx])
864  continue;
865  auto tx = basex + int(cx);
866  if (tx < 0) continue;
867  if (tx >= int(GetWidth())) continue;
868  auto ty = basey + int(cy);
869  if (ty < 0) continue;
870  if (ty >= int(GetHeight())) continue;
871 
872  pixel_type tmppix = pixels[tx + ty * width];
873  if (!cmp(tmppix, pix))
874  pix = tmppix;
875  }
876  newpix[x + y * width] = pix;
877  }
878 
879  pixels.swap(newpix);
880  }
881 
890  template<typename T> template<typename CMP> void Image<T>::FastDilate(size_t halfwin, size_t index, CMP cmp)
891  {
892  if (halfwin == 0)
893  return;
894  if ((index == 0) && (halfwin < 8))
895  {
896  Dilate(MatrixInt{halfwin + halfwin + 1, halfwin + halfwin + 1, 1}, cmp);
897  return;
898  }
899  // initialize windows
900  auto wins = std::vector<std::multiset<pixel_type, CMP>>(width, std::multiset<pixel_type, CMP>(cmp));
901  for (size_t x = 0; x < width; ++x)
902  for (size_t y = 0; y < halfwin; ++y)
903  wins[x].insert(At(x, y));
904 
905  auto newpix = std::vector<pixel_type>(width * height, pixels.front());
906  for (size_t y = 0; y < height; ++y)
907  {
908  // update windows
909  if (y > halfwin)
910  {
911  const auto ry = y - halfwin - 1;
912  for (size_t x = 0; x < width; ++x)
913  {
914  auto it = wins[x].find(At(x, ry));
915  if (it != wins[x].end())
916  wins[x].erase(it);
917  }
918  }
919  if (y + halfwin < height)
920  {
921  const auto ay = y + halfwin;
922  for (size_t x = 0; x < width; ++x)
923  wins[x].insert(At(x, ay));
924  }
925  // compute mins
926  auto mins = std::vector<pixel_type>(width, pixels.front());
927  for (size_t x = 0; x < width; ++x)
928  {
929  auto cmin = wins[x].begin();
930  std::advance(cmin, Min(index, wins[x].size() - 1));
931  mins[x] = *cmin;
932  }
933 
934  // compute
935  for (size_t x = 0; x < width; ++x)
936  {
937  const auto bx = (x >= halfwin) ? x - halfwin : 0;
938  const auto ex = Min(x + halfwin, width - 1);
939  auto minval = mins[bx];
940  for (auto tx = bx + 1; tx <= ex; ++tx)
941  if (cmp(mins[tx], minval)) minval = mins[tx];
942  newpix[x + y * width] = minval;
943  } // for x
944  } // for y
945  pixels.swap(newpix);
946  }
947 
956  template<typename T> template<typename CMP> void Image<T>::FastErode(size_t halfwin, size_t index, CMP cmp)
957  {
958  if (halfwin == 0)
959  return;
960  if ((index == 0) && (halfwin < 8))
961  {
962  Erode(MatrixInt{halfwin + halfwin + 1, halfwin + halfwin + 1, 1}, cmp);
963  return;
964  }
965  // initialize windows
966  auto wins = std::vector<std::multiset<pixel_type, CMP>>(width, std::multiset<pixel_type, CMP>(cmp));
967  for (size_t x = 0; x < width; ++x)
968  for (size_t y = 0; y < halfwin; ++y)
969  wins[x].insert(At(x, y));
970 
971  auto newpix = std::vector<pixel_type>(width * height, pixels.front());
972  for (size_t y = 0; y < height; ++y)
973  {
974  // update windows
975  if (y > halfwin)
976  {
977  const auto ry = y - halfwin - 1;
978  for (size_t x = 0; x < width; ++x)
979  {
980  auto it = wins[x].find(At(x, ry));
981  if (it != wins[x].end())
982  wins[x].erase(it);
983  }
984  }
985  if (y + halfwin < height)
986  {
987  const auto ay = y + halfwin;
988  for (size_t x = 0; x < width; ++x)
989  wins[x].insert(At(x, ay));
990  }
991  // compute maxes
992  auto maxes = std::vector<pixel_type>(width, pixels.front());
993  for (size_t x = 0; x < width; ++x)
994  {
995  auto cmax = wins[x].rbegin();
996  std::advance(cmax, Min(index, wins[x].size() - 1));
997  maxes[x] = *cmax;
998  }
999 
1000  // compute
1001  for (size_t x = 0; x < width; ++x)
1002  {
1003  const auto bx = (x >= halfwin) ? x - halfwin : 0;
1004  const auto ex = Min(x + halfwin, width - 1);
1005  auto maxval = maxes[bx];
1006  for (auto tx = bx + 1; tx <= ex; ++tx)
1007  if (!cmp(maxes[tx], maxval)) maxval = maxes[tx];
1008  newpix[x + y * width] = maxval;
1009  } // for x
1010  } // for y
1011  pixels.swap(newpix);
1012  }
1013 
1014  /****************************************************************************/
1023  template<typename T> void Image<T>::Convolve(const MatrixDouble &mat)
1024  {
1025  if (!(mat.GetRows() & 1) || !(mat.GetCols() & 1))
1026  throw ExceptionDimension("void Image::Convolve(MatrixDouble &mat): even matrix dimensions.");
1027 
1028  auto halfw = mat.GetCols() / 2;
1029  auto halfh = mat.GetRows() / 2;
1030  if ((halfh > height) || (halfw > width))
1031  throw ExceptionDomain("void Image::Convolve(MatrixDouble &mat): matrix bigger than the image!");
1032 
1033  auto oldimg = pixels;
1034 
1035  const auto nullvalue = pixels.front() - pixels.front();
1036  for (size_t y = 0; y < halfh; y++)
1037  for (size_t x = 0; x < width; x++)
1038  { // top
1039  auto basey = int(y) - int(halfh);
1040  auto basex = int(x) - int(halfw);
1041  auto sum = DecimalType<pixel_type>(nullvalue);
1042  for (size_t cy = 0; cy < mat.GetRows(); cy++)
1043  for (size_t cx = 0; cx < mat.GetCols(); cx++)
1044  {
1045  auto tx = basex + int(cx);
1046  if (tx < 0) tx = 0;
1047  if (tx >= int(width)) tx = int(width) - 1;
1048  auto ty = basey + int(cy);
1049  if (ty < 0) ty = 0;
1050  if (ty >= int(height)) ty = int(height) - 1;
1051  sum += oldimg[tx + ty * width] * mat[cy][cx];
1052  }
1053  pixels[x + y * width] = pixel_type(sum);
1054  }
1055  for (int y = int(height) - int(halfh); y < int(height); y++)
1056  for (size_t x = 0; x < width; x++)
1057  { // bottom
1058  auto basey = y - int(halfh);
1059  auto basex = int(x) - int(halfw);
1060  auto sum = DecimalType<pixel_type>(nullvalue);
1061  for (size_t cy = 0; cy < mat.GetRows(); cy++)
1062  for (size_t cx = 0; cx < mat.GetCols(); cx++)
1063  {
1064  auto tx = basex + int(cx);
1065  if (tx < 0) tx = 0;
1066  if (tx >= int(width)) tx = int(width) - 1;
1067  auto ty = basey + int(cy);
1068  if (ty < 0) ty = 0;
1069  if (ty >= int(height)) ty = int(height) - 1;
1070  sum += oldimg[tx + ty * width] * mat[cy][cx];
1071  }
1072  pixels[x + y * width] = pixel_type(sum);
1073  }
1074  for (int y = int(halfh); y < int(height) - int(halfh); y++)
1075  for (size_t x = 0; x < halfw; x++)
1076  { // left
1077  auto basey = y - int(halfh);
1078  auto basex = int(x) - int(halfw);
1079  auto sum = DecimalType<pixel_type>(nullvalue);
1080  for (size_t cy = 0; cy < mat.GetRows(); cy++)
1081  for (size_t cx = 0; cx < mat.GetCols(); cx++)
1082  {
1083  auto tx = basex + int(cx);
1084  if (tx < 0) tx = 0;
1085  if (tx >= int(width)) tx = int(width) - 1;
1086  auto ty = basey + int(cy);
1087  if (ty < 0) ty = 0;
1088  if (ty >= int(height)) ty = int(height) - 1;
1089  sum += oldimg[tx + ty * width] * mat[cy][cx];
1090  }
1091  pixels[x + y * width] = pixel_type(sum);
1092  }
1093  for (int y = int(halfh); y < int(height) - int(halfh); y++)
1094  for (int x = int(width) - int(halfw); x < int(width); x++)
1095  { // right
1096  auto basey = y - int(halfh);
1097  auto basex = x - int(halfw);
1098  auto sum = DecimalType<pixel_type>(nullvalue);
1099  for (size_t cy = 0; cy < mat.GetRows(); cy++)
1100  for (size_t cx = 0; cx < mat.GetCols(); cx++)
1101  {
1102  auto tx = basex + int(cx);
1103  if (tx < 0) tx = 0;
1104  if (tx >= int(width)) tx = int(width) - 1;
1105  auto ty = basey + int(cy);
1106  if (ty < 0) ty = 0;
1107  if (ty >= int(height)) ty = int(height) - 1;
1108  sum += oldimg[tx + ty * width] * mat[cy][cx];
1109  }
1110  pixels[x + y * width] = pixel_type(sum);
1111  }
1112  for (int y = int(halfh); y < int(height) - int(halfh); y++)
1113  for (int x = int(halfw); x < int(width) - int(halfw); x++)
1114  {
1115  auto basey = y - int(halfh);
1116  auto basex = x - int(halfw);
1117  auto sum = DecimalType<pixel_type>(nullvalue);
1118  for (size_t cy = 0; cy < mat.GetRows(); cy++)
1119  for (size_t cx = 0; cx < mat.GetCols(); cx++)
1120  {
1121  sum += oldimg[basex + cx + (basey + cy) * width] * mat[cy][cx];
1122  }
1123  pixels[x + y * width] = pixel_type(sum);
1124  }
1125  }
1126 
1132  template<typename T> void Image<T>::GaussianBlur(double sigma)
1133  {
1134  // Backup... Just in case...
1135  /*
1136  auto gauss = MatrixDouble::NewGaussianLine(sigma);
1137  gauss.NormalizeForConvolution();
1138  Convolve(gauss); // may throw
1139  gauss.Transpose();
1140  Convolve(gauss); // may throw
1141  */
1142 
1143  auto kernel = MatrixDouble::NewGaussianLine(sigma);
1144 
1145  kernel.NormalizeForConvolution();
1146 
1147  const auto kernel_size(kernel.GetCols());
1148  const auto radius(kernel_size / 2);
1149  std::vector<double> coeffs(radius + 1);
1150  auto index = size_t(0);
1151 
1152  for (auto k = radius; k < kernel_size; ++k)
1153  {
1154  coeffs[index] = kernel[0][k];
1155  ++index;
1156  }
1157 
1158  /********* "Horizontal" convolution *********/
1159 
1160  index = 0;
1161  auto left_bound = size_t(0);
1162  auto right_bound = width - 1;
1163 
1164  std::vector<pixel_type> tmp_row(width, pixels.front());
1165  auto back_shifted_index = size_t(0);
1166 
1167  for (auto r = 0; r < height; ++r)
1168  {
1169  auto tmp_id = 0;
1170 
1171  // Left margin
1172 
1173  for (auto c = size_t(0); c < radius; ++c)
1174  {
1175  auto weight = coeffs[0];
1176  auto local_cumul = pixels[index] * weight;
1177 
1178  auto index_left = int(index); // Because of "--"
1179  auto index_right = index_left;
1180 
1181  for (auto offset = size_t(1); offset <= radius; ++offset)
1182  {
1183  --index_left;
1184  index_left = Max(int(left_bound), index_left);
1185  ++index_right;
1186  weight = coeffs[offset];
1187 
1188  local_cumul += (pixels[index_left] + pixels[index_right]) * weight;
1189  }
1190 
1191  tmp_row[tmp_id] = pixel_type(local_cumul);
1192  ++tmp_id;
1193  ++index;
1194  }
1195 
1196  left_bound += width;
1197 
1198  // Central region
1199 
1200  for (auto c = radius; c < width - radius; ++c)
1201  {
1202  auto weight = coeffs[0];
1203  auto local_cumul = pixels[index] * weight;
1204 
1205  auto index_left = int(index); // Because of "--"
1206  auto index_right = index_left;
1207 
1208  for (auto offset = size_t(1); offset <= radius; ++offset)
1209  {
1210  --index_left;
1211  ++index_right;
1212  weight = coeffs[offset];
1213 
1214  local_cumul += (pixels[index_left] + pixels[index_right]) * weight;
1215  }
1216 
1217  tmp_row[tmp_id] = pixel_type(local_cumul);
1218  ++tmp_id;
1219  ++index;
1220  }
1221 
1222  // Right margin
1223 
1224  for (auto c = width - radius; c < width; ++c)
1225  {
1226  auto weight = coeffs[0];
1227  auto local_cumul = pixels[index] * weight;
1228 
1229  auto index_left = int(index); // Because of "--"
1230  auto index_right = index_left;
1231 
1232  for (auto offset = 1; offset <= radius; ++offset)
1233  {
1234  --index_left;
1235  ++index_right;
1236  index_right = Min(int(right_bound), int(index_right));
1237  weight = coeffs[offset];
1238 
1239  local_cumul += (pixels[index_left] + pixels[index_right]) * weight;
1240  }
1241 
1242  tmp_row[tmp_id] = pixel_type(local_cumul);
1243  ++tmp_id;
1244  ++index;
1245  }
1246 
1247  right_bound += width;
1248 
1249  // Finalize current row
1250 
1251  for (auto k = 0; k < width; ++k)
1252  {
1253  pixels[back_shifted_index] = tmp_row[k];
1254  ++back_shifted_index;
1255  }
1256  }
1257 
1258  /********* "Vertical" convolution *********/
1259 
1260  std::vector<size_t> indices(kernel_size, 0);
1261  std::vector<pixel_type> tmp_col(height, pixels.front());
1262  auto weight = 0.0;
1263 
1264  for (auto c = 0; c < width; ++c)
1265  {
1266  index = 0;
1267 
1268  // Upper margin
1269 
1270  for (auto k = size_t(0); k <= radius; ++k)
1271  indices[k] = c;
1272 
1273  for (auto k = radius + 1; k < kernel_size; ++k)
1274  indices[k] = indices[k - 1] + width;
1275 
1276  for (auto r = size_t(0); r < radius; ++r)
1277  {
1278  weight = coeffs[0];
1279  auto local_cumul = pixels[indices[radius]] * weight;
1280 
1281  for (auto offset = 1; offset <= radius; ++offset)
1282  {
1283  auto id_upper = indices[radius - offset];
1284  auto id_lower = indices[radius + offset];
1285 
1286  weight = coeffs[offset];
1287  local_cumul += (pixels[id_upper] + pixels[id_lower]) * weight;
1288  }
1289 
1290  tmp_col[index] = pixel_type(local_cumul);
1291  ++index;
1292 
1293  for (auto k = 0; k < kernel_size - 1; ++k)
1294  indices[k] = indices[k + 1];
1295 
1296  indices[kernel_size - 1] += width;
1297 
1298  }
1299 
1300  // Central region
1301 
1302  for (auto r = radius; r < height - radius; ++r)
1303  {
1304  weight = coeffs[0];
1305  auto local_cumul = pixels[indices[radius]] * weight;
1306 
1307  for (auto offset = size_t(1); offset <= radius; ++offset)
1308  {
1309  auto id_upper = indices[radius - offset];
1310  auto id_lower = indices[radius + offset];
1311 
1312  weight = coeffs[offset];
1313  local_cumul += (pixels[id_upper] + pixels[id_lower]) * weight;
1314  }
1315 
1316  tmp_col[index] = pixel_type(local_cumul);
1317  ++index;
1318 
1319  for (auto k = 0; k < kernel_size - 1; ++k)
1320  indices[k] = indices[k + 1];
1321 
1322  indices[kernel_size - 1] += width;
1323  }
1324 
1325  // Lower margin
1326 
1327  indices[kernel_size - 1] = c + (height - 1) * width;
1328 
1329  for (auto r = height - radius; r < height; ++r)
1330  {
1331  weight = coeffs[0];
1332  auto local_cumul = pixels[indices[radius]] * weight;
1333 
1334  for (auto offset = size_t(1); offset <= radius; ++offset)
1335  {
1336  auto id_upper = indices[radius - offset];
1337  auto id_lower = indices[radius + offset];
1338 
1339  weight = coeffs[offset];
1340  local_cumul += (pixels[id_upper] + pixels[id_lower]) * weight;
1341  }
1342 
1343  tmp_col[index] = pixel_type(local_cumul);
1344  ++index;
1345 
1346  for (auto k = size_t(0); k < kernel_size - 1; ++k)
1347  indices[k] = indices[k + 1];
1348  }
1349 
1350  // Finalize current column
1351 
1352  index = c;
1353 
1354  for (auto k = size_t(0); k < height; ++k)
1355  {
1356  pixels[index] = tmp_col[k];
1357  index += width;
1358  }
1359  }
1360  }
1361 }
1362 
1364 {
1365  if ((i1.GetWidth() != i2.GetWidth()) || (i1.GetHeight() != i2.GetHeight()))
1366  throw crn::ExceptionDimension("operator+(Image, Image): images do not have the same sizes.");
1368  auto res = crn::Image<R>{i1.GetWidth(), i1.GetHeight()};
1369  for (auto i : crn::Range(i1))
1370  res.At(i) = R(i1.At(i)) + R(i2.At(i));
1371  return res;
1372 }
1373 
1375 {
1376  if ((i1.GetWidth() != i2.GetWidth()) || (i1.GetHeight() != i2.GetHeight()))
1377  throw crn::ExceptionDimension("operator-(Image, Image): images do not have the same sizes.");
1379  auto res = crn::Image<R>{i1.GetWidth(), i1.GetHeight()};
1380  for (auto i : crn::Range(i1))
1381  res.At(i) = R(i1.At(i)) - R(i2.At(i));
1382  return res;
1383 }
1384 
1386 {
1387  if ((i1.GetWidth() != i2.GetWidth()) || (i1.GetHeight() != i2.GetHeight()))
1388  throw crn::ExceptionDimension("operator*(Image, Image): images do not have the same sizes.");
1390  auto res = crn::Image<R>{i1.GetWidth(), i1.GetHeight()};
1391  for (auto i : crn::Range(i1))
1392  res.At(i) = R(i1.At(i)) * R(i2.At(i));
1393  return res;
1394 }
1395 
1397 {
1399  auto res = crn::Image<R>{i.GetWidth(), i.GetHeight()};
1400  for (auto tmp : crn::Range(i))
1401  res.At(tmp) = R(i.At(tmp)) * R(d);
1402  return res;
1403 }
1404 
1406 {
1407  return d * i;
1408 }
1409 
1411 {
1412  if ((i1.GetWidth() != i2.GetWidth()) || (i1.GetHeight() != i2.GetHeight()))
1413  throw crn::ExceptionDimension("operator/(Image, Image): images do not have the same sizes.");
1415  auto res = crn::Image<R>{i1.GetWidth(), i1.GetHeight()};
1416  for (auto i : Range(i1))
1417  res.At(i) = R(i1.At(i)) / R(i2.At(i));
1418  return res;
1419 }
1420 
1421 namespace crn
1422 {
1423  /**************************************************************************************
1424  * Characterization
1425  *************************************************************************************/
1426 
1431  template<typename T> bool IsBitonal(const Image<T> &img)
1432  {
1433  auto values = std::set<T>{};
1434  for (const auto &p : img)
1435  {
1436  values.insert(p);
1437  if (values.size() > 2)
1438  return false;
1439  }
1440  return true;
1441  }
1442 
1447  template<typename T, typename CMP> std::pair<T, T> MinMax(const Image<T> &img, CMP cmp)
1448  {
1449  auto mM = std::minmax_element(img.begin(), img.end(), cmp);
1450  return std::make_pair(*mM.first, *mM.second);
1451  }
1452 
1459  template<typename T> Rect AutoCrop(const Image<T> &img, const T &bgval)
1460  {
1461  auto bbox = img.GetBBox();
1462  // left
1463  for (int x = bbox.GetLeft(); x <= bbox.GetRight(); ++x)
1464  {
1465  bool bl = false;
1466  for (int y = bbox.GetTop(); y <= bbox.GetBottom(); ++y)
1467  if (img.At(x, y) != bgval)
1468  {
1469  bl = true;
1470  break;
1471  }
1472  if (bl)
1473  {
1474  bbox.SetLeft(x);
1475  break;
1476  }
1477  }
1478  // right
1479  for (int x = bbox.GetRight(); x >= bbox.GetLeft(); --x)
1480  {
1481  bool bl = false;
1482  for (int y = bbox.GetTop(); y <= bbox.GetBottom(); ++y)
1483  if (img.At(x, y) != bgval)
1484  {
1485  bl = true;
1486  break;
1487  }
1488  if (bl)
1489  {
1490  bbox.SetRight(x);
1491  break;
1492  }
1493  }
1494  // top
1495  for (int y = bbox.GetTop(); y <= bbox.GetBottom(); ++y)
1496  {
1497  bool bl = false;
1498  for (int x = bbox.GetLeft(); x <= bbox.GetRight(); ++x)
1499  if (img.At(x, y) != bgval)
1500  {
1501  bl = true;
1502  break;
1503  }
1504  if (bl)
1505  {
1506  bbox.SetTop(y);
1507  break;
1508  }
1509  }
1510  // bottom
1511  for (int y = bbox.GetBottom(); y >= bbox.GetTop(); --y)
1512  {
1513  bool bl = false;
1514  for (int x = bbox.GetLeft(); x <= bbox.GetRight(); ++x)
1515  if (img.At(x, y) != bgval)
1516  {
1517  bl = true;
1518  break;
1519  }
1520  if (bl)
1521  {
1522  bbox.SetBottom(y);
1523  break;
1524  }
1525  }
1526  return bbox;
1527  }
1528 
1535  template<typename T> Image<T> MakeAutoCrop(const Image<T> &img, const T &bgval)
1536  {
1537  return Image<T>(img, AutoCrop(img, bgval));
1538  }
1539 
1550  template<typename T, typename Y> Point2DInt CrossCorrelation(const Image<T> &img1, const Image<Y> &img2, T fill1, Y fill2)
1551  {
1552  // find nearest power of 2 dimensions
1553  auto w = crn::Max(img1.GetWidth(), img2.GetWidth());
1554  auto h = crn::Max(img1.GetHeight(), img2.GetHeight());
1555  double logs = log2(int(w));
1556  while (logs != ceil(logs))
1557  {
1558  w += 1;
1559  logs = log2(int(w));
1560  }
1561  if (w < 2) w = 2;
1562  logs = log2(int(h));
1563  while (logs != ceil(logs))
1564  {
1565  h += 1;
1566  logs = log2(int(h));
1567  }
1568  if (h < 2) h = 2;
1569  // compute cross correlation
1570  auto c1 = MatrixComplex(h, w, std::complex<double>(fill1));
1571  for (size_t r = 0; r < img1.GetHeight(); ++r)
1572  for (size_t c = 0; c < img1.GetWidth(); ++c)
1573  c1[r][c] = std::complex<double>(img1.At(c, r));
1574  c1.FFT(true);
1575  auto c2 = MatrixComplex(h, w, std::complex<double>(fill2));
1576  for (size_t r = 0; r < img2.GetHeight(); ++r)
1577  for (size_t c = 0; c < img2.GetWidth(); ++c)
1578  c2[r][c] = std::complex<double>(img2.At(c, r));
1579  c2.FFT(true);
1580  for (size_t r = 0; r < h; ++r)
1581  for (size_t c = 0; c < w; ++c)
1582  c1[r][c] *= std::conj(c2[r][c]);
1583  c2 = MatrixComplex(1, 1); // free memory
1584  c1.FFT(true);
1585  // find best match
1586  Point2DInt p;
1587  double maxc = 0;
1588  for (size_t r = 0; r < h; ++r)
1589  for (size_t c = 0; c < w; ++c)
1590  {
1591  double corr = std::norm(c1[r][c]);
1592  if (corr > maxc)
1593  {
1594  p.X = int(c);
1595  p.Y = int(r);
1596  maxc = corr;
1597  }
1598  }
1599  int hw = (int)w/2;
1600  int hh = (int)h/2;
1601  if (p.X >= hw)
1602  p.X = (p.X % hw) - hw;
1603  if (p.Y >= hh)
1604  p.Y = (p.Y % hh) - hh;
1605  return p;
1606  }
1607 
1608  /**************************************************************************************
1609  * Transformation
1610  *************************************************************************************/
1611 
1613  template<typename T> SummedAreaTable<SumType<T>> MakeSummedAreaTable(const Image<T> &img)
1614  {
1615  using SumType = SumType<T>;
1616  auto sum = SummedAreaTable<SumType>(img.GetWidth(), img.GetHeight());
1617  // first element
1618  sum.SetValue(0, 0, SumType(img.At(0)));
1619  // first line
1620  for (size_t x = 1; x < img.GetWidth(); ++x)
1621  sum.SetValue(x, 0, sum.GetValue(x - 1, 0) + SumType(img.At(x)));
1622  // first column
1623  for (size_t y = 1; y < img.GetHeight(); ++y)
1624  sum.SetValue(0, y, sum.GetValue(0, y - 1) + SumType(img.At(0, y)));
1625  // table: I(x,y) = i(x,y) + I(x-1,y) + I(x,y-1) - I(x-1,y-1)
1626  for (size_t y = 1; y < img.GetHeight(); ++y)
1627  for (size_t x = 1; x < img.GetWidth(); ++x)
1628  sum.SetValue(x, y, SumType(img.At(x, y)) + sum.GetValue(x - 1, y) + sum.GetValue(x, y - 1) - sum.GetValue(x - 1, y - 1));
1629  return sum;
1630  }
1631 
1632  namespace impl
1633  {
1642  template<typename T> void ShiftCopyRow(Image<T> &dest, const Image<T> &src, size_t row, int offset, double prevWeight, const T &bgColor)
1643  {
1644  size_t rowoffset = row * dest.GetWidth();
1645 
1646  // fill the possible left gap
1647  for (int x = 0; x < offset; ++x)
1648  dest.At(x + rowoffset) = bgColor;
1649 
1650  auto oldPrevPix = DecimalType<T>(bgColor) * prevWeight;
1651  auto srcrowoffset = row * src.GetWidth();
1652  for (size_t x = 0; x < src.GetWidth(); ++x)
1653  { // scan original pixels
1654  T pix = src.At(x + srcrowoffset); // do not use auto because of vector<bool>!
1655  auto prevPix = DecimalType<T>(pix) * prevWeight;
1656  // compute value to copy
1657  auto tmp = DecimalType<T>(pix);
1658  tmp += oldPrevPix;
1659  tmp -= prevPix;
1660  pix = T(tmp);
1661  // copy the new value
1662  if ((x + offset >= 0) && (x + offset < dest.GetWidth()))
1663  {
1664  dest.At(x + offset + rowoffset) = pix;
1665  }
1666  // remember the weighted previous pixel
1667  oldPrevPix = prevPix;
1668  }
1669  // fill the possible right gap
1670  size_t x = src.GetWidth() + offset;
1671  if (x < dest.GetWidth())
1672  {
1673  auto prevPix = DecimalType<T>(bgColor) * prevWeight;
1674  auto pix = bgColor;
1675  auto tmp = DecimalType<T>(pix);
1676  tmp += oldPrevPix;
1677  tmp -= prevPix;
1678  pix = T(tmp);
1679  dest.At(x + rowoffset) = pix;
1680  }
1681  for (++x; x < dest.GetWidth(); ++x)
1682  {
1683  dest.At(x + rowoffset) = bgColor;
1684  }
1685  }
1686 
1695  template<typename T> void ShiftCopyColumn(Image<T> &dest, const Image<T> &src, size_t col, int offset, double prevWeight, const T &bgColor)
1696  {
1697  // fill the possible upper gap
1698  for (int y = 0; y < offset; ++y)
1699  {
1700  dest.At(col, y) = bgColor;
1701  }
1702 
1703  auto oldPrevPix = DecimalType<T>(bgColor) * prevWeight;
1704  for (size_t y = 0; y < src.GetHeight(); ++y)
1705  { // scan original pixels
1706  T pix = src.At(col, y); // do not use auto because of vector<bool>!
1707  auto prevPix = DecimalType<T>(pix) * prevWeight;
1708  // compute value to copy
1709  auto tmp = DecimalType<T>(pix);
1710  tmp -= prevPix;
1711  tmp += oldPrevPix;
1712  pix = T(tmp);
1713  // copy the new value
1714  if ((y + offset >= 0) && (y + offset < dest.GetHeight()))
1715  {
1716  dest.At(col, y + offset) = pix;
1717  }
1718  // remember the weighted previous pixel
1719  oldPrevPix = prevPix;
1720  }
1721  // fill the possible bottom gap
1722  size_t y = src.GetHeight() + offset;
1723  if (y < dest.GetHeight())
1724  {
1725  auto prevPix = DecimalType<T>(bgColor) * prevWeight;
1726  auto pix = bgColor;
1727  auto tmp = DecimalType<T>(pix);
1728  tmp += oldPrevPix;
1729  tmp -= prevPix;
1730  pix = T(tmp);
1731  dest.At(col, y) = pix;
1732  }
1733  for (++y; y < dest.GetHeight(); ++y)
1734  {
1735  dest.At(col, y) = bgColor;
1736  }
1737  }
1738 
1744  template<typename T> Image<T> MakeSmallRotation(const Image<T> &img, const Angle<Degree> &angle, const T &bgColor)
1745  {
1746  if (angle.value == 0.0)
1747  { // if no rotation, simply copy
1748  return img;
1749  }
1750  const auto radAngle = angle.Get<Radian>();
1751  const auto rotCos = cos(radAngle);
1752  const auto rotSin = sin(radAngle);
1753  const auto rotTan = tan(radAngle / 2.0);
1754 
1755  // first shear (horizontal)
1756  auto shear1 = Image<T>(img.GetWidth() + int(double(img.GetHeight()) * Abs(rotTan)), img.GetHeight());
1757  for (size_t y = 0; y < shear1.GetHeight(); ++y)
1758  {
1759  auto shear = 0.0;
1760  if (rotTan >= 0.0)
1761  { // angle >= 0
1762  shear = (double(y) + 0.5) * rotTan;
1763  }
1764  else
1765  { // angle < 0
1766  shear = (double(y) - double(shear1.GetHeight()) + 0.5) * rotTan;
1767  }
1768  ShiftCopyRow(shear1, img, y, int(floor(shear)), shear - floor(shear), bgColor);
1769  }
1770 
1771  // second shear (vertical)
1772  auto shear2 = Image<T>(shear1.GetWidth(), size_t(double(img.GetWidth()) * Abs(rotSin) + double(img.GetHeight()) * rotCos) + 1);
1773  auto shear = 0.0;
1774  if (rotSin > 0.0)
1775  { // angle > 0
1776  shear = (double(img.GetWidth()) - 1.0) * rotSin;
1777  }
1778  else
1779  { // angle <= 0
1780  shear = double(shear1.GetWidth() - img.GetWidth()) * rotSin;
1781  }
1782  for (size_t x = 0; x < shear2.GetWidth(); ++x, shear -= rotSin)
1783  {
1784  ShiftCopyColumn(shear2, shear1, x, int(floor(shear)), shear - floor(shear), bgColor);
1785  }
1786  shear1 = Image<T>{}; // free a bit of memory
1787 
1788  // third sear (horizontal)
1789  auto shear3 = Image<T>(size_t(double(img.GetHeight()) * Abs(rotSin) + double(img.GetWidth()) * rotCos) + 1, shear2.GetHeight());
1790  if (rotSin >= 0.0)
1791  { // angle >= 0
1792  shear = (1.0 - double(img.GetWidth())) * rotSin * rotTan;
1793  }
1794  else
1795  {
1796  shear = ((1.0 - double(img.GetWidth())) * rotSin + 1.0 - double(shear3.GetHeight())) * rotTan;
1797  }
1798  for (size_t y = 0; y < shear3.GetHeight(); ++y, shear += rotTan)
1799  {
1800  ShiftCopyRow(shear3, shear2, y, int(floor(shear)), shear - floor(shear), bgColor);
1801  }
1802 
1803  return shear3;
1804  }
1805  }
1806 
1813  template<typename T> Image<T> MakeRotation(const Image<T> &img, const Angle<Degree> &angle, const T &bgColor)
1814  {
1815  auto rot = angle;
1816  while (rot.value >= 360.0)
1817  rot.value -= 360.0;
1818  while (rot.value < 0.0)
1819  rot.value += 360.0;
1820  if (rot.value == 0)
1821  return img;
1822  else if (rot.value == 90)
1823  return Make90Rotation(img);
1824  else if (rot.value == 180)
1825  return Make180Rotation(img);
1826  else if (rot.value == 270)
1827  return Make270Rotation(img);
1828 
1829  auto midimg = &img;
1830  auto tmpbuff = Image<T>{};
1831  if (rot.value > 225.0)
1832  { // in ]225, 360[, rotate 270 then small rotation
1833  tmpbuff = Make270Rotation(img);
1834  midimg = &tmpbuff;
1835  rot.value -= 270.0;
1836  }
1837  else if (rot.value > 135.0)
1838  { // in ]135, 225], rotate 180 then small rotation
1839  tmpbuff = Make180Rotation(img);
1840  midimg = &tmpbuff;
1841  rot.value -= 180.0;
1842  }
1843  else if (rot.value > 45.0)
1844  { // in ]45, 135], rotate 90 then small rotation
1845  tmpbuff = Make90Rotation(img);
1846  midimg = &tmpbuff;
1847  rot.value -= 90.0;
1848  }
1849  // do the final rotation
1850  return impl::MakeSmallRotation(*midimg, rot, bgColor);
1851  }
1852 
1857  template<typename T> Image<T> Make90Rotation(const Image<T> &img)
1858  {
1859  auto newi = Image<T>(img.GetHeight(), img.GetWidth());
1860  FOREACHPIXEL(x, y, img)
1861  newi.At(y, img.GetWidth() - 1 - x) = img.At(x, y);
1862  return newi;
1863  }
1864 
1869  template<typename T> Image<T> Make180Rotation(const Image<T> &img)
1870  {
1871  auto newi = Image<T>(img.GetWidth(), img.GetHeight());
1872  const auto size = img.Size();
1873  for (size_t tmp = 0; tmp < size; ++tmp)
1874  newi.At(size - tmp) = img.At(tmp);
1875  return newi;
1876  }
1877 
1882  template<typename T> Image<T> Make270Rotation(const Image<T> &img)
1883  {
1884  auto newi = Image<T>(img.GetHeight(), img.GetWidth());
1885  FOREACHPIXEL(x, y, img)
1886  newi.At(img.GetHeight() - 1 - y, x) = img.At(x, y);
1887  return newi;
1888  }
1889 
1890 
1891 }
1892 
1893 #endif
void ShiftCopyRow(Image< T > &dest, const Image< T > &src, size_t row, int offset, double prevWeight, const T &bgColor)
Copies a row with a shift.
Definition: CRNImage.hpp:1642
Image< T > MakeSmallRotation(const Image< T > &img, const Angle< Degree > &angle, const T &bgColor)
Definition: CRNImage.hpp:1744
typename TypeInfo< T >::SumType SumType
Definition: CRNType.h:185
void GaussianBlur(double sigma)
Gaussian blur.
Definition: CRNImage.hpp:1132
Abstract class for images.
Definition: CRNImage.h:141
void Dilate(const MatrixInt &strel, CMP cmp=std::less< pixel_type >{})
Morphological dilatation.
Definition: CRNImage.hpp:761
crn::Image< crn::SumType< typename std::common_type< T1, T2 >::type > > operator*(const crn::Image< T1 > &i1, const crn::Image< T2 > &i2)
Definition: CRNImage.hpp:1385
ScalarRange< T > Range(T b, T e)
Creates a range [[b, e[[.
Definition: CRNType.h:257
size_t GetRows() const noexcept
Returns the number of rows.
Definition: CRNMatrix.h:157
size_t GetCols() const noexcept
Returns the number of columns.
Definition: CRNMatrix.h:163
void ShiftCopyColumn(Image< T > &dest, const Image< T > &src, size_t col, int offset, double prevWeight, const T &bgColor)
Copies a column with a shift.
Definition: CRNImage.hpp:1695
void ScanFill(size_t x, size_t y, const pixel_type &val, crn::DistanceType dist=crn::DistanceType::D4)
Fills a portion of the image.
Definition: CRNImage.hpp:363
Orientation
An enumeration of orientations.
Definition: CRNMath.h:152
void Assign(const Image< Y > &img)
Force copy operator (pixel cast)
Definition: CRNImage.hpp:147
Radian angle unit.
std::vector< pixel_type >::reference At(size_t x, size_t y) noexcept
Returns a reference to a pixel.
Definition: CRNImage.h:224
Image< T > MakeRotation(const Image< T > &img, const Angle< Degree > &angle, const T &bgColor)
Creates a rotated version of the image.
Definition: CRNImage.hpp:1813
size_t GetHeight() const noexcept
Definition: CRNImage.h:74
std::vector< pixel_type > pixels
Definition: CRNImage.h:295
static MatrixDouble NewGaussianLine(double sigma)
Creates a line matrix with a centered Gaussian.
const T & Max(const T &a, const T &b)
Returns the max of two values.
Definition: CRNMath.h:47
Rect AutoCrop(const Image< T > &img, const T &bgval)
Estimates the ideal crop for the image.
Definition: CRNImage.hpp:1459
void Flip(const Orientation &ori)
Flips the image.
Definition: CRNImage.hpp:729
void SavePNG(const ImageBW &img, const Path &fname)
Saves as PNG file.
crn::Image< crn::SumType< typename std::common_type< T1, T2 >::type > > operator/(const crn::Image< T1 > &i1, const crn::Image< T2 > &i2)
Definition: CRNImage.hpp:1410
std::vector< pixel_type >::iterator end()
Definition: CRNImage.h:213
#define FOREACHPIXEL(x, y, img)
Convenience macro to sweep an image.
Definition: CRNImage.h:37
int GetTop() const
Returns the topmost coordinate.
Definition: CRNRect.h:107
int GetLeft() const
Returns the leftmost coordinate.
Definition: CRNRect.h:99
void Erode(const MatrixInt &strel, CMP cmp=std::less< pixel_type >{})
Morphological erosion.
Definition: CRNImage.hpp:827
A convenience class for angles units.
Image & operator*=(double f)
Multiplies all pixels.
Definition: CRNImage.hpp:221
Integer matrix class.
Definition: CRNMatrixInt.h:40
crn::Image< crn::DiffType< typename std::common_type< T1, T2 >::type > > operator-(const crn::Image< T1 > &i1, const crn::Image< T2 > &i2)
Definition: CRNImage.hpp:1374
A generic domain error.
Definition: CRNException.h:83
Image< T > Make270Rotation(const Image< T > &img)
Creates a rotated version of the image.
Definition: CRNImage.hpp:1882
void FastDilate(size_t halfwin, size_t index=0, CMP cmp=std::less< pixel_type >{})
Morphological dilatation.
Definition: CRNImage.hpp:890
size_t width
Definition: CRNImage.h:92
int SetWidth(int wid)
Changes the width of the rectangle.
Definition: CRNRect.h:210
Image & operator/=(const Image &img)
Divides by another image's pixels.
Definition: CRNImage.hpp:247
Image & operator=(const Image &img)=default
Copy operator.
void SetValue(size_t x, size_t y, T val)
Alters the value of a bin.
A convenience class for file paths.
Definition: CRNPath.h:39
U::type Get() const noexcept(std::is_nothrow_constructible< typename U::type >::value &&std::is_nothrow_constructible< typename TypeInfo< typename Unit::type >::SumType >::value)
Conversion to any other angle unit.
bool IsBitonal(const Image< T > &img)
Is the image binary (black & white)?
Definition: CRNImage.hpp:1431
Image & operator+=(const Image &img)
Adds another image.
Definition: CRNImage.hpp:194
SummedAreaTable< SumType< T > > MakeSummedAreaTable(const Image< T > &img)
Creates a summed area table of the image.
Definition: CRNImage.hpp:1613
value_type X
Definition: CRNPoint2D.h:63
crn::Image< crn::SumType< typename std::common_type< T1, T2 >::type > > operator+(const crn::Image< T1 > &i1, const crn::Image< T2 > &i2)
Definition: CRNImage.hpp:1363
Image & operator-=(const Image &img)
Subtracts another image.
Definition: CRNImage.hpp:208
Image< T > Make90Rotation(const Image< T > &img)
Creates a rotated version of the image.
Definition: CRNImage.hpp:1857
Point2DInt CrossCorrelation(const Image< T > &img1, const Image< Y > &img2, T fill1=T(0), Y fill2=Y(0))
Best match between two images.
Definition: CRNImage.hpp:1550
A dimension error.
Definition: CRNException.h:119
double matrix class
virtual void SavePNG(const Path &fname) const override
Saves as PNG file.
Definition: CRNImage.hpp:160
void Convolve(const MatrixDouble &mat)
Convolution.
Definition: CRNImage.hpp:1023
A summed area table.
typename TypeInfo< T >::DecimalType DecimalType
Definition: CRNType.h:187
void DrawLine(size_t x1, size_t y1, size_t x2, size_t y2, pixel_type color)
Draws a line using a specified color.
Definition: CRNImage.hpp:543
size_t height
Definition: CRNImage.h:92
virtual void SaveJPEG(const Path &fname, unsigned int qual) const override
Saves as JPEG file.
Definition: CRNImage.hpp:169
Unit::type value
void Abs(Image< T > &img, typename std::enable_if< std::is_arithmetic< T >::value >::type *dummy=nullptr) noexcept
Replaces each pixel by its absolute value.
Definition: CRNImageGray.h:47
typename TypeInfo< T >::DiffType DiffType
Definition: CRNType.h:186
DistanceType
An enumeration of distances.
Definition: CRNMath.h:159
std::vector< pixel_type >::iterator begin()
Definition: CRNImage.h:211
const T & Min(const T &a, const T &b)
Returns the min of two values.
Definition: CRNMath.h:49
Image()
Default constructor.
Definition: CRNImage.h:153
value_type Y
Definition: CRNPoint2D.h:63
Rect GetBBox() const noexcept
Definition: CRNImage.cpp:56
void Complement(pixel_type maxval=std::numeric_limits< pixel_type >::max())
Complement.
Definition: CRNImage.hpp:270
void Blit(const Image< Y > &src, const Rect &srczone, size_t dx, size_t dy)
Copies a part of an image.
Definition: CRNImage.hpp:282
std::pair< T, T > MinMax(const Image< T > &img, CMP cmp=CMP{})
Returns min and max pixel values.
Definition: CRNImage.hpp:1447
void SaveJPEG(const ImageBW &img, const Path &fname, unsigned int qual)
Saves as JPEG file.
Complex matrix class.
size_t GetWidth() const noexcept
Definition: CRNImage.h:72
void Negative()
Negative.
Definition: CRNImage.hpp:261
size_t Size() const noexcept
Definition: CRNImage.h:78
Image< T > MakeAutoCrop(const Image< T > &img, const T &bgval)
Creates a new image as the ideal crop for the image.
Definition: CRNImage.hpp:1535
T pixel_type
Definition: CRNImage.h:144
Image< T > Make180Rotation(const Image< T > &img)
Creates a rotated version of the image.
Definition: CRNImage.hpp:1869
void DrawRect(const Rect &r, pixel_type color, bool filled=false)
Draws a rectangle using a specified color.
Definition: CRNImage.hpp:515
A 2D point class.
Definition: CRNPoint2DInt.h:39
virtual void ScaleToSize(size_t w, size_t h) override
Scales the image.
Definition: CRNImage.hpp:614
bool operator==(const Image &other) const
Tests equality.
Definition: CRNImage.hpp:178
Base class for images.
Definition: CRNImage.h:46
void FastErode(size_t halfwin, size_t index=0, CMP cmp=std::less< pixel_type >{})
Morphological erosion.
Definition: CRNImage.hpp:956
Invalid argument error (e.g.: nullptr pointer)
Definition: CRNException.h:107
void FloodFill(size_t x, size_t y, const pixel_type &val, crn::DistanceType dist=crn::DistanceType::D4)
Flood fills a portion of the image.
Definition: CRNImage.hpp:311
int SetLeft(int begX) noexcept
Changes the leftmost coordinate.
Definition: CRNRect.h:152
A rectangle class.
Definition: CRNRect.h:46