00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef MLN_CLUSTERING_KMEAN_RGB_HH
00027 # define MLN_CLUSTERING_KMEAN_RGB_HH
00028
00093
00094 # include <limits.h>
00095 # include <iostream>
00096
00097 # include <mln/accu/stat/histo3d_rgb.hh>
00098
00099 # include <mln/algebra/vec.hh>
00100
00101 # include <mln/core/concept/image.hh>
00102 # include <mln/core/contract.hh>
00103 # include <mln/core/image/image2d.hh>
00104 # include <mln/core/macros.hh>
00105
00106 # include <mln/data/compute.hh>
00107 # include <mln/data/fill.hh>
00108 # include <mln/data/transform.hh>
00109
00110 # include <mln/debug/println.hh>
00111
00112 # include <mln/io/ppm/save.hh>
00113 # include <mln/io/pgm/save.hh>
00114
00115 # include <mln/labeling/colorize.hh>
00116 # include <mln/labeling/mean_values.hh>
00117
00118 # include <mln/literal/zero.hh>
00119 # include <mln/literal/one.hh>
00120
00121 # include <mln/math/min.hh>
00122 # include <mln/math/sqr.hh>
00123
00124 # include <mln/norm/l2.hh>
00125
00126 # include <mln/opt/at.hh>
00127
00128 # include <mln/trace/entering.hh>
00129 # include <mln/trace/exiting.hh>
00130
00131 # include <mln/trait/value_.hh>
00132
00133 # include <mln/util/array.hh>
00134
00135 # include <mln/value/int_u.hh>
00136 # include <mln/value/rgb8.hh>
00137 # include <mln/value/label_8.hh>
00138
00139
00140
00141
00142
00143
00144
00145 namespace mln
00146 {
00147
00148 namespace clustering
00149 {
00168 template <typename T, unsigned n, typename I>
00169 inline
00170 mln_ch_value(I,value::label_8)
00171 kmean_rgb(const Image<I>& point,
00172 const unsigned k_center,
00173 const unsigned watch_dog,
00174 const unsigned n_times);
00175
00176 }
00177
00178 namespace clustering
00179 {
00180
00181
00182 # ifndef MLN_INCLUDE_ONLY
00183
00184
00185
00186
00187
00188 namespace internal
00189 {
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483 }
00484
00485
00486
00487
00488
00489
00490 namespace impl
00491 {
00492
00493
00494
00495
00496
00497
00498
00499
00500 template <unsigned n>
00501 struct rgbn_to_lbl8 : Function_v2v< rgbn_to_lbl8<n> >
00502 {
00503 typedef value::rgb<n> argument;
00504 typedef value::label_8 result;
00505 typedef value::label_8 t_label;
00506 typedef image3d<t_label> t_group_img;
00507
00508 t_group_img _group;
00509
00510 rgbn_to_lbl8(t_group_img group) : _group(group) {}
00511
00512 result operator()(const argument& c) const
00513 {
00514 value::label_8 tmp = opt::at(_group, c.blue(), c.red(), c.green());
00515
00516
00517 return ++tmp;
00518 }
00519 };
00520
00521 template <typename T, unsigned n>
00522 struct rgb_to_dist : Function_v2v< rgb_to_dist<T,n> >
00523 {
00524 typedef value::rgb<n> argument;
00525 typedef T result;
00526 typedef T t_result1d;
00527 typedef algebra::vec<3,T> t_result3d;
00528 typedef image3d<unsigned> t_histo_img;
00529
00530 t_result3d _mean;
00531 t_histo_img _histo;
00532
00533 rgb_to_dist(t_result3d mean, t_histo_img histo) : _mean(mean),
00534 _histo(histo) {}
00535
00536 result operator()(const argument& c) const
00537 {
00538 t_result1d diff2_row = math::sqr(c.row() - _mean[0]);
00539 t_result1d diff2_col = math::sqr(c.col() - _mean[1]);
00540 t_result1d diff2_sli = math::sqr(c.sli() - _mean[2]);
00541 t_result1d tmp = _histo(c)*(diff2_row + diff2_col + diff2_sli);
00542
00543 return tmp;
00544 }
00545 };
00546
00547 template <typename T, unsigned n, typename I>
00548 inline
00549 mln_ch_value(I,value::label_8)
00550 kmean_image2d_rgb(const Image<I>& point__,
00551 const unsigned k_center,
00552 const unsigned watch_dog = 10,
00553 const unsigned n_times = 10)
00554 {
00555 trace::entering("mln::clustering::impl::kmean_image2d_rgb");
00556
00557 const I& point = exact(point__);
00558 typedef mln_value(I) V;
00559 mlc_is(V, value::rgb<n>)::check();
00560 mlc_bool(mln_site_(I)::dim == 2u)::check();
00561 mln_precondition(point.is_valid());
00562
00563
00564 typedef value::rgb<8> t_rgb;
00565 typedef value::label<8> t_label;
00566 typedef value::rgb<n> t_value;
00567 typedef mln_trait_value_comp(t_value,0) t_value_comp0;
00568 typedef mln_trait_value_comp(t_value,1) t_value_comp1;
00569 typedef mln_trait_value_comp(t_value,2) t_value_comp2;
00570 typedef T t_result1d;
00571 typedef algebra::vec<3,T> t_result3d;
00572
00573 typedef I t_point_img;
00574 typedef image3d<unsigned> t_histo_img;
00575 typedef util::array<t_result1d> t_number_img;
00576 typedef util::array<t_result3d> t_mean_img;
00577 typedef util::array<t_result1d> t_variance_img;
00578
00579 typedef image3d<t_label> t_group_img;
00580 typedef image3d<t_result1d> t_distance_val;
00581 typedef util::array<t_distance_val> t_distance_img;
00582
00583 typedef mln_ch_value(I,t_label) t_label_dbg;
00584 typedef image2d<t_rgb> t_color_dbg;
00585 typedef image2d<t_value> t_mean_dbg;
00586
00587 typedef image1d<t_result3d> t_mean_val;
00588 typedef util::array<t_mean_val> t_mean_set;
00589 typedef util::array<t_mean_set> t_mean_cnv;
00590 typedef image1d<t_result1d> t_variance_val;
00591 typedef util::array<t_variance_val> t_variance_cnv;
00592
00593
00594
00595 mln_precondition(point.is_valid());
00596
00597 static const unsigned _N_TRIES = 3;
00598
00599 typedef accu::meta::stat::histo3d_rgb t_histo3d_rgb;
00600
00601 t_result1d _within_variance;
00602
00603 unsigned _k_center = k_center;
00604 unsigned _watch_dog = watch_dog;
00605 unsigned _n_times = n_times;
00606 t_point_img _point = point;
00607
00608
00609 t_histo_img _histo = data::compute(t_histo3d_rgb(),
00610 _point);
00611
00612
00613 t_number_img _number;
00614 t_mean_img _mean;
00615 t_variance_img _variance;
00616
00617 for (unsigned i = 0; i < _k_center; ++i)
00618 {
00619 _number.append(literal::zero);
00620 _mean.append(literal::zero);
00621 _variance.append(literal::zero);
00622 }
00623
00624
00625 unsigned _current_step = 0;
00626 unsigned _current_launching = 0;
00627 bool _is_number_valid = false;
00628
00629 unsigned _launching_min;
00630 t_result1d _variance_min;
00631 t_mean_img _mean_min;
00632
00633
00634
00635 t_group_img _group;
00636 t_distance_img _distance;
00637
00638
00639 t_label_dbg _label_dbg;
00640 t_color_dbg _color_dbg;
00641 t_mean_dbg _mean_dbg;
00642
00643
00644 t_mean_cnv _mean_cnv;
00645 t_variance_cnv _variance_cnv;
00646
00647
00648
00649
00650 _group.init_(box3d(point3d(mln_min(t_value_comp2),
00651 mln_min(t_value_comp0),
00652 mln_min(t_value_comp1)),
00653 point3d(mln_max(t_value_comp2),
00654 mln_max(t_value_comp0),
00655 mln_max(t_value_comp1))));
00656
00657 for (unsigned i = 0; i < _k_center; ++i)
00658 {
00659 t_distance_val img(box3d(point3d(mln_min(t_value_comp2),
00660 mln_min(t_value_comp0),
00661 mln_min(t_value_comp1)),
00662 point3d(mln_max(t_value_comp2),
00663 mln_max(t_value_comp0),
00664 mln_max(t_value_comp1))));
00665
00666 _distance.append(img);
00667 }
00668
00669
00670 initialize(_label_dbg, _point);
00671 initialize(_color_dbg, _point);
00672 initialize(_mean_dbg, _point);
00673
00674
00675
00676 for (unsigned i = 0; i < _n_times; ++i)
00677 {
00678 t_variance_val img(box1d(point1d(0), point1d(_watch_dog-1)));
00679
00680 data::fill(img, literal::zero);
00681
00682 _variance_cnv.append(img);
00683 }
00684
00685 for (unsigned i = 0; i < _k_center; ++i)
00686 {
00687 t_mean_set mean_set;
00688
00689 for (unsigned j = 0; j < _n_times; ++j)
00690 {
00691 t_mean_val img(box1d(point1d(0), point1d(_watch_dog-1)));
00692
00693 data::fill(img, literal::zero);
00694
00695 mean_set.append(img);
00696 }
00697
00698 _mean_cnv.append(mean_set);
00699 }
00700
00701
00702
00703 {
00704 unsigned tries = 0;
00705 _variance_min = mln_max(t_result1d);
00706 _current_launching = 0;
00707
00708 while (_current_launching < _n_times)
00709 {
00710
00711 trace::entering("Launch one time");
00712 {
00713 t_result1d old_variance = mln_max(t_result1d);
00714 _within_variance = mln_max(t_result1d);
00715 _current_step = 0;
00716
00717
00718 trace::entering("init mean");
00719 {
00720 t_value_comp0 min_comp0 = mln_min(t_value_comp0);
00721 t_value_comp0 max_comp0 = mln_max(t_value_comp0);
00722 t_value_comp1 min_comp1 = mln_min(t_value_comp1);
00723 t_value_comp1 max_comp1 = mln_max(t_value_comp1);
00724 t_value_comp2 min_comp2 = mln_min(t_value_comp2);
00725 t_value_comp2 max_comp2 = mln_max(t_value_comp2);
00726 mln_eiter(t_mean_img) l(_mean);
00727
00728 for_all(l)
00729 {
00730 _mean[l.index_()][0]=(rand()%(max_comp0-min_comp0))+min_comp0;
00731 _mean[l.index_()][1]=(rand()%(max_comp1-min_comp1))+min_comp1;
00732 _mean[l.index_()][2]=(rand()%(max_comp2-min_comp2))+min_comp2;
00733 }
00734 }
00735 trace::exiting("init mean");
00736
00737
00738
00739
00740 trace::entering("update distance");
00741
00742 for (unsigned i = 0; i < _k_center; ++i)
00743 {
00744
00745
00746
00747
00748
00749 mln_piter(t_distance_val) d(_distance[i].domain());
00750
00751 for_all(d)
00752 {
00753 t_result1d diff2_row = math::sqr(d.row() - _mean[i][0]);
00754 t_result1d diff2_col = math::sqr(d.col() - _mean[i][1]);
00755 t_result1d diff2_sli = math::sqr(d.sli() - _mean[i][2]);
00756 _distance[i](d) = _histo(d)*
00757 (diff2_row + diff2_col + diff2_sli);
00758 }
00759 }
00760
00761 trace::exiting("update distance");
00762
00763
00764 do
00765 {
00766 old_variance = _within_variance;
00767
00768
00769 trace::entering("update group");
00770 {
00771 mln_piter(t_group_img) rgb(_group.domain());
00772
00773 for_all(rgb)
00774 {
00775 mln_eiter(t_distance_img) l(_distance);
00776 t_result1d min = mln_max(t_result1d);
00777 t_label label = mln_max(t_label);
00778
00779 for_all(l)
00780 {
00781 if (min > _distance[l.index_()](rgb))
00782 {
00783 min = _distance[l.index_()](rgb);
00784 label = l.index_();
00785 }
00786 }
00787
00788 _group(rgb) = label;
00789 }
00790
00791 }
00792 trace::exiting("update group");
00793
00794
00795
00796 trace::entering("update mean");
00797 {
00798 mln_eiter(t_number_img) en(_number);
00799 mln_eiter(t_mean_img) em(_mean);
00800
00801 for_all_2(en,em)
00802 {
00803 _number[en.index_()] = literal::zero;
00804 _mean[em.index_()] = literal::zero;
00805 }
00806
00807 mln_piter(t_group_img) rgb(_group.domain());
00808
00809 for_all(rgb)
00810 {
00811 _mean[_group(rgb)][0] += rgb.row() * _histo(rgb);
00812 _mean[_group(rgb)][1] += rgb.col() * _histo(rgb);
00813 _mean[_group(rgb)][2] += rgb.sli() * _histo(rgb);
00814 _number(_group(rgb)) += _histo(rgb);
00815 }
00816
00817 mln_eiter(t_mean_img) l(_mean);
00818
00819 for_all(l)
00820 {
00821 _is_number_valid = (0 != _number[l.index_()]);
00822
00823 if (!_is_number_valid)
00824 break;
00825
00826 _mean[l.index_()] /= _number[l.index_()];
00827 }
00828 }
00829 trace::exiting("update mean");
00830
00831
00832
00833
00834 if (!_is_number_valid)
00835 break;
00836
00837
00838 trace::entering("update distance");
00839
00840 for (unsigned i = 0; i < _k_center; ++i)
00841 {
00842 mln_piter(t_distance_val) d(_distance[i].domain());
00843
00844 for_all(d)
00845 {
00846
00847 t_result1d diff2_row = math::sqr(d.row() - _mean[i][0]);
00848 t_result1d diff2_col = math::sqr(d.col() - _mean[i][1]);
00849 t_result1d diff2_sli = math::sqr(d.sli() - _mean[i][2]);
00850 _distance[i](d) = _histo(d)*
00851 (diff2_row + diff2_col + diff2_sli);
00852 }
00853 }
00854 trace::exiting("update distance");
00855
00856
00857
00858 trace::entering("update variance");
00859 {
00860 _within_variance = literal::zero;
00861 mln_eiter(t_variance_img) l(_variance);
00862
00863 for_all(l)
00864 {
00865 _variance[l.index_()] = literal::zero;
00866
00867 mln_piter(t_group_img) rgb(_group.domain());
00868
00869 for_all(rgb)
00870 {
00871 if (l.index_() == _group(rgb))
00872 _variance[l.index_()] += _distance[l.index_()](rgb);
00873 }
00874
00875 _within_variance += _variance[l.index_()];
00876 }
00877
00878 }
00879 trace::exiting("update variance");
00880
00881
00882
00883
00884 ++_current_step;
00885 }
00886 while (_current_step < _watch_dog &&
00887 _within_variance < old_variance);
00888
00889
00890
00891 }
00892 trace::exiting("Launch one time");
00893
00894
00895 if ((_is_number_valid && (_current_step < _watch_dog))||
00896 _N_TRIES < tries)
00897 {
00898 if (_within_variance < _variance_min)
00899 {
00900 _variance_min = _within_variance;
00901 _mean_min = _mean;
00902 _launching_min = _current_launching;
00903 }
00904
00905
00906 tries = 0;
00907
00908
00909
00910
00911
00912
00913
00914 ++_current_launching;
00915 }
00916 else
00917 ++tries;
00918 }
00919
00920
00921
00922
00923 }
00924
00925
00926
00927 _label_dbg = data::transform(_point, rgbn_to_lbl8<n>(_group));
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943 trace::exiting("mln::clustering::impl::kmean_image2d_rgb");
00944
00945 return _label_dbg;
00946
00947 }
00948
00949 }
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959 namespace internal
00960 {
00961
00962 template <typename T, unsigned n, typename I>
00963 inline
00964 mln_ch_value(I,value::label_8)
00965 kmean_rgb_dispatch(const Image<I>& img,
00966 const unsigned k_center,
00967 const unsigned watch_dog,
00968 const unsigned n_times,
00969 const value::rgb<n>&,
00970 const point2d&)
00971 {
00972 return impl::kmean_image2d_rgb<T,n>(img, k_center, watch_dog, n_times);
00973 }
00974
00975 template <typename T, unsigned n, typename I, typename V, typename P>
00976 inline
00977 mln_ch_value(I,value::label_8)
00978 kmean_rgb_dispatch(const Image<I>& img,
00979 const unsigned k_center,
00980 const unsigned watch_dog,
00981 const unsigned n_times,
00982 const V&,
00983 const P&)
00984 {
00985
00986 mlc_abort(I)::check();
00987
00988 typedef mln_ch_value(I, value::label_8) output_t;
00989 return output_t();
00990 }
00991
00992
00993 template <typename T, unsigned n, typename I>
00994 inline
00995 mln_ch_value(I,value::label_8)
00996 kmean_rgb_dispatch(const Image<I>& img,
00997 const unsigned k_center,
00998 const unsigned watch_dog,
00999 const unsigned n_times)
01000 {
01001 typedef mln_value(I) V;
01002 typedef mln_site(I) P;
01003 return kmean_rgb_dispatch<T,n>(img, k_center, watch_dog,
01004 n_times, V(), P());
01005 }
01006
01007
01008 }
01009
01010
01011
01012
01013
01014
01015 template <typename T, unsigned n, typename I>
01016 inline
01017 mln_ch_value(I,value::label_8)
01018 kmean_rgb(const Image<I>& point,
01019 const unsigned k_center,
01020 const unsigned watch_dog,
01021 const unsigned n_times)
01022 {
01023 trace::entering("mln::clustering::kmean_rgb");
01024
01025 mln_ch_value(I, value::label_8)
01026 output = internal::kmean_rgb_dispatch<T,n>(point, k_center,
01027 watch_dog, n_times);
01028 trace::exiting("mln::clustering::kmean_rgb");
01029
01030 return output;
01031 }
01032
01033
01034 # endif // ! MLN_INCLUDE_ONLY
01035
01036 }
01037
01038 }
01039
01040 #endif // ! MLN_CLUSTERING_KMEAN_RGB_HH