URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-05 Fabien Freling <fabien.freling@lrde.epita.fr> Add source for dicom2dump tool. * fabien/TODO: Simple TODO reminder. * fabien/bin/Makefile: Makefile to compile dicom binaries. * fabien/bin/dicom2dump.cc: Binary to convert dicom to dump. * fabien/bin/dicom_mask.cc: Binary to extract simple mask from dicom images. * fabien/igr/Makefile: Update. * fabien/igr/launch2d.sh: Update. * fabien/igr/launch3d.sh: Update. * fabien/igr/seg2d.cc: Update. * fabien/igr/seg3d.cc: Update. * fabien/igr/seg_vol_irm.hh: Implement different threshold techniques: double and deviation. --- TODO | 13 +++++++- bin/Makefile | 23 ++++++++++++++ bin/dicom2dump.cc | 32 ++++++++++++++++++++ bin/dicom_mask.cc | 38 ++++++++++++++++++++++++ igr/Makefile | 6 +++ igr/launch2d.sh | 14 ++++++--- igr/launch3d.sh | 12 +++++-- igr/seg2d.cc | 17 ++++++++++ igr/seg3d.cc | 17 ++++++++++ igr/seg_vol_irm.hh | 82 +++++++++++++++++++++++++++++++++++++++-------------- 10 files changed, 221 insertions(+), 33 deletions(-) Index: trunk/milena/sandbox/fabien/igr/seg_vol_irm.hh =================================================================== --- trunk/milena/sandbox/fabien/igr/seg_vol_irm.hh (revision 3481) +++ trunk/milena/sandbox/fabien/igr/seg_vol_irm.hh (revision 3482) @@ -39,6 +39,7 @@ #include <mln/value/int_u12.hh> #include <mln/value/int_u16.hh> #include <mln/value/rgb8.hh> +#include <mln/value/label_16.hh> #include <mln/io/pgm/all.hh> #include <mln/io/pbm/all.hh> @@ -56,6 +57,7 @@ #include <mln/labeling/blobs.hh> #include <mln/labeling/compute.hh> +#include <mln/labeling/foreground.hh> #include <mln/labeling/fill_holes.hh> #include <mln/labeling/n_max.hh> @@ -85,7 +87,6 @@ #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> @@ -149,12 +150,30 @@ +template <typename I> +void +save_regions_color(const Image<I>& ima, metal::int_<2>) +{ + io::ppm::save(ima, "regions_color.ppm"); +} + + + +template <typename I> +void +save_regions_color(const Image<I>& ima, metal::int_<3>) +{ + io::dump::save(ima, "regions_color.dump"); +} + + + 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; + int bg_thres = 20; + int obj_thres = 10; mln_ch_value(I, bool) ima_bg; initialize(ima_bg, input); @@ -166,15 +185,16 @@ unsigned result = 0; - histo::array<mln_value(I)> arr_histo = histo::compute(input); + // We remove the 0 value because it is not part of the image. + histo::array<mln_value(I)> arr_histo = histo::compute(input | pw::value(input) != 0); 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; - + std::ofstream fout("histo.plot"); + fout << "0 0" << std::endl; for (unsigned int i = 1; i < ima_histo.nelements(); ++i) { + fout << i << " " << ima_histo(point1d(i)) << std::endl; ima_histo(point1d(i)) += ima_histo(point1d(i - 1)); } accu::max<unsigned> max_accu; @@ -207,26 +227,41 @@ 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); + ima_bg = close_threshold(ima_bg, metal::int_<I::site::dim>(), 9, 15); + ima_obj = close_threshold(ima_obj, metal::int_<I::site::dim>(), 9, 11); // 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); + save_regions_color(out, metal::int_<I::site::dim>()); 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"); } + // Labeling + /*label_16 nlabels = 0; + + mln_ch_value(I, label_16) bg_labels = labeling::foreground(ima_bg, nbh, nlabels); + mln_ch_value(I, label_16) obj_labels = labeling::foreground(ima_obj, nbh, nlabels); + + accu::count<int_u8> a_; + util::array<unsigned> arr_label = labeling::compute(a_, ima_bg, bg_labels, nlabels); + util::array<label_16> arr_big = labeling::n_max<label_16>(arr_label, 1); + data::fill((ima_bg | (pw::value(bg_labels) != pw::cst(arr_big[1]))).rw(), false); + + arr_label = labeling::compute(a_, ima_obj, obj_labels, nlabels); + arr_big = labeling::n_max<label_16>(arr_label, 1); + data::fill((ima_obj | (pw::value(obj_labels) != pw::cst(arr_big[1]))).rw(), false);*/ + + // Histo 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); @@ -236,23 +271,27 @@ 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"); + std::ofstream fout_bg("bg_histo.plot"); + std::ofstream fout_p_bg("bg_histo_norm.plot"); for (unsigned int i = 0; i < ima_bg_histo.nelements(); ++i) { + fout_bg << i << " " << ima_bg_histo(point1d(i)) << std::endl; ima_bg_histo(point1d(i)) *= 10000; ima_bg_histo(point1d(i)) /= bg_sum; - fout_bg << i << " " << ima_bg_histo(point1d(i)) << std::endl; + fout_p_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"); + std::ofstream fout_obj("obj_histo.plot"); + std::ofstream fout_p_obj("obj_histo_norm.plot"); for (unsigned int i = 0; i < ima_obj_histo.nelements(); ++i) { + fout_obj << i << " " << ima_obj_histo(point1d(i)) << std::endl; ima_obj_histo(point1d(i)) *= 10000; ima_obj_histo(point1d(i)) /= obj_sum; - fout_obj << i << " " << ima_obj_histo(point1d(i)) << std::endl; + fout_p_obj << i << " " << ima_obj_histo(point1d(i)) << std::endl; } // Search for the index with the min distance between histogrammes. @@ -275,7 +314,7 @@ unsigned find_threshold_mean(const Image<I>& input, const Neighborhood<N>& nbh) { - unsigned result = 0; + unsigned coef = 1; accu::mean<unsigned> mean_accu; unsigned mean = level::compute(mean_accu, (input | (pw::value(input) != 0))); @@ -283,13 +322,13 @@ 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); + std::cout << "[mean = " << mean << " | deviation = " << deviation << "]"; + return floor(mean + coef * deviation); } -template <typename I, typename N, typename L> +/*template <typename I, typename N, typename L> mln_ch_value(I, bool) igr_seg(const Image<I>& input_, const Neighborhood<N>& nbh_, L& nlabels) { @@ -301,9 +340,10 @@ mln_ch_value(I, bool) ima_thres; initialize(ima_thres, input); data::fill(ima_thres, false); - //unsigned threshold_value = find_threshold_value(input, nbh); + std::cout << "double threshold value = " << find_threshold_value(input, nbh) << std::endl; unsigned threshold_value = find_threshold_mean(input, nbh); + std::cout << " deviation threshold value = " << threshold_value << std::endl; 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 3481) +++ trunk/milena/sandbox/fabien/igr/seg2d.cc (revision 3482) @@ -34,7 +34,22 @@ label_16 nlabels; image2d<int_u12> ima; io::dicom::load(ima, argv[1]); - io::pbm::save(igr_seg(ima, c4(), nlabels), "result.pbm"); + + image2d<bool> ima_thres; + initialize(ima_thres, ima); + data::fill(ima_thres, false); + unsigned threshold_value = find_threshold_value(ima, c4()); + std::cout << "double threshold value = " << threshold_value << std::endl; + data::fill((ima_thres | pw::value(ima) < pw::cst(threshold_value)).rw(), true); + + io::pbm::save(ima_thres, "result_double.pbm"); + + data::fill(ima_thres, false); + threshold_value = find_threshold_mean(ima, c4()); + std::cout << " deviation threshold value = " << threshold_value << std::endl; + data::fill((ima_thres | pw::value(ima) < pw::cst(threshold_value)).rw(), true); + + io::pbm::save(ima_thres, "result_deviation.pbm"); return 0; } Index: trunk/milena/sandbox/fabien/igr/launch2d.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/launch2d.sh (revision 3481) +++ trunk/milena/sandbox/fabien/igr/launch2d.sh (revision 3482) @@ -3,8 +3,9 @@ process_file () { ./seg2d $1 - if [ -f result.pbm ]; then - mv result.pbm results/${2}_06_result.pbm + if [ -f result_double.pbm ]; then + mv result_double.pbm results/${2}_06_result_double.pbm + mv result_deviation.pbm results/${2}_07_result_deviation.pbm fi if [ -f bg.pbm ]; then @@ -14,10 +15,15 @@ mv bg_closed.pbm results/${2}_02_bg_closed.pbm mv obj_closed.pbm results/${2}_04_obj_closed.pbm + if [ -f regions_color.ppm ]; then mv regions_color.ppm results/${2}_05_regions_color.ppm + fi - mv bg_histo_norm.plot results/${2}_bg.plot - mv obj_histo_norm.plot results/${2}_obj.plot + mv histo.plot results/${2}.plot + mv bg_histo.plot results/${2}_bg.plot + mv obj_histo.plot results/${2}_obj.plot + mv bg_histo_norm.plot results/${2}_p_bg.plot + mv obj_histo_norm.plot results/${2}_p_obj.plot fi } Index: trunk/milena/sandbox/fabien/igr/seg3d.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/seg3d.cc (revision 3481) +++ trunk/milena/sandbox/fabien/igr/seg3d.cc (revision 3482) @@ -36,7 +36,22 @@ label_16 nlabels; image3d<int_u12> ima; io::dicom::load(ima, argv[1]); - io::dump::save(igr_seg(ima, c6(), nlabels), "result.dump"); + + image3d<bool> ima_thres; + initialize(ima_thres, ima); + data::fill(ima_thres, false); + unsigned threshold_value = find_threshold_value(ima, c6()); + std::cout << "double threshold value = " << threshold_value << std::endl; + data::fill((ima_thres | pw::value(ima) < pw::cst(threshold_value)).rw(), true); + + io::dump::save(ima_thres, "result_double.dump"); + + data::fill(ima_thres, false); + threshold_value = find_threshold_mean(ima, c6()); + std::cout << " deviation threshold value = " << threshold_value << std::endl; + data::fill((ima_thres | pw::value(ima) < pw::cst(threshold_value)).rw(), true); + + io::dump::save(ima_thres, "result_deviation.dump"); return 0; } Index: trunk/milena/sandbox/fabien/igr/launch3d.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/launch3d.sh (revision 3481) +++ trunk/milena/sandbox/fabien/igr/launch3d.sh (revision 3482) @@ -3,7 +3,8 @@ process_file () { ./seg3d $1 - ../bin/dump2pbm result.dump results/${2}_06_result.pbm + ../bin/dump2pbm result_double.dump results/${2}_06_result_double.pbm + ../bin/dump2pbm result_deviation.dump results/${2}_07_result_deviation.pbm ../bin/dump2pbm bg.dump results/${2}_01_bg.pbm ../bin/dump2pbm obj.dump results/${2}_03_obj.pbm @@ -11,12 +12,15 @@ ../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 + ../bin/dump2ppm regions_color.dump results/${2}_05_regions_colors.ppm rm *.dump - mv bg_histo_norm.plot results/${2}_bg.plot - mv obj_histo_norm.plot results/${2}_obj.plot + mv histo.plot results/${2}.plot + mv bg_histo.plot results/${2}_bg.plot + mv obj_histo.plot results/${2}_obj.plot + mv bg_histo_norm.plot results/${2}_p_bg.plot + mv obj_histo_norm.plot results/${2}_p_obj.plot } process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0052.dcm" "52" Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3481) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3482) @@ -16,5 +16,11 @@ 3d: seg_vol_irm.hh seg3d.cc g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} seg3d.cc -o seg3d +grad: grad_clo_and_wshd.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o grad_clo + +wst: wst_rag.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o wst_rag + clean: rm -rf *.dump *.p?m *.plot *.log *.csv Index: trunk/milena/sandbox/fabien/TODO =================================================================== --- trunk/milena/sandbox/fabien/TODO (revision 3481) +++ trunk/milena/sandbox/fabien/TODO (revision 3482) @@ -9,5 +9,14 @@ [X] Create binaries [X] Create scripts shell -[ ] Check standard deviation - +[X] Generate color images for regions +[X] Generate histograms (normal, bg, obj, p_bg, p_obj) +[ ] Test processing chain on US +[X] Create README file for special images +[ ] Implement watershed +[X] Create tool for dicom mask +[X] Check conversion for mask (everything except int_u12(0)) +[ ] Extract ROI (projection or fill holes) +[ ] Create routine for region colors +[ ] Create dicom2dump (2d, 3d, int_u8, int_u12) +[ ] After threshold, take biggest object and then erode, dilate Index: trunk/milena/sandbox/fabien/bin/dicom2dump.cc =================================================================== --- trunk/milena/sandbox/fabien/bin/dicom2dump.cc (revision 0) +++ trunk/milena/sandbox/fabien/bin/dicom2dump.cc (revision 3482) @@ -0,0 +1,32 @@ +#include <mln/core/concept/image.hh> +#include <mln/core/image/image3d.hh> + +#include <mln/value/int_u12.hh> +#include <mln/io/dicom/load.hh> +#include <mln/io/dump/save.hh> + + + +int usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.dcm output.dump" << std::endl; + std::cerr << "\t work for 3D images encoded in int_u12" << std::endl; + return 1; +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u12; + + if (argc != 3) + return usage(argv); + + image3d<int_u12> ima; + io::dicom::load(ima, argv[1]); + io::dump::save(ima, argv[2]); + + return 0; +} Index: trunk/milena/sandbox/fabien/bin/dicom_mask.cc =================================================================== --- trunk/milena/sandbox/fabien/bin/dicom_mask.cc (revision 0) +++ trunk/milena/sandbox/fabien/bin/dicom_mask.cc (revision 3482) @@ -0,0 +1,38 @@ +#include <mln/core/concept/image.hh> +#include <mln/core/image/image2d.hh> + +#include <mln/value/int_u8.hh> +#include <mln/io/dicom/load.hh> +#include <mln/io/pbm/save.hh> + +#include <mln/literal/colors.hh> + +#include <mln/level/transform.hh> +#include <mln/fun/v2b/threshold.hh> + +#include <mln/data/fill.hh> +#include <mln/pw/all.hh> + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.dcm output.pbm" << std::endl; + abort(); +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 3) + usage(argv); + + image2d<int_u8> input; + io::dicom::load(input, argv[1]); + + image2d<bool> ima = level::transform(input, fun::v2b::threshold<int_u8>(1)); + io::pbm::save(ima, argv[2]); +} Index: trunk/milena/sandbox/fabien/bin/Makefile =================================================================== --- trunk/milena/sandbox/fabien/bin/Makefile (revision 0) +++ trunk/milena/sandbox/fabien/bin/Makefile (revision 3482) @@ -0,0 +1,23 @@ +GDCM_SRC_DIR = /Users/HiSoKa/Downloads/gdcm-2.0.10 +GDCM_BIN_DIR = /Users/HiSoKa/Downloads/gdcmbin + +DICOM_INC = -I${GDCM_SRC_DIR}/Source/Common/ \ + -I${GDCM_BIN_DIR}/Source/Common/ \ + -I${GDCM_SRC_DIR}/Source/DataDictionary/ \ + -I${GDCM_SRC_DIR}/Source/MediaStorageAndFileFormat/ \ + -I${GDCM_SRC_DIR}/Source/DataStructureAndEncodingDefinition/ + +# "-framework CoreFoundation" is a Mac OS X specific flag +DICOM_LIB = -L${GDCM_BIN_DIR}/bin \ + -lgdcmCommon -lgdcmDICT -lgdcmDSED -lgdcmIOD -lgdcmMSFF -lgdcmexpat -lgdcmjpeg12 -lgdcmjpeg16 -lgdcmjpeg8 -lgdcmopenjpeg -lgdcmuuid -lgdcmzlib \ + -framework CoreFoundation + +CXXFLAGS = -DNDEBUG -O1 + +all: mask dump + +mask: dicom_mask.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o dicom_mask + +dump: dicom2dump.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o dicom2dump