[Olena-patches] r3471: Create small tools for IGR

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-03 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Create small tools for IGR. * fabien/TODO: New. * fabien/igr/dump2pbm.cc: Move this... * fabien/bin/dump2pbm.cc: ...here. * fabien/igr/dump2ppm.cc: Move this... * fabien/bin/dump2ppm.cc: ...here. * fabien/igr/Makefile: Update. * fabien/igr/launch.sh: New tool. * fabien/igr/launch2d.sh: New tool. * fabien/igr/launch3d.sh: New tool. * fabien/igr/seg2d.cc: New tool. * fabien/igr/seg3d.cc: New tool. * fabien/igr/seg_vol_irm.cc: Rename this... * fabien/igr/seg_vol_irm.hh: ...into this. --- TODO | 13 ++ bin/dump2pbm.cc | 35 ++++++ bin/dump2ppm.cc | 35 ++++++ igr/Makefile | 23 ++- igr/launch.sh | 4 igr/launch2d.sh | 26 ++++ igr/launch3d.sh | 24 ++++ igr/seg2d.cc | 40 ++++++ igr/seg3d.cc | 42 +++++++ igr/seg_vol_irm.hh | 309 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 10 files changed, 542 insertions(+), 9 deletions(-) Index: trunk/milena/sandbox/fabien/igr/dump2ppm.cc (deleted) =================================================================== Index: trunk/milena/sandbox/fabien/igr/seg_vol_irm.cc (deleted) =================================================================== Index: trunk/milena/sandbox/fabien/igr/dump2pbm.cc (deleted) =================================================================== Index: trunk/milena/sandbox/fabien/igr/launch.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/launch.sh (revision 0) +++ trunk/milena/sandbox/fabien/igr/launch.sh (revision 3471) @@ -0,0 +1,4 @@ +#!/bin/zsh + +./launch2d.sh +./launch3d.sh Property changes on: trunk/milena/sandbox/fabien/igr/launch.sh ___________________________________________________________________ Name: svn:executable + * Index: trunk/milena/sandbox/fabien/igr/seg_vol_irm.hh =================================================================== --- trunk/milena/sandbox/fabien/igr/seg_vol_irm.hh (revision 0) +++ trunk/milena/sandbox/fabien/igr/seg_vol_irm.hh (revision 3471) @@ -0,0 +1,309 @@ +// Copyright (C) 2007, 2008, 2009 EPITA Research and Development +// Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +#include <mln/core/concept/image.hh> +#include <mln/core/concept/neighborhood.hh> +#include <mln/core/alias/neighb1d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/core/alias/neighb3d.hh> + +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> + +#include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> +#include <mln/value/int_u16.hh> +#include <mln/value/rgb8.hh> + +#include <mln/io/pgm/all.hh> +#include <mln/io/pbm/all.hh> +#include <mln/io/ppm/all.hh> +#include <mln/io/cloud/all.hh> +#include <mln/io/dump/all.hh> +#include <mln/io/dicom/load.hh> + +#include <mln/labeling/flat_zones.hh> +#include <mln/labeling/background.hh> +#include <mln/literal/colors.hh> +#include <mln/norm/l1.hh> + +#include <mln/geom/bbox.hh> + +#include <mln/labeling/blobs.hh> +#include <mln/labeling/compute.hh> +#include <mln/labeling/fill_holes.hh> +#include <mln/labeling/n_max.hh> + +#include <mln/level/compare.hh> +#include <mln/level/compute.hh> +#include <mln/level/convert.hh> +#include <mln/level/stretch.hh> +#include <mln/level/transform.hh> + +#include <mln/fun/internal/selector.hh> + +#include <mln/fun/v2b/threshold.hh> +#include <mln/level/transform.hh> + +#include <mln/accu/count.hh> +#include <mln/accu/center.hh> +#include <mln/accu/max.hh> +#include <mln/accu/sum.hh> +#include <mln/accu/mean.hh> +#include <mln/accu/stat/deviation.hh> + +#include <mln/histo/compute.hh> + +#include <mln/set/compute.hh> +#include <mln/value/label_16.hh> +#include <mln/data/fill.hh> +#include <mln/pw/all.hh> + +#include <mln/morpho/elementary/gradient_internal.hh> +#include <mln/morpho/closing.hh> +#include <mln/morpho/dilation.hh> +#include <mln/morpho/erosion.hh> + +#include <mln/win/disk2d.hh> +#include <mln/win/sphere3d.hh> + +#include <mln/math/diff_abs.hh> + +#include <mln/convert/from_to.hh> + +#include <mln/metal/int.hh> + +#include <iostream> +#include <fstream> +#include <mln/debug/println.hh> + +using namespace mln; +using value::int_u8; +using value::int_u12; +using value::int_u16; +using value::rgb8; +using value::label_16; + + + +template <typename T> +struct L_to_int_u8 +: mln::fun::internal::selector_<int_u8, T, L_to_int_u8<T> >::ret +{ + typedef int_u8 result; + int_u8 operator()(const T& t) const; +}; + + + +template <typename T> +inline +int_u8 +L_to_int_u8<T>::operator()(const T& t) const +{ + return static_cast<int_u8>(t == 0 ? 0 : 1 + ((unsigned) t - 1) % 255); +} + + + +template <typename I> +I +close_threshold(const Image<I>& input, metal::int_<2>, int dil, int ero) +{ + return morpho::erosion(morpho::dilation(input, win::disk2d(dil)), win::disk2d(ero)); +} + + + +template <typename I> +I +close_threshold(const Image<I>& input, metal::int_<3>, int dil, int ero) +{ + return morpho::erosion(morpho::dilation(input, win::sphere3d(dil)), win::sphere3d(ero)); +} + + + +template <typename I, typename N> +unsigned +find_threshold_value(const Image<I>& input, const Neighborhood<N>& nbh) +{ + int bg_thres = 30; + int obj_thres = 15; + + mln_ch_value(I, bool) ima_bg; + initialize(ima_bg, input); + data::fill(ima_bg, false); + + mln_ch_value(I, bool) ima_obj; + initialize(ima_obj, input); + data::fill(ima_obj, false); + + unsigned result = 0; + + histo::array<mln_value(I)> arr_histo = histo::compute(input); + image1d<unsigned> ima_histo; + convert::from_to(arr_histo, ima_histo); + + // We remove the 0 value because it is not part of the image. + ima_histo(point1d(0)) = 0; + + for (unsigned int i = 1; i < ima_histo.nelements(); ++i) + { + ima_histo(point1d(i)) += ima_histo(point1d(i - 1)); + } + accu::max<unsigned> max_accu; + unsigned max = level::compute(max_accu, ima_histo); + bool low_done = false; + bool high_done = false; + for (unsigned int i = 0; i < ima_histo.nelements(); ++i) + { + if (!low_done && ((ima_histo(point1d(i)) * 100) / max) > bg_thres) + { + data::fill((ima_bg | pw::value(input) < pw::cst(i)).rw(), true); + low_done = true; + } + if (!high_done && ((ima_histo(point1d(i)) * 100) / max) > (100 - obj_thres)) + { + data::fill((ima_obj | pw::value(input) > pw::cst(i)).rw(), true); + high_done = true; + } + } + + // Debug output images + if (I::site::dim == 2) + { + io::pbm::save(ima_bg, "bg.pbm"); + io::pbm::save(ima_obj, "obj.pbm"); + } + if (I::site::dim == 3) + { + io::dump::save(ima_bg, "bg.dump"); + io::dump::save(ima_obj, "obj.dump"); + } + + ima_bg = close_threshold(ima_bg, metal::int_<I::site::dim>(), 3, 5); // 5, 7? + ima_obj = close_threshold(ima_obj, metal::int_<I::site::dim>(), 3, 5); + + // Debug output images + mln_ch_value(I, rgb8) out = level::convert(rgb8(), level::stretch(int_u8(), input)); + data::fill((out | pw::value(morpho::elementary::gradient_internal(ima_bg, nbh)) == true).rw(), literal::red); + data::fill((out | pw::value(morpho::elementary::gradient_internal(ima_obj, nbh)) == true).rw(), literal::green); + if (I::site::dim == 2) + { + io::pbm::save(ima_bg, "bg_closed.pbm"); + io::pbm::save(ima_obj, "obj_closed.pbm"); + io::ppm::save(out, "regions_color.ppm"); + } + if (I::site::dim == 3) + { + io::dump::save(ima_bg, "bg_closed.dump"); + io::dump::save(ima_obj, "obj_closed.dump"); + io::dump::save(out, "regions_color.dump"); + } + + histo::array<mln_value(I)> bg_histo = histo::compute(input | pw::value(ima_bg) == true); + histo::array<mln_value(I)> obj_histo = histo::compute(input | pw::value(ima_obj) == true); + + accu::sum<unsigned> sum_accu; + image1d<unsigned> ima_bg_histo; + convert::from_to(bg_histo, ima_bg_histo); + ima_bg_histo(point1d(0)) = 0; + unsigned bg_sum = level::compute(sum_accu, ima_bg_histo); + // We remove the 0 value because it is not part of the image. + std::ofstream fout_bg("bg_histo_norm.plot"); + for (unsigned int i = 0; i < ima_bg_histo.nelements(); ++i) + { + ima_bg_histo(point1d(i)) *= 10000; + ima_bg_histo(point1d(i)) /= bg_sum; + fout_bg << i << " " << ima_bg_histo(point1d(i)) << std::endl; + } + + image1d<unsigned> ima_obj_histo; + convert::from_to(obj_histo, ima_obj_histo); + unsigned obj_sum = level::compute(sum_accu, ima_obj_histo); + std::ofstream fout_obj("obj_histo_norm.plot"); + for (unsigned int i = 0; i < ima_obj_histo.nelements(); ++i) + { + ima_obj_histo(point1d(i)) *= 10000; + ima_obj_histo(point1d(i)) /= obj_sum; + fout_obj << i << " " << ima_obj_histo(point1d(i)) << std::endl; + } + + // Search for the index with the min distance between histogrammes. + unsigned min = math::diff_abs<unsigned>(ima_bg_histo(point1d(1)), ima_obj_histo(point1d(1))); + for (unsigned int i = 1; i < ima_bg_histo.nelements(); ++i) + { + if (math::diff_abs<unsigned>(ima_bg_histo(point1d(i)), ima_obj_histo(point1d(i))) < min) + { + min = math::diff_abs<unsigned>(ima_bg_histo(point1d(i)), ima_obj_histo(point1d(i))); + result = i; + } + } + + return result; +} + + + +template <typename I, typename N> +unsigned +find_threshold_mean(const Image<I>& input, const Neighborhood<N>& nbh) +{ + unsigned result = 0; + + accu::mean<unsigned> mean_accu; + unsigned mean = level::compute(mean_accu, (input | (pw::value(input) != 0))); + + accu::stat::deviation<unsigned, unsigned, float> dev_accu(mean); + float deviation = level::compute(dev_accu, (input | pw::value(input) != 0)); + + std::cout << "mean = " << mean << std::endl << "deviation = " << deviation << std::endl; + return floor(deviation); +} + + + +template <typename I, typename N, typename L> +mln_ch_value(I, bool) +igr_seg(const Image<I>& input_, const Neighborhood<N>& nbh_, L& nlabels) +{ + const I& input = exact(input_); + const N& nbh = exact(nbh_); + + // Threshold. + + mln_ch_value(I, bool) ima_thres; + initialize(ima_thres, input); + data::fill(ima_thres, false); + //unsigned threshold_value = find_threshold_value(input, nbh); + unsigned threshold_value = find_threshold_mean(input, nbh); + data::fill((ima_thres | pw::value(input) < pw::cst(threshold_value)).rw(), true); + + return ima_thres; +} Index: trunk/milena/sandbox/fabien/igr/seg2d.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/seg2d.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/seg2d.cc (revision 3471) @@ -0,0 +1,40 @@ +#include <mln/core/concept/image.hh> +#include <mln/core/concept/neighborhood.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/core/image/image2d.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/pbm/all.hh> +#include <mln/io/dicom/load.hh> + +#include <mln/value/label_16.hh> + +#include "seg_vol_irm.hh" + +using namespace mln; +using value::int_u12; +using value::label_16; + + + +int usage() +{ + std::cout << "Usage: ./seg2d image.dcm" << std::endl; + + return 1; +} + +int main(int argc, char* argv[]) +{ + if (argc < 2) + return usage(); + + label_16 nlabels; + image2d<int_u12> ima; + io::dicom::load(ima, argv[1]); + io::pbm::save(igr_seg(ima, c4(), nlabels), "result.pbm"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/launch2d.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/launch2d.sh (revision 0) +++ trunk/milena/sandbox/fabien/igr/launch2d.sh (revision 3471) @@ -0,0 +1,26 @@ +#!/bin/zsh + +process_file () +{ + ./seg2d $1 + if [ -f result.pbm ]; then + mv result.pbm results/${2}_06_result.pbm + fi + + if [ -f bg.pbm ]; then + mv bg.pbm results/${2}_01_bg.pbm + mv obj.pbm results/${2}_03_obj.pbm + + mv bg_closed.pbm results/${2}_02_bg_closed.pbm + mv obj_closed.pbm results/${2}_04_obj_closed.pbm + + mv regions_color.ppm results/${2}_05_regions_color.ppm + + mv bg_histo_norm.plot results/${2}_bg.plot + mv obj_histo_norm.plot results/${2}_obj.plot + fi +} + +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0049.dcm" "49" +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0055.dcm" "55" +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0058.dcm" "58" Property changes on: trunk/milena/sandbox/fabien/igr/launch2d.sh ___________________________________________________________________ Name: svn:executable + * Index: trunk/milena/sandbox/fabien/igr/seg3d.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/seg3d.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/seg3d.cc (revision 3471) @@ -0,0 +1,42 @@ +#include <mln/core/concept/image.hh> +#include <mln/core/concept/neighborhood.hh> +#include <mln/core/alias/neighb3d.hh> + +#include <mln/core/image/image3d.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dump/all.hh> +#include <mln/io/dicom/load.hh> + +#include <mln/value/label_16.hh> + +#include "seg_vol_irm.hh" + +using namespace mln; +using value::int_u12; +using value::label_16; + +template <typename I, typename N, typename L> +mln_ch_value(I, bool) +igr_seg(const Image<I>& input_, const mln::Neighborhood<N>& nbh_, L& nlabels); + +int usage() +{ + std::cout << "Usage: ./seg3d image.dcm" << std::endl; + + return 1; +} + +int main(int argc, char* argv[]) +{ + if (argc < 2) + return usage(); + + label_16 nlabels; + image3d<int_u12> ima; + io::dicom::load(ima, argv[1]); + io::dump::save(igr_seg(ima, c6(), nlabels), "result.dump"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/launch3d.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/launch3d.sh (revision 0) +++ trunk/milena/sandbox/fabien/igr/launch3d.sh (revision 3471) @@ -0,0 +1,24 @@ +#!/bin/sh + +process_file () +{ + ./seg3d $1 + ../bin/dump2pbm result.dump results/${2}_06_result.pbm + + ../bin/dump2pbm bg.dump results/${2}_01_bg.pbm + ../bin/dump2pbm obj.dump results/${2}_03_obj.pbm + + ../bin/dump2pbm bg_closed.dump results/${2}_02_bg_closed.pbm + ../bin/dump2pbm obj_closed.dump results/${2}_04_obj_closed.pbm + +#../bin/dump2ppm regions_color.dump results/${2}_05_colors.ppm + + rm *.dump + + mv bg_histo_norm.plot results/${2}_bg.plot + mv obj_histo_norm.plot results/${2}_obj.plot +} + +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0052.dcm" "52" +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0061.dcm" "61" +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0064.dcm" "64" Property changes on: trunk/milena/sandbox/fabien/igr/launch3d.sh ___________________________________________________________________ Name: svn:executable + * Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3470) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3471) @@ -1,15 +1,20 @@ -all: seg_vol_irm.cc - g++ -I../../../ \ - -I/Users/HiSoKa/Downloads/gdcm-2.0.10/Source/Common/ \ +DICOM_INC = -I/Users/HiSoKa/Downloads/gdcm-2.0.10/Source/Common/ \ -I/Users/HiSoKa/Downloads/gdcmbin/Source/Common/ \ -I/Users/HiSoKa/Downloads/gdcm-2.0.10/Source/DataDictionary/ \ -I/Users/HiSoKa/Downloads/gdcm-2.0.10/Source/MediaStorageAndFileFormat/ \ - -I/Users/HiSoKa/Downloads/gdcm-2.0.10/Source/DataStructureAndEncodingDefinition/ \ - -L/Users/HiSoKa/Downloads/gdcmbin/bin \ + -I/Users/HiSoKa/Downloads/gdcm-2.0.10/Source/DataStructureAndEncodingDefinition/ +DICOM_LIB = -L/Users/HiSoKa/Downloads/gdcmbin/bin \ -lgdcmCommon -lgdcmDICT -lgdcmDSED -lgdcmIOD -lgdcmMSFF -lgdcmexpat -lgdcmjpeg12 -lgdcmjpeg16 -lgdcmjpeg8 -lgdcmopenjpeg -lgdcmuuid -lgdcmzlib \ - -framework CoreFoundation \ - -DNDEBUG -O1 \ - seg_vol_irm.cc -o seg_vol_irm + -framework CoreFoundation +CXXFLAGS = -DNDEBUG -O1 + +all: 2d 3d + +2d: seg_vol_irm.hh seg2d.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} seg2d.cc -o seg2d + +3d: seg_vol_irm.hh seg3d.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} seg3d.cc -o seg3d clean: - rm -rf *.dump *.p?m + rm -rf *.dump *.p?m *.plot *.log *.csv Index: trunk/milena/sandbox/fabien/TODO =================================================================== --- trunk/milena/sandbox/fabien/TODO (revision 0) +++ trunk/milena/sandbox/fabien/TODO (revision 3471) @@ -0,0 +1,13 @@ + ___________ +< TODO list > + ----------- + \ ^__^ + \ (oo)\_______ + (__)\ )\/\ + ||----w | + || || + +[X] Create binaries +[X] Create scripts shell +[ ] Check standard deviation + Index: trunk/milena/sandbox/fabien/bin/dump2ppm.cc =================================================================== --- trunk/milena/sandbox/fabien/bin/dump2ppm.cc (revision 0) +++ trunk/milena/sandbox/fabien/bin/dump2ppm.cc (revision 3471) @@ -0,0 +1,35 @@ + +#include <mln/core/image/image2d.hh> +#include <mln/make/image3d.hh> +#include <mln/debug/slices_2d.hh> + +#include <mln/value/rgb8.hh> +#include <mln/io/dump/load.hh> +#include <mln/io/ppm/save.hh> + +#include <mln/literal/colors.hh> + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.dump output.ppm" << std::endl; + abort(); +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::rgb8; + + if (argc != 3) + usage(argv); + + image3d<rgb8> vol; + io::dump::load(vol, argv[1]); + + rgb8 bg = literal::black; + image2d<rgb8> ima = debug::slices_2d(vol, 1.f, bg); + io::ppm::save(ima, argv[2]); +} Index: trunk/milena/sandbox/fabien/bin/dump2pbm.cc =================================================================== --- trunk/milena/sandbox/fabien/bin/dump2pbm.cc (revision 0) +++ trunk/milena/sandbox/fabien/bin/dump2pbm.cc (revision 3471) @@ -0,0 +1,35 @@ + +#include <mln/core/image/image2d.hh> +#include <mln/make/image3d.hh> +#include <mln/debug/slices_2d.hh> + +#include <mln/value/int_u16.hh> +#include <mln/io/dump/load.hh> +#include <mln/io/pbm/save.hh> + +#include <mln/literal/colors.hh> + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.dump output.pbm" << std::endl; + abort(); +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u16; + + if (argc != 3) + usage(argv); + + image3d<bool> vol; + io::dump::load(vol, argv[1]); + + int_u16 bg = 0; + image2d<bool> ima = debug::slices_2d(vol, 1.f, bg); + io::pbm::save(ima, argv[2]); +}