CVD 0.8
cvd/morphology.h
00001 #ifndef CVD_INCLUDE_MORPHOLOGY_H
00002 #define CVD_INCLUDE_MORPHOLOGY_H
00003 
00004 #include <cvd/vision_exceptions.h>
00005 #include <cvd/vision.h>
00006 #include <cvd/image.h>
00007 #include <map>
00008 #include <vector>
00009 
00010 
00011 namespace CVD
00012 {
00013     #ifndef DOXYGEN_IGNORE_INTERNAL
00014     namespace Internal
00015     {
00016         namespace MorphologyHelpers
00017         {
00018             using namespace std;
00019             
00020             //Compute pointer offsets for a bunch of ImageRef offsets.
00021             template<class T> vector<ptrdiff_t> offsets(const vector<ImageRef>& v, const SubImage<T>& s)
00022             {
00023                 vector<ptrdiff_t> off;
00024 
00025                 for(unsigned int i=0; i < v.size(); i++)
00026                     off.push_back(v[i].x + v[i].y * s.row_stride()-1);
00027                 return off;
00028             }
00029             
00030             //Split a list of ImageRefs up in to rows.
00031             vector<vector<ImageRef> > row_split(const vector<ImageRef>& v, int y_lo, int y_hi);
00032         }
00033     }
00034     #endif
00035 
00075     template<class Accumulator, class T>
00076     void morphology(const SubImage<T>& in, const std::vector<ImageRef>& selem, const Accumulator& a_, SubImage<T>& out)
00077     {
00078         using Internal::MorphologyHelpers::offsets; 
00079         using Internal::MorphologyHelpers::row_split;   
00080         using std::min;
00081         using std::max;
00082         using std::vector;
00083 
00084         if(in.size() != out.size())
00085             throw Exceptions::Vision::IncompatibleImageSizes(__FUNCTION__);
00086 
00087         //Cases are:
00088         //
00089         // Small selem compared to image:
00090         //   Topleft corner, top row, top right corner
00091         //   left edge, centre, right edge
00092         //   etc
00093         
00095         //Find the extents of the structuring element
00096         int x_lo = selem[0].x;
00097         int x_hi = selem[0].x;
00098         int y_lo = selem[0].y;
00099         int y_hi = selem[0].y;
00100 
00101         for(unsigned int i=0; i < selem.size(); i++)
00102         {
00103             x_lo = min(x_lo, selem[i].x);
00104             x_hi = max(x_hi, selem[i].x);
00105             y_lo = min(y_lo, selem[i].y);
00106             y_hi = max(y_hi, selem[i].y);
00107         }
00108 
00110         //Shift the  structure element by one and find the differeneces
00111         vector<ImageRef> structure_element = selem;
00112         vector<ImageRef> shifted;
00113 
00114         sort(structure_element.begin(), structure_element.end());
00115         for(unsigned int i=0; i < structure_element.size(); i++)
00116             shifted.push_back(structure_element[i] + ImageRef(1, 0));
00117         
00118         vector<ImageRef> add, remove;
00119         set_difference(shifted.begin(), shifted.end(), structure_element.begin(), structure_element.end(), back_inserter(add));
00120         set_difference(structure_element.begin(), structure_element.end(), shifted.begin(), shifted.end(), back_inserter(remove));
00121 
00123         //  
00124         //Compute the integer offsets to pixels for speed;
00125         vector<ptrdiff_t> add_off = offsets(add, in);
00126         vector<ptrdiff_t> remove_off = offsets(remove, in);
00127         
00128 
00130         //  
00131         // Split by rows, to make the top and bottom edges easier.
00132         //
00133         //Because of set operations, the ImageRefs are ordered within each row.
00134         vector<vector<ImageRef> > split_selem = row_split(structure_element, y_lo, y_hi);
00135         vector<vector<ImageRef> > split_add   = row_split(add, y_lo, y_hi);
00136         vector<vector<ImageRef> > split_remove= row_split(remove, y_lo, y_hi);
00137         
00138         Accumulator acc(a_);
00139         //If the image is at least as wide as the structuring element
00140         if(x_hi - x_lo + 1 <= in.size().x)
00141             for(int y=0; y < in.size().y; y++)
00142             {
00143                 //Find the rows which overlap with the image. Only work with these rows.
00144                 int startrow = max(0, - y_lo - y);
00145                 int endrow =   split_selem.size() - max(0, y + y_hi - in.size().y+1);
00146                 
00147                 //Figure out the range of the "easy" bit.
00148                 int x_first_full = max(0, -x_lo);                               //This is the first position at which we have a full kernel in the image
00149                 int x_after_last_full = min(in.size().x, in.size().x - x_hi);   //This is one beyone the end of the position where the last kernel fits in the image.
00150 
00151                 //Clear the  accumulator
00152                 acc.clear();
00153 
00154                 //Fill in the accumulator sitting up against the left hand side of the image.
00155                 for(int i=startrow; i < endrow; i++)
00156                     for(int j=(int)split_selem[i].size()-1; j >=0 && split_selem[i][j].x >=0; j--)
00157                         acc.insert(in[y + split_selem[i][0].y][split_selem[i][j].x]);
00158                 
00159                 out[y][0] = acc.get();
00160 
00161                 //Shift the kernel until we get to the point where
00162                 //we can start shifting the kernel without testing to 
00163                 //see it fits withing the image width.
00164                 for(int x=1; x <= x_first_full ; x++)
00165                 {
00166                     for(int i=startrow; i < endrow; i++)
00167                         for(int j=(int)split_remove[i].size()-1; j >=0 && split_remove[i][j].x+x-1 >=0; j--)
00168                             acc.remove(in[y + split_remove[i][0].y][x+split_remove[i][j].x-1]);
00169 
00170                     for(int i=startrow; i < endrow; i++)
00171                         for(int j=(int)split_add[i].size()-1; j >=0 && split_add[i][j].x+x-1 >=0; j--)
00172                             acc.insert(in[y + split_add[i][0].y][x+split_add[i][j].x-1]);
00173 
00174                     out[y][x] = acc.get();
00175                 }
00176 
00177                 //Go through the two incremental kernels to figure out which
00178                 //indices are fit within the image. This removes a test from
00179                 //the following shift section.
00180                 int add_start=0, add_end=0, remove_start=0, remove_end=0;
00181                 for(int i=0; i < startrow; i++)
00182                 {
00183                     add_start+=split_add[i].size();
00184                     remove_start+=split_remove[i].size();
00185                 }
00186                 for(int i=0; i < endrow; i++)
00187                 {
00188                     add_end+=split_add[i].size();
00189                     remove_end+=split_remove[i].size();
00190                 }
00191 
00192                 //Shift the kernel in the area which requires no tests.
00193                 for(int x=max(0, -x_lo+1); x < x_after_last_full; x++)
00194                 {
00195                     for(int i=remove_start; i < remove_end; i++)
00196                         acc.remove(*(in[y] + x + remove_off[i]));
00197                     
00198                     for(int i=add_start; i < add_end; i++)
00199                         acc.insert(*(in[y] + x + add_off[i]));
00200                     
00201                     out[y][x] = acc.get();
00202                 }
00203                 
00204                 //Now perform the right hand edge
00205                 for(int x=x_after_last_full; x < in.size().x ; x++)
00206                 {
00207                     for(int i=startrow; i < endrow; i++)
00208                         for(int j=0; j < (int)split_remove[i].size() && split_remove[i][j].x+x-1 < in.size().x; j++)
00209                             acc.remove(in[y + split_remove[i][0].y][x+split_remove[i][j].x-1]);
00210 
00211                     for(int i=startrow; i < endrow; i++)
00212                         for(int j=0; j < (int)split_add[i].size() && split_add[i][j].x+x-1 < in.size().x; j++)
00213                             acc.insert(in[y + split_add[i][0].y][x+split_add[i][j].x-1]);
00214 
00215                     out[y][x] = acc.get();
00216                 }
00217             }
00218         else
00219         {
00220             //The image is too narrow to have a clear area in the middle.
00221 
00222             for(int y=0; y < in.size().y; y++)
00223             {
00224                 //Find the rows which overlap with the image. Only work with these rows.
00225                 int startrow = max(0, - y_lo - y);
00226                 int endrow =   split_selem.size() - max(0, y + y_hi - in.size().y+1);
00227                 
00228                 //Clear the accumulator
00229                 acc.clear();
00230 
00231                 //Fill in the accumulator sitting up against the left hand side of the image.
00232                 for(int i=startrow; i < endrow; i++)
00233                     for(int j=0; j < (int)split_selem[i].size(); j++)
00234                     {
00235                         int xp = split_selem[i][j].x;
00236                         if(xp >= 0 && xp < in.size().x)
00237                             acc.insert(in[y + split_selem[i][0].y][xp]);
00238                     }
00239                 
00240                 out[y][0] = acc.get();
00241 
00242                 //Shift the kernel using the incrementals
00243                 for(int x=1; x < in.size().x ; x++)
00244                 {
00245                     for(int i=startrow; i < endrow; i++)
00246                         for(int j=0; j < (int)split_remove[i].size(); j++)
00247                         {
00248                             int xp = x + split_remove[i][j].x - 1;
00249                             if(xp >= 0 && xp < in.size().x)
00250                                 acc.remove(in[y + split_add[i][0].y][xp]);
00251                         }
00252 
00253                     for(int i=startrow; i < endrow; i++)
00254                         for(int j=0; j < (int)split_add[i].size(); j++)
00255                         {
00256                             int xp = x + split_add[i][j].x - 1;
00257                             if(xp >= 0 && xp < in.size().x)
00258                                 acc.insert(in[y + split_add[i][0].y][xp]);
00259                         }
00260 
00261                     out[y][x] = acc.get();
00262                 }
00263             }
00264         }
00265     }
00266 
00267     #ifndef DOXYGEN_IGNORE_INTERNAL
00268         namespace Internal
00269         {
00270             template<class C, class D> class PerformMorphology{};
00271 
00272             template<class C, class D>  struct ImagePromise<PerformMorphology<C, D> >
00273             {
00274                 ImagePromise(const SubImage<C>& im, const D& acc, const std::vector<ImageRef>& s_)
00275                 :i(im),a(acc),s(s_)
00276                 {}
00277 
00278                 const SubImage<C>& i;
00279                 const D& a;
00280                 const std::vector<ImageRef>& s;
00281 
00282                 template<class E> void execute(Image<E>& j)
00283                 {
00284                     j.resize(i.size());
00285                     morphology(i, s, a, j);
00286                 }
00287             };
00288         };
00289 
00290         template<class C, class D> Internal::ImagePromise<Internal::PerformMorphology<C, D> > morphology(const SubImage<C>& c, const std::vector<ImageRef>& selem, const D& a)
00291         {
00292             return Internal::ImagePromise<Internal::PerformMorphology<C, D> >(c, a, selem);
00293         }
00294     #else
00295         
00303         Image<T> morphology(const SubImage<T>& in, const std::vector<ImageRef>& selem, const Accumulator& a_);
00304 
00305     #endif
00306 
00309     namespace Morphology
00310     {
00315         template<class T, template<class> class Cmp> struct BasicGray
00316         {
00317             private:
00318                 std::map<T, int, Cmp<T> > pix;
00319             
00320             public:
00321                 void clear()
00322                 {
00323                     pix.clear();
00324                 }
00325                 void insert(const T& t)
00326                 {
00327                     pix[t]++;
00328                 }
00329 
00330                 void remove(const T& t)
00331                 {
00332                     --pix[t];
00333                 }
00334 
00335                 T get()
00336                 {
00337                     typedef typename std::map<T, int, Cmp<T> >::iterator it;
00338 
00339                     for(it i=pix.begin(); i != pix.end();)
00340                     {
00341                         it old = i;
00342                         i++;
00343 
00344                         if(old->second == 0)
00345                             pix.erase(old);
00346                         else
00347                             return old->first;
00348                     }
00349 
00350                     assert(0);
00351                     return 0;
00352                 }
00353         };
00354         
00355 
00358         template<class T> class Erode: public BasicGray<T, std::less>
00359         {
00360         };
00361 
00362 
00365         template<class T> class Dilate: public BasicGray<T, std::greater>
00366         {
00367         };
00368 
00369 
00372         template<class T> class Percentile;
00373 
00376         template<class T> class Median;
00377 
00382         struct BasicGrayByte
00383         {
00384             protected:
00385                 int histogram[256];
00386                 int total;
00387             
00388             public:
00389                 BasicGrayByte()
00390                 {
00391                     clear();
00392                 }
00393 
00394                 void clear()
00395                 {
00396                     total=0;
00397                     for(int i=0; i < 256; i++)
00398                         histogram[i] = 0;
00399                 }
00400 
00401                 void insert(byte t)
00402                 {
00403                     total++;
00404                     histogram[t]++;
00405                 }
00406 
00407                 void remove(byte t)
00408                 {
00409                     total--;
00410                     histogram[t]--;
00411                 }
00412         };
00413 
00414 
00417         template<> class Erode<byte>: public BasicGrayByte
00418         {
00419             public:
00420                 byte get()
00421                 {
00422                     for(int j=0; j < 256; j++)
00423                         if(histogram[j])
00424                             return j;
00425                     
00426                     assert(0);
00427                     return 0;
00428                 }
00429         };
00430 
00433         template<> class Dilate<byte>: public BasicGrayByte
00434         {
00435             public:
00436                 byte get()
00437                 {
00438                     for(int j=255; j >=0 ; j--)
00439                         if(histogram[j])
00440                             return j;
00441                     
00442                     assert(0);
00443                     return 0;
00444                 }
00445         };
00446 
00449         template<> class Percentile<byte>: public BasicGrayByte
00450         {
00451             private:
00452                 double ptile;
00453 
00454             public:
00455                 Percentile(double p)
00456                 :ptile(p)
00457                 {
00458                 }
00459 
00460                 byte get()
00461                 {
00462                     using std::max;
00463 
00464                     if(ptile < 0.5)
00465                     {
00466                         int sum=0;
00467                         //because we use a > test below (to work for the 0th ptile)
00468                         //we have to use the scaled threshold -1 otherwise it will
00469                         //not work for the 100th percentile.
00470                         int threshold = max(0, (int)floor(total * ptile+.5)- 1);
00471 
00472                         for(int j=0; j < 255; j++)
00473                         {
00474                             sum += histogram[j];
00475 
00476                             if(sum > threshold)
00477                                 return j;
00478                         }
00479                         
00480                         return 255;
00481                     }
00482                     else
00483                     {
00484                         //Approach from the top for high percentiles
00485                         int sum=0;
00486                         int threshold = max(0, (int)floor(total * (1-ptile)+.5)- 1);
00487 
00488                         for(int j=255; j > 0; j--)
00489                         {
00490                             sum += histogram[j];
00491 
00492                             if(sum > threshold)
00493                                 return j;
00494                         }
00495                         
00496                         return 0;
00497                     }
00498                 }
00499         };
00500 
00501 
00504         template<> class Median<byte>: public Percentile<byte>
00505         {
00506             public:
00507                 Median()
00508                 :Percentile<byte>(0.5)
00509                 {
00510                 }
00511         };
00512 
00516         template<class T> struct BasicBinary
00517         {
00518                 int t, f;
00519                 BasicBinary()
00520                 :t(0),f(0)
00521                 {}
00522 
00523                 void insert(bool b)
00524                 {
00525                     if(b)
00526                         t++;
00527                     else
00528                         f++;
00529                 }
00530 
00531                 void remove(bool b)
00532                 {
00533                     if(b)
00534                         t--;
00535                     else
00536                         f--;
00537                 }
00538 
00539                 void clear()
00540                 {
00541                     t=f=0;
00542                 }
00543         };
00544 
00547         template<class T=bool> struct BinaryErode : public BasicBinary<T>
00548         {
00549             using BasicBinary<T>::f;
00550             T get()
00551             {
00552                 return !f;
00553             }
00554         };
00555 
00558         template<class T=bool> struct BinaryDilate : public BasicBinary<T>
00559         {
00560             using BasicBinary<T>::t;
00561             T get()
00562             {
00563                 return t!= 0;
00564             }
00565         };
00566 
00569         template<class T=bool> struct BinaryMedian : public BasicBinary<T>
00570         {
00571             using BasicBinary<T>::t;
00572             using BasicBinary<T>::f;
00573             T get()
00574             {
00575                 return (t > f);
00576             }
00577         };
00578     }
00579 
00580     namespace median{
00581         //Some helper classes for median
00582         template<class T> T median4(T a, T b, T c, T d)
00583         {
00584             int v[4] = {a, b, c, d};
00585             std::nth_element(v, v+2, v+4);
00586             return v[2];
00587         }
00588 
00589         template<class T> T median4(const SubImage<T>& im, int r, int c)
00590         {
00591             return median4(im[r][c], im[r][c+1], im[r+1][c], im[r+1][c+1]);
00592         }
00593 
00594         template<class T> T median6(T a, T b, T c, T d, T e, T f)
00595         {
00596             int v[6] = {a, b, c, d, e, f};
00597             std::nth_element(v, v+3, v+6);
00598             return v[3];
00599         }
00600 
00601         template<class T> T median6_row(const SubImage<T>& im, int r, int c)
00602         {
00603             return median6(im[r][c], im[r][c+1], im[r][c+2], im[r+1][c], im[r+1][c+1], im[r+1][c+2]);
00604         }
00605         template<class T> T median6_col(const SubImage<T>& im, int r, int c)
00606         {
00607             return median6(im[r][c], im[r][c+1], im[r+1][c], im[r+1][c+1], im[r+2][c], im[r+2][c+1]);
00608         }
00609 
00610     };
00611     
00612     #ifndef DOXYGEN_IGNORE_INTERNAL
00613         //Overload for median filtering of byte images, to special-case 3x3
00614         void morphology(const SubImage<byte>& in, const std::vector<ImageRef>& selem, const Morphology::Median<byte>& m, SubImage<byte>& out);
00615     #endif
00616 
00617 }
00618 
00619 #endif