/* -*- mona-c++ -*-
* Copyright (c) Leipzig, Madrid 2004 - 2008
* Max-Planck-Institute for Human Cognitive and Brain Science
* Max-Planck-Institute for Evolutionary Anthropology
* BIT, ETSI Telecomunicacion, UPM
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program 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 program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
/**
This program implements some automatic region growing
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <mona.hh>
#include <memory>
using namespace mona;
static const char *revision = "$Revision: 0$";
enum ENeighbourhood { nn_6, nn_18, nn_26 };
popt::option::Dictionary NeighbourOptions[] = {
{"n6", nn_6},
{"n18", nn_18},
{"n26", nn_26},
{NULL, 0}
};
class C3DRegionGrow : public TUnaryImageFilter<C3DImageWrap> {
float _M_thresh;
float _M_maxgradient;
int _M_neighborhood;
public:
C3DRegionGrow(float thresh, float maxgradient, int neighborhood);
template <class Data3D>
C3DRegionGrow::result_type operator()(const Data3D& data);
private:
template <class Data3D>
int region_grow(const C3DBounds& start_loc, int class_nr, C3DULImage& result, const Data3D& image,
typename Data3D::value_type value);
};
C3DRegionGrow::C3DRegionGrow(float thresh, float maxgradient,int neighborhood):
_M_thresh(thresh),
_M_maxgradient(maxgradient),
_M_neighborhood(neighborhood)
{
}
template <typename T>
struct location {
C3DBounds loc;
T value;
};
template <typename T>
bool operator < (const location<T>& a, const location<T>& b)
{
return (a.loc.x < b.loc.x ||
(a.loc.x == b.loc.x) && ((a.loc.y < b.loc.y) ||
(a.loc.y == b.loc.y) && (a.loc.z < b.loc.z)));
}
template <class Data3D>
int C3DRegionGrow::region_grow(const C3DBounds& start_loc, int class_nr, C3DULImage& result, const Data3D& image,
typename Data3D::value_type value)
{
typedef std::set<location<typename Data3D::value_type> > loc_stack;
loc_stack free_ends;
int npixels = 0;
location<typename Data3D::value_type> l;
l.loc = start_loc;
l.value = value;
free_ends.insert(l);
while (!free_ends.empty()) {
typename loc_stack::iterator pos = free_ends.begin();
result(pos->loc) = class_nr;
++npixels;
#define CHECK(A,B,C) l.loc = C3DBounds(pos->loc.x + A, pos->loc.y + B, pos->loc.z +C); \
if (l.loc.x >=0 &&l.loc.y >=0 && l.loc.z >= 0 && \
l.loc.x < image.get_size().x && l.loc.y < image.get_size().y && l.loc.z < image.get_size().z) { \
l.value = image(l.loc); \
if (l.value > _M_thresh) \
if (result(l.loc) == 0 && fabs(l.value - pos->value) < _M_maxgradient) \
free_ends.insert(l); }
switch (_M_neighborhood) {
case nn_26: // with corners
CHECK( 1, 1, 1);
CHECK( 1, 1,-1);
CHECK( 1,-1, 1);
CHECK( 1,-1,-1);
CHECK(-1, 1, 1);
CHECK(-1, 1,-1);
CHECK(-1,-1, 1);
CHECK(-1,-1,-1);
case nn_18: // with edges
CHECK( 1, 1, 0);
CHECK( 1,-1, 0);
CHECK(-1,-1, 0);
CHECK(-1, 1, 0);
CHECK( 0, 1,-1);
CHECK( 0, 1, 1);
CHECK( 0,-1, 1);
CHECK( 0,-1, 1);
CHECK( 1, 0, 1);
CHECK( 1, 0,-1);
CHECK(-1, 0, 1);
CHECK(-1, 0,-1);
default: // only centers
CHECK( 1, 0, 0);
CHECK( 0, 1, 0);
CHECK( 0, 0, 1);
CHECK(-1, 0, 0);
CHECK( 0,-1, 0);
CHECK( 0, 0,-1);
}
#undef CHECK
free_ends.erase(pos);
}
return npixels > 10 ? 1 : 0;
}
template <class Data3D>
C3DRegionGrow::result_type C3DRegionGrow::operator()(const Data3D& data)
{
int class_nr = 1;
C3DULImage *result = new C3DULImage(data.get_size());
typename Data3D::const_iterator inp = data.begin();
C3DULImage::iterator start_pixel = result->begin();
C3DBounds start_loc;
for (start_loc.z = 0; start_loc.z < data.get_size().z; ++start_loc.z)
for (start_loc.y = 0; start_loc.y < data.get_size().y; ++start_loc.y)
for (start_loc.x = 0; start_loc.x < data.get_size().x; ++start_loc.x, ++start_pixel, ++inp)
if (*start_pixel == 0 && *inp > _M_thresh) {
region_grow(start_loc, class_nr, *result, data, *inp);
++class_nr;
cvmsg() << "class nr = " << class_nr << std::endl;
}
return C3DImageWrap(result);
}
int main(int argc, const char *argv[])
{
std::string in_filename;
std::string out_filename;
float threshold;
float maxgradient;
int neighbourhood;
try {
popt::COptions opts;
opts.push_back(popt::option( in_filename, "in-file", 'i', "image to be segmented", NULL ));
opts.push_back(popt::option( out_filename, "b0-file", 'o', "output label image", NULL ));
opts.push_back(popt::option( threshold, "thresh", 't', "black threshold", "64"));
opts.push_back(popt::option( maxgradient, "maxgradient", 'g', "maximum gradient value", "100"));
opts.push_back(popt::option( neighbourhood, NeighbourOptions, "neighbourhood", 'n',
"neighbourhood", "n18"));
std::vector<std::string> non_options;
popt::parse_options(argc, argv, opts, non_options);
// what to do with unrecognized options
if ( non_options.size() > 0 )
throw std::invalid_argument("unknown options");
// required options (anything that has no default value)
if ( in_filename.empty() )
throw mona_runtime_error("'--in-file' ('i') option required\n");
if ( in_filename.empty() )
throw mona_runtime_error("'--cls-file' ('c') option required\n");
#ifdef DEBUG
cverb.set_verbosity( vstream::ml_debug );
#endif
// initialize the functor
C3DRegionGrow grower(threshold, maxgradient, neighbourhood);
C3DImageIOPluginHandler imageio;
// read data
std::auto_ptr<C3DImageList> inImage_list(imageio.load(in_filename));
if (!inImage_list.get() || !inImage_list->size() ) {
std::string not_found = std::string("No supported data found in ") +
in_filename;
throw mona_runtime_error(not_found);
}
for (C3DImageList::iterator image = inImage_list->begin();
image != inImage_list->end(); ++image)
*image = wrap_filter(grower, *image);
CHistory::instance().append(argv[0], revision, opts);
if ( !imageio.save(inImage_list->get_sourceformat(), *inImage_list, out_filename) ){
std::string not_save = ("unable to save result to ") + out_filename;
throw mona_runtime_error(not_save);
}
return EXIT_SUCCESS;
}
catch (const mona_runtime_error& e){
cverr() << argv[0] << " error: " << e.what() << std::endl;
}
catch (const mona_fatal_error& e){
cverr() << argv[0] << " fatal: " << e.what() << std::endl;
}
catch (const mona_exception& e){
cverr() << argv[0] << " error: " << e.what() << std::endl;
}
catch (const std::invalid_argument &e){
cverr() << argv[0] << " error: " << e.what() << std::endl;
}
catch (const std::exception& e){
cverr() << argv[0] << " error: " << e.what() << std::endl;
}
catch (...){
cverr() << argv[0] << " unknown exception" << std::endl;
}
return EXIT_FAILURE;
}