* green/exp/annotating/error: New directory.
* green/exp/annotating/error/Makefile.am: New Makefile.
* green/exp/annotating/error/error.cc: New source.
---
.../annotating/{achromastism => error}/Makefile.am | 0
milena/sandbox/green/exp/annotating/error/error.cc | 700 ++++++++++++++++++++
2 files changed, 700 insertions(+), 0 deletions(-)
copy milena/sandbox/green/exp/annotating/{achromastism => error}/Makefile.am (100%)
create mode 100644 milena/sandbox/green/exp/annotating/error/error.cc
diff --git a/milena/sandbox/green/exp/annotating/achromastism/Makefile.am
b/milena/sandbox/green/exp/annotating/error/Makefile.am
similarity index 100%
copy from milena/sandbox/green/exp/annotating/achromastism/Makefile.am
copy to milena/sandbox/green/exp/annotating/error/Makefile.am
diff --git a/milena/sandbox/green/exp/annotating/error/error.cc
b/milena/sandbox/green/exp/annotating/error/error.cc
new file mode 100644
index 0000000..5bc67ba
--- /dev/null
+++ b/milena/sandbox/green/exp/annotating/error/error.cc
@@ -0,0 +1,700 @@
+// HSV TEST CF MILLET 2008
+
+#include <iostream>
+#include <sstream>
+#include <boost/filesystem.hpp>
+
+#include <mln/algebra/vec.hh>
+
+#include <mln/img_path.hh>
+
+#include <mln/accu/stat/mean.hh>
+#include <mln/accu/stat/histo1d.hh>
+
+#include <mln/arith/minus.hh>
+#include <mln/arith/times.hh>
+#include <mln/arith/diff_abs.hh>
+
+#include <mln/core/image/image1d.hh>
+#include <mln/core/image/image2d.hh>
+#include <mln/core/image/dmorph/image_if.hh>
+#include <mln/core/alias/point1d.hh>
+#include <mln/core/alias/box1d.hh>
+
+#include <mln/data/transform.hh>
+#include <mln/data/compute.hh>
+#include <mln/data/stretch.hh>
+#include <mln/data/fill.hh>
+
+#include <mln/fun/v2v/component.hh>
+
+#include <mln/io/ppm/load.hh>
+#include <mln/io/pgm/save.hh>
+#include <mln/io/plot/save_image_sh.hh>
+
+#include <mln/literal/zero.hh>
+
+#include <mln/math/ceil.hh>
+#include <mln/math/floor.hh>
+
+#include <mln/opt/at.hh>
+
+#include <mln/trait/value_.hh>
+
+#include <mln/value/rgb8.hh>
+#include <mln/value/int_u8.hh>
+
+//============================================================================//
+// CLASSIFICATION DE FISHER EN 2 CLASSES SUR UN HISTO 1D
+//============================================================================//
+
+template <typename I>
+mln_value(I) cnt_histo(const mln::Image<I>& histo_)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) sum = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ {
+ sum += histo(p);
+ }
+
+ return sum;
+}
+
+template <typename I>
+mln_value(I) sum_histo(const mln::Image<I>& histo_)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) pos = mln::literal::zero;
+ mln_value(I) sum = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ {
+ pos = p.ind();
+ sum += pos*histo(p);
+ }
+
+ return sum;
+}
+
+template <typename I>
+mln_value(I) avg_histo(const mln::Image<I>& histo_)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) pos = mln::literal::zero;
+ mln_value(I) sum = mln::literal::zero;
+ mln_value(I) cnt = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ {
+ pos = p.ind();
+ cnt += histo(p);
+ sum += pos*histo(p);
+ }
+
+ return (0 == cnt)? 0 : sum/cnt;
+}
+
+template <typename I>
+mln_value(I) var3_histo(const mln::Image<I>& histo_, float mean)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) pos = mln::literal::zero;
+ mln_value(I) sum = mln::literal::zero;
+ mln_value(I) cnt = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ {
+ cnt += histo(p);
+ sum += (mln::math::sqr(p.ind()-mean)*histo(p));
+ }
+
+ return (0 == cnt)? 0 : sum/cnt;
+}
+
+template <typename I>
+mln_value(I) var2_histo(const mln::Image<I>& histo_, float mean)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) pos = mln::literal::zero;
+ mln_value(I) sum = mln::literal::zero;
+ mln_value(I) sqr = mln::literal::zero;
+ mln_value(I) cnt = mln::literal::zero;
+ mln_value(I) dlt = mln::literal::zero;
+ mln_value(I) mxt = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ {
+ pos = p.ind();
+ cnt += (histo(p));
+ sum += (histo(p)*pos);
+ mxt += (histo(p)*pos*mean);
+ sqr += (histo(p)*mln::math::sqr(pos));
+ dlt += (histo(p)*mln::math::sqr(pos - mean));
+
+ std::cout << "p = " << pos
<< std::endl;
+ std::cout << "p² = " << mln::math::sqr(pos)
<< std::endl;
+ std::cout << "p*mean = " << (pos*mean)
<< std::endl;
+ std::cout << "p-mean = " << (pos-mean)
<< std::endl;
+ std::cout << "(p-mean)² = " << mln::math::sqr(pos-mean)
<< std::endl;
+ std::cout << "histo(p) = " << histo(p)
<< std::endl;
+ std::cout << "histo(p)*p = " << (pos*histo(p))
<< std::endl;
+ std::cout << "histo(p)*p²= " << (mln::math::sqr(pos)*histo(p))
+ << std::endl;
+ std::cout << "cnt = " << cnt << std::endl;
+ std::cout << "sum = " << sum << std::endl;
+ std::cout << "sqr = " << sqr << std::endl;
+ std::cout << "dlt = " << dlt << std::endl;
+ std::cout << "mxt = " << mxt << std::endl;
+ std::cout << std::endl;
+ }
+
+ std::cout << "sqr/cnt = " << (sqr/cnt)
<< std::endl;
+ std::cout << "sum/cnt = " << (sum/cnt)
<< std::endl;
+ std::cout << "(sum/cnt)² = " << mln::math::sqr(sum/cnt)
<< std::endl;
+ std::cout << "dlt/cnt = " << dlt/cnt
<< std::endl;
+ std::cout << "mxt/cnt = " << mxt/cnt
<< std::endl;
+ std::cout << std::endl;
+
+ std::cout << "sqr = "
+ << (sqr) << std::endl;
+ std::cout << "dlt = "
+ << (dlt) << std::endl;
+ std::cout << "cnt*avg² = "
+ << (mln::math::sqr(sum/cnt)*cnt) << std::endl;
+ std::cout << "2*mxt = "
+ << (2*mxt) << std::endl;
+ std::cout << "sqr - cnt*avg² = "
+ << (sqr/cnt - mln::math::sqr(sum/cnt)) << std::endl;
+ std::cout << "(sqr -2*mxt + cnt*avg²) = "
+ << ((sqr - 2*mxt + mln::math::sqr(sum/cnt))/cnt) << std::endl;
+ std::cout << std::endl;
+ return (0 == cnt)? 0 : sqr/cnt - mln::math::sqr(sum/cnt);
+}
+
+template <typename I>
+mln_value(I) var_histo(const mln::Image<I>& histo_)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) pos = mln::literal::zero;
+ mln_value(I) sum = mln::literal::zero;
+ mln_value(I) sqr = mln::literal::zero;
+ mln_value(I) cnt = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ {
+ pos = p.ind();
+ cnt += (histo(p));
+ sum += (histo(p)*pos);
+ sqr += (histo(p)*mln::math::sqr(pos));
+ }
+
+ return (0 == cnt)? 0 : sqr/cnt - mln::math::sqr(sum/cnt);
+}
+
+//============================================================================//
+// CLASSIFIEUR
+//============================================================================//
+
+
+// Linear discriminant analysis in 1d
+// With same variance, threshold = (m1+m2)/2
+// With different variance, (m1*sqrt(v1)+m2*sqrt(v2))/(sqrt(v1)+sqrt(v2))
+float threshold_histo(float avg1, float var1, float avg2, float var2)
+{
+ float sigma1 = mln::math::sqrt(var1);
+ float sigma2 = mln::math::sqrt(var2);
+ float threshold = (avg1*sigma1+avg2*sigma2)/(sigma1+sigma2);
+
+ return threshold;
+}
+
+float threshold3_histo(float avg1, float var1, float avg2, float var2)
+{
+ float threshold = (avg1*var1+avg2*var2)/(var1+var2);
+
+ return threshold;
+}
+
+
+// if gaussian without the same variance
+float threshold2_histo(float avg1, float var1, float avg2, float var2)
+{
+ float a = var2 - var1;
+ float b = -2 * (avg1 * var2 - avg2 * var1);
+ float c = var2 * mln::math::sqr(avg1) - var1 * mln::math::sqr(avg2);
+ float d = mln::math::sqr(b) - 4 * a * c;
+
+ if (d < 0)
+ std::cout << "delta negatif" << std::endl;
+
+ float x1 = (-b+mln::math::sqrt(d))/(2*a);
+ float x2 = (-b-mln::math::sqrt(d))/(2*a);
+
+ std::cout << "a = " << a << std::endl;
+ std::cout << "b = " << b << std::endl;
+ std::cout << "c = " << c << std::endl;
+ std::cout << "d = " << d << std::endl;
+ std::cout << "x1 = " << x1 << std::endl;
+ std::cout << "x2 = " << x2 << std::endl;
+
+ return x2;
+}
+
+template <typename I>
+mln_value(I) sqr_histo(const mln::Image<I>& histo_)
+{
+ const I& histo = exact(histo_);
+
+ mln_precondition(histo.is_valid());
+
+ mln_value(I) sum = mln::literal::zero;
+ mln_piter(I) p(histo.domain());
+
+ for_all(p)
+ sum += (mln::math::sqr(p.ind())*histo(p));
+
+ return sum;
+}
+
+
+short min_error(const mln::image1d<float> histo_grp1,
+ const mln::image1d<float> histo_grp2)
+{
+ float c00[256];
+ float c10[256];
+ float c01[256];
+ float c11[256];
+ float err[256];
+
+ for (short p = 0; p < 256; p++)
+ {
+ c00[p] = cnt_histo(histo_grp1|mln::box1d(mln::point1d(0),
+ mln::point1d(p)));
+
+ c10[p] = cnt_histo(histo_grp1|mln::box1d(mln::point1d(p+1),
+ mln::point1d(255)));
+
+ c01[p] = cnt_histo(histo_grp2|mln::box1d(mln::point1d(0),
+ mln::point1d(p)));
+
+ c11[p] = cnt_histo(histo_grp2|mln::box1d(mln::point1d(p+1),
+ mln::point1d(255)));
+ }
+
+ short threshold = 0;
+ float error = c01[0] + c01[0] + c00[0] + c11[0];
+
+ for(short p = 0; p < 256; p++)
+ {
+ err[p] = c10[p] + c01[p];
+
+ std::cout << " p = " << p
+ << ";c00 = " << c00[p]
+ << ";c10 = " << c10[p]
+ << ";c01 = " << c01[p]
+ << ";c11 = " << c11[p]
+ << std::endl;
+// std::cout << "err[" << p << "] = " <<
err[p] << std::endl;
+
+ if (error > err[p])
+ {
+ error = err[p];
+ threshold = p;
+ }
+ }
+
+ return threshold;
+}
+
+// return the threshold
+short fisher_analysis(const mln::image1d<float> histo)
+{
+ typedef mln::value::int_u8 t_int_u8;
+
+ // FIXE ME SIZE const short a = mln_min(t_int_u8);
+ // float cnt1[a];
+
+ float cnt1[256];
+ float sum1[256];
+ float sqr1[256];
+ float avg1[256];
+ float var1[256];
+
+ float cnt2[256];
+ float sum2[256];
+ float sqr2[256];
+ float avg2[256];
+ float var2[256];
+
+ float cnt0[256]; // count of points
+ float sum0[256]; // sum of population
+ float sqr0[256]; // sqr of population
+ float avg0[256]; // average of the population
+ float v_in[256]; // variance with-in class
+ float v_bw[256]; // variance between class
+ float var0[256]; // variance of the population
+ short threshold;
+ float pos;
+ float min_var;
+
+ // Assign the first elements
+ cnt1[0] = 0;
+ sum1[0] = 0;
+ sqr1[0] = 0;
+ avg1[0] = 0;
+ var1[0] = 0;
+
+ // Compute the stats of the wall histogram
+ cnt2[0] = 0;
+ sum2[0] = 0;
+ sqr2[0] = 0;
+
+ for(short p = 0; p < 256; p++)
+ {
+ pos = p;
+ cnt2[0] += mln::opt::at(histo,p);
+ sum2[0] += (pos * mln::opt::at(histo,p));
+ sqr2[0] += (mln::math::sqr(pos) * mln::opt::at(histo,p));
+ }
+
+ avg2[0] = (0 == cnt2[0])? 0 : sum2[0] / cnt2[0];
+ var2[0] = (0 == cnt2[0])? 0 : sqr2[0] / cnt2[0] - mln::math::sqr(avg2[0]);
+
+ // watch the array limits
+ for (short p = 1; p < 256; p++)
+ {
+ pos = p;
+
+ // Assign the statistics to the primary class
+ cnt1[p] = cnt1[p-1] + mln::opt::at(histo, p);
+ sum1[p] = sum1[p-1] + pos * mln::opt::at(histo, p);
+ sqr1[p] = sqr1[p-1] + mln::math::sqr(pos) * mln::opt::at(histo, p);
+ avg1[p] = (0 == cnt1[p])? 0 : (sum1[p] / cnt1[p]);
+ var1[p] = (0 == cnt1[p])? 0 : ((sqr1[p] / cnt1[p])-mln::math::sqr(avg1[p]));
+
+ // Assign the statistics to the second class
+ cnt2[p] = cnt2[p-1] - mln::opt::at(histo, p);;
+ sum2[p] = sum2[p-1] - p * mln::opt::at(histo, p);;
+ sqr2[p] = sqr2[p-1] - mln::math::sqr(p) * mln::opt::at(histo, p);;
+ avg2[p] = (0 == cnt2[p])? 0 : (sum2[p] / cnt2[p]);
+ var2[p] = (0 == cnt2[p])? 0 : ((sqr2[p] / cnt2[p])-mln::math::sqr(avg2[p]));
+
+ // Lets compute the invariants
+ cnt0[p] = cnt1[p] + cnt2[p];
+ sum0[p] = sum1[p] + sum2[p];
+ sqr0[p] = sqr1[p] + sqr2[p];
+ avg0[p] = (cnt1[p] * avg1[p] + cnt2[p] * avg2[p])/cnt0[p];
+ v_in[p] = (cnt1[p] * var1[p] + cnt2[p] * var2[p])/cnt0[p];
+ v_bw[p] = (cnt1[p] * mln::math::sqr(avg1[p]-avg0[p]) +
+ cnt2[p] * mln::math::sqr(avg2[p]-avg0[p]))/cnt0[p];
+ var0[p] = v_in[p] + v_bw[p];
+ }
+
+ // Find the threshold that minimizes the intra-class variance
+ min_var = cnt2[0]*var2[0];
+ threshold = 0;
+
+ for(short p = 0; p < 256; p++)
+ {
+ // Compute the intra-class variance
+ v_in[p] = cnt1[p]*var1[p] + cnt2[p]*var2[p];
+// std::cout << "var intra[" << p << "]= "
<< v_in[p] << std::endl;
+
+ if (min_var > v_in[p])
+ {
+ min_var = v_in[p];
+ threshold = p;
+ }
+ }
+
+ return threshold;
+}
+
+
+//============================================================================//
+// ERROR (MSE, PNSNR) compression p278 Handbook Color
+//============================================================================//
+
+float error_test(const std::string original,
+ const std::string reduced)
+
+{
+
+ typedef mln::value::rgb8 t_rgb8;
+ typedef mln::value::int_u8 t_int_u8;
+ typedef mln::fun::v2v::component<t_rgb8,0> t_component_red;
+ typedef mln::fun::v2v::component<t_rgb8,1> t_component_green;
+ typedef mln::fun::v2v::component<t_rgb8,2> t_component_blue;
+ typedef mln::accu::meta::stat::mean t_mean;
+ typedef mln::accu::meta::stat::histo1d t_histo;
+ typedef mln::image2d<t_int_u8> t_img;
+ typedef mln_trait_op_minus_(t_img,t_img) t_minus;
+ typedef mln_trait_op_times_(t_minus,t_minus) t_times;
+
+ mln::image2d<mln::value::rgb8> original_rgb8;
+ mln::image2d<mln::value::rgb8> reduced_rgb8;
+
+ mln::image2d<mln::value::int_u8> original_red;
+ mln::image2d<mln::value::int_u8> original_green;
+ mln::image2d<mln::value::int_u8> original_blue;
+
+ mln::image2d<mln::value::int_u8> reduced_red;
+ mln::image2d<mln::value::int_u8> reduced_green;
+ mln::image2d<mln::value::int_u8> reduced_blue;
+
+// mln::image2d<mln::value::int_u8> map_red;
+// mln::image2d<mln::value::int_u8> map_green;
+// mln::image2d<mln::value::int_u8> map_blue;
+
+// mln::image1d<unsigned> histo_red;
+// mln::image1d<unsigned> histo_green;
+// mln::image1d<unsigned> histo_blue;
+
+ t_minus minus_red;
+ t_minus minus_green;
+ t_minus minus_blue;
+
+ t_times times_red;
+ t_times times_green;
+ t_times times_blue;
+
+ float error_red;
+ float error_green;
+ float error_blue;
+
+ float error;
+
+
+ mln::io::ppm::load(original_rgb8, original.c_str());
+ mln::io::ppm::load(reduced_rgb8, reduced.c_str());
+
+ original_red = mln::data::transform(original_rgb8, t_component_red());
+ original_green = mln::data::transform(original_rgb8, t_component_green());
+ original_blue = mln::data::transform(original_rgb8, t_component_blue());
+
+ reduced_red = mln::data::transform(reduced_rgb8, t_component_red());
+ reduced_green = mln::data::transform(reduced_rgb8, t_component_green());
+ reduced_blue = mln::data::transform(reduced_rgb8, t_component_blue());
+
+ minus_red = (original_red - reduced_red);
+ times_red = minus_red * minus_red;
+
+ minus_green = (original_green - reduced_green);
+ times_green = minus_green * minus_green;
+
+ minus_blue = (original_blue - reduced_blue);
+ times_blue = minus_blue * minus_blue;
+
+ error_red = mln::data::compute(t_mean(), times_red);
+ error_green = mln::data::compute(t_mean(), times_green);
+ error_blue = mln::data::compute(t_mean(), times_blue);
+
+// map_red = mln::data::stretch(t_int_u8(), times_red);
+// map_green = mln::data::stretch(t_int_u8(), times_blue);
+// map_blue = mln::data::stretch(t_int_u8(), times_green);
+
+// histo_red = mln::data::compute(t_histo(), map_red);
+// histo_green = mln::data::compute(t_histo(), map_green);
+// histo_blue = mln::data::compute(t_histo(), map_blue);
+
+// mln::io::plot::save_image_sh(histo_red, "histo_red.sh");
+// mln::io::plot::save_image_sh(histo_green,"histo_green.sh");
+// mln::io::plot::save_image_sh(histo_blue, "histo_blue.sh");
+
+// mln::io::pgm::save(map_red, "red.pgm");
+// mln::io::pgm::save(map_green,"green.pgm");
+// mln::io::pgm::save(map_blue, "blue.pgm");
+
+ error = (error_red + error_green + error_blue)/3.0;
+ error = mln::math::sqrt(error);
+ error = 20 * log(255/error);
+
+// Le PNSNR semble offrir plus d'espace pour la discrimination
+// Si les images sont identiques ==> PNSNR = +inf
+// Si les images sont très différentes ==> PNSNR = 0
+
+ return error;
+}
+
+
+//============================================================================//
+// MAIN
+//============================================================================//
+
+
+int main2()
+{
+ float val = error_test(AFP_PPM_IMG_PATH"/000_APP2003011515775.ppm",
+ AFP_GMP10_IMG_PATH"/000_APP2003011515775.ppm");
+
+ std::cout << val << std::endl;
+
+ return 0;
+}
+
+int main()
+{
+ typedef boost::filesystem::path t_path;
+ typedef boost::filesystem::directory_iterator t_iter_path;
+
+ mln::image1d<float> histo(256);
+ mln::image1d<float> histo_grp[2]; // histo by group
+
+ histo_grp[0].init_(mln::box1d(mln::point1d(0),mln::point1d(255)));
+ histo_grp[1].init_(mln::box1d(mln::point1d(0),mln::point1d(255)));
+
+ mln::data::fill(histo, mln::literal::zero);
+ mln::data::fill(histo_grp[0], mln::literal::zero);
+ mln::data::fill(histo_grp[1], mln::literal::zero);
+
+ t_path original_path[] = {ICDAR_20P_INPUT_IMG_PATH,
+ AFP_PPM_IMG_PATH};
+
+// t_path reduced1_path[] = {ICDAR_20P_MGK30_IMG_PATH,
+// AFP_MGK30_IMG_PATH};
+
+// t_path reduced1_path[] = {ICDAR_20P_MGK20_IMG_PATH,
+// AFP_MGK20_IMG_PATH};
+
+ t_path reduced1_path[] = {ICDAR_20P_MGK10_IMG_PATH,
+ AFP_MGK10_IMG_PATH};
+
+// t_path reduced2_path[] = {ICDAR_20P_GMP30_IMG_PATH,
+// AFP_GMP30_IMG_PATH};
+
+// t_path reduced2_path[] = {ICDAR_20P_GMP20_IMG_PATH,
+// AFP_GMP20_IMG_PATH};
+
+ t_path reduced2_path[] = {ICDAR_20P_GMP10_IMG_PATH,
+ AFP_GMP10_IMG_PATH};
+
+
+ std::cout << "#!/usr/bin/gnuplot" <<
std::endl;
+ std::cout << "set terminal x11 persist 1" <<
std::endl;
+ std::cout << "ERROR" <<
std::endl;
+ std::cout << "plot '-' using 1 with point notitle,\\"
<< std::endl;
+ std::cout << " '-' using 1 with point notitle"
<< std::endl;
+
+ for (int i = 0; i < 2; i++)
+ {
+ if (boost::filesystem::exists(original_path[i]) &&
+ boost::filesystem::exists(reduced1_path[i]) &&
+ boost::filesystem::exists(reduced2_path[i]) &&
+ boost::filesystem::is_directory(original_path[i]) &&
+ boost::filesystem::is_directory(reduced1_path[i]) &&
+ boost::filesystem::is_directory(reduced2_path[i]))
+ {
+ boost::filesystem::system_complete(original_path[i]);
+ const t_iter_path end_iter;
+ float error1 = 0.0;
+ float error2 = 0.0;
+ t_path leaf;
+ t_path reduced1_file;
+ t_path reduced2_file;
+
+ for (t_iter_path dir_iter(original_path[i]);end_iter!=dir_iter;++dir_iter)
+ {
+ leaf = dir_iter->path().leaf();
+ reduced1_file = reduced1_path[i] / leaf;
+ reduced2_file = reduced2_path[i] / leaf;
+
+ error1 = error_test(dir_iter->path().string(), reduced1_file.string());
+ error2 = error_test(dir_iter->path().string(), reduced2_file.string());
+
+ float a1 = 1;
+ short v1 = (short)mln::math::floor(error2+0.4999);
+ mln::opt::at(histo,v1) += a1;
+ mln::opt::at(histo_grp[i],v1) += a1;
+
+// float a1 = error2 - mln::math::floor(error2);
+// float a2 = mln::math::ceil(error2) - error2;
+// short v1 = (short)mln::math::floor(error2);
+// short v2 = (short)mln::math::ceil(error2);
+// mln::opt::at(histo,v1) += a1;
+// mln::opt::at(histo,v2) += a2;
+// mln::opt::at(histo_grp[i],v1) += a1;
+// mln::opt::at(histo_grp[i],v2) += a2;
+
+ std::cout << error1 << " ";
+ std::cout << error2 << " ";
+ std::cout << "# " << dir_iter->path().leaf() <<
std::endl;
+ }
+ std::cout << "e" << std::endl;
+ }
+ }
+
+ mln::io::plot::save_image_sh(histo, "histo.sh");
+ mln::io::plot::save_image_sh(histo_grp[1], "histo_grp1.sh");
+ mln::io::plot::save_image_sh(histo_grp[0], "histo_grp2.sh");
+
+ float threshold = fisher_analysis(histo);
+ float threshold2 = threshold_histo(avg_histo(histo_grp[1]),
+ var_histo(histo_grp[1]),
+ avg_histo(histo_grp[0]),
+ var_histo(histo_grp[0]));
+ float threshold3 = threshold2_histo(avg_histo(histo_grp[1]),
+ var_histo(histo_grp[1]),
+ avg_histo(histo_grp[0]),
+ var_histo(histo_grp[0]));
+ float threshold4 = min_error(histo_grp[1],histo_grp[0]);
+
+ std::cout << "threshold = " << threshold << std::endl;
+ std::cout << "threshold2 = " << threshold2 << std::endl;
+ std::cout << "threshold3 = " << threshold3 << std::endl;
+ std::cout << "threshold4 = " << threshold4 << std::endl;
+ std::cout << "avg_grp1 = " << avg_histo(histo_grp[1]) <<
std::endl;
+ std::cout << "avg_grp2 = " << avg_histo(histo_grp[0]) <<
std::endl;
+
+ // compute the classification matrix
+ // for each sub population
+
+ float c00 = cnt_histo(histo_grp[1] | mln::box1d(mln::point1d(0),
+ mln::point1d(threshold)));
+
+ float c10 = cnt_histo(histo_grp[1] | mln::box1d(mln::point1d(threshold+1),
+ mln::point1d(255)));
+
+ float c01 = cnt_histo(histo_grp[0] | mln::box1d(mln::point1d(0),
+ mln::point1d(threshold)));
+
+ float c11 = cnt_histo(histo_grp[0] | mln::box1d(mln::point1d(threshold+1),
+ mln::point1d(255)));
+
+
+ std::cout << "pop0 = " << cnt_histo(histo_grp[1]) <<
std::endl;
+ std::cout << "pop1 = " << cnt_histo(histo_grp[0]) <<
std::endl;
+ std::cout << std::endl;
+
+ std::cout << "c00 = " << c00 << std::endl;
+ std::cout << "c10 = " << c10 << std::endl;
+ std::cout << "c01 = " << c01 << std::endl;
+ std::cout << "c11 = " << c11 << std::endl;
+
+ return 0;
+}
--
1.5.6.5