CVD 0.8
|
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