/* -*- 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
*
*/
// $Id: stackdt.cc,v 1.5 2006-08-30 14:43:04 wollny Exp $
/*! \brief eva-2dimagefilter
\sa 3va-2dimagefilter.cc
\file mask.cc
\author G. Wollny, wollny eva.mpg.de, 2005
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <iostream>
#include <iomanip>
#include <string>
#include <stdexcept>
#include <dlfcn.h>
#include <libmona/2DImageWrap.hh>
#include <libmona/monaException.hh>
#include <libmona/monaHistory.hh>
#include <libmona/filetools.hh>
#include <libmona/fnameana.hh>
#include <boost/cast.hpp>
#include <boost/lambda/lambda.hpp>
using namespace std;
using namespace mona;
/* Revision string */
const char revision[] = "not specified";
using namespace boost::lambda;
const float g_far = 1e+20f;
class C3DDT {
public:
C3DDT(float scale, bool float_output);
void initialise(const C2DImageWrap& wimage);
void read( const C2DImageWrap& wimage, size_t q);
C2DImageWrap get_slice(size_t s) const;
bool test();
bool test_initialise(const C2DBounds& size);
private:
float d(float fq, float q, float fv, float v)const;
void dt1d(vector<float>& f)const;
void dt2d(C2DFImage& image)const;
struct SParabola {
int k;
int v;
float z;
float fv;
};
C2DBounds _M_size;
vector<size_t> _M_k;
vector< vector<SParabola> > _M_zdt;
float _M_scale;
bool _M_float_output;
};
C3DDT::C3DDT(float scale, bool float_output):_M_scale(scale), _M_float_output(float_output)
{
}
float C3DDT::d(float fq, float q, float fv, float v) const
{
const float d = ( fq - fv + q * q - v * v) / (q - v) * 0.5;
return d;
}
void C3DDT::dt1d(vector<float>& r)const
{
vector<float> f(r);
vector<int> v(f.size());
vector<float> z(f.size() + 1);
int k = 0;
v[0] = 0;
z[0] = -numeric_limits<float>::max();
z[1] = +numeric_limits<float>::max();
for (size_t q = 1; q < f.size(); q++) {
float s = d(f[q], q, f[v[k]], v[k]);
while (s <= z[k]) {
--k;
s = d(f[q], q, f[v[k]], v[k]);
}
++k;
v[k] = q;
z[k] = s;
z[k+1] = numeric_limits<float>::max();
}
k = 0;
for (size_t q = 0; q < f.size(); ++q) {
while (z[k+1] < q)
++k;
float delta = float(q) - v[k];
r[q] = delta * delta + f[v[k]];
}
}
void C3DDT::initialise(const C2DImageWrap& wimage)
{
const C2DBitImage *src = wimage.getC2DBitImage();
_M_size = src->get_size();
_M_k.resize(_M_size.x * _M_size.y);
fill(_M_k.begin(), _M_k.end(), 0);
_M_zdt.resize(_M_size.x * _M_size.y);
C2DBitImage::const_iterator pixel = src->begin();
vector< vector<SParabola> >::iterator k = _M_zdt.begin();
SParabola p = {0, 0, -numeric_limits<float>::max(), 0.0f};
for (size_t i = 0; i < _M_zdt.size(); ++i, ++pixel, ++k) {
p.fv = *pixel ? 0.0f : g_far;
k->push_back(p);
}
}
void C3DDT::read( const C2DImageWrap& wimage, size_t q)
{
cvdebug() << "push image " << q << "\n\n";
const C2DBitImage *src = wimage.getC2DBitImage();
if (!src)
throw runtime_error("input image is not a bit image!");
if (src->size() != _M_zdt.size())
throw runtime_error("input image has different size!");
C2DBitImage::const_iterator si = src->begin();
C2DBitImage::const_iterator ei = src->end();
vector<size_t>::iterator k = _M_k.begin();
vector< vector<SParabola> >::iterator p = _M_zdt.begin();
while (si != ei) {
float f = *si ? 0.0f : numeric_limits<float>::max();
SParabola& parabola = (*p)[*k];
float s = d (f, q, parabola.fv, parabola.v);
while (s <= parabola.z) {
--(*k);
parabola = (*p)[*k];
s = d (f, q, parabola.fv, parabola.v);
}
++(*k);
if (*k > p->size()) {
cverr() << "k = " << *k << " but p->size() = " << p->size() <<"\n";
assert(!"can't do");
}
SParabola new_p = {*k, q, s, f};
if ( *k == p->size() ) {
p->push_back(new_p);
}else {
//new_p.fv = (*p)[*k].fv;
(*p)[*k] = new_p;
if (*k < p->size() - 1)
p->resize(*k + 1);
}
++si;
++k;
++p;
}
}
void C3DDT::dt2d(C2DFImage& image)const
{
vector<float> buffer(image.get_size().x);
for (size_t y = 0; y < image.get_size().y; ++y) {
image.get_data_line_x(y, buffer);
dt1d(buffer);
image.put_data_line_x(y, buffer);
}
buffer.resize(image.get_size().y);
for (size_t x = 0; x < image.get_size().x; ++x) {
image.get_data_line_y(x, buffer);
dt1d(buffer);
image.put_data_line_y(x, buffer);
}
}
struct Convert {
Convert(float scale):_M_scale(scale){}
unsigned short operator ()(float x) {
float v = _M_scale * sqrt(x);
if (v < numeric_limits<unsigned short>::max())
return boost::numeric_cast<unsigned short>(v);
else
return numeric_limits<unsigned short>::max();
}
private:
float _M_scale;
};
C2DImageWrap C3DDT::get_slice(size_t s) const
{
cvdebug() << "get slice " << s << "\n";
C2DFImage slice_tmp(_M_size);
C2DFImage::iterator i = slice_tmp.begin();
C2DFImage::iterator e = slice_tmp.end();
vector< vector<SParabola> >::const_iterator p = _M_zdt.begin();
while (i != e) {
size_t k = 0;
while ( k < p->size() - 1 && (*p)[ k + 1 ].z < s) {
++k;
}
float delta = float(s) - (*p)[k].v;
*i = delta * delta + (*p)[k].fv;
++i;
++p;
}
dt2d(slice_tmp);
if (_M_float_output) {
C2DFImage *result = new C2DFImage(_M_size);
copy(slice_tmp.begin(), slice_tmp.end(), result->begin());
return C2DImageWrap(result);
}else{
C2DUSImage *result = new C2DUSImage(_M_size);
transform(slice_tmp.begin(), slice_tmp.end(), result->begin(), Convert(_M_scale));
return C2DImageWrap(result);
}
}
bool C3DDT::test_initialise(const C2DBounds& size)
{
if (_M_zdt.size() != size.x * size.y) {
cverr() << " buffer size should be " << size.x * size.y << " but is " << _M_zdt.size() <<"\n";
return false;
}
vector< vector<SParabola> >::const_iterator ki = _M_zdt.begin();
vector< vector<SParabola> >::const_iterator ke = _M_zdt.end();
bool result = true;
while (ki != ke) {
if (ki->size() != 1) {
cverr() << "not all slots are properly initialised in _M_zdt\n";
return false;
}
const SParabola& p = (*ki)[0];
if (p.k != 0) {
cverr() << "got a wrong k-initializer\n";
result = false;
}
if (p.v != 0) {
cverr() << "got a wrong v-initializer\n";
result = false;
}
if (p.z != -numeric_limits<float>::max()) {
cverr() << "got a wrong z-initializer\n";
result = false;
}
if ( ki == _M_zdt.begin() ) {
if (p.fv != 0.0 ) {
cverr() << "got a wrong (0,0,0) fv-initializer\n";
result = false;
}
} else {
if (p.fv != g_far) {
cverr() << "got a wrong general fv-initializer: " << p.fv <<" != " << g_far<<"\n";
result = false;
}
}
++ki;
}
return result;
}
static void print_diff(const C2DUSImage& image, unsigned short *ref)
{
size_t k = 0;
for (C2DUSImage::const_iterator i = image.begin();
i != image.end(); ++i, ++k)
cverr() << *i << " vs. ref:" << ref[k] << "\n";
}
bool C3DDT::test()
{
float in_1d[16] = { 1, 1, 1, 0,
1, 1, 1, 1,
1, 1, 0, 1,
1, 1, 1, 1 };
float out_1d[16] = { 9, 4, 1, 0,
1, 4, 9, 9,
4, 1, 0, 1,
4, 9, 16, 25 };
vector<float> src(16);
transform(&in_1d[0], &in_1d[16],src.begin(), numeric_limits<float>::max() * _1 );
dt1d(src);
if (!equal(src.begin(), src.end(), &out_1d[0])) {
cverr() << "1D distance transform failed\n";
for (size_t i = 0; i < 16; ++i) {
cverr() << out_1d[i] << " vs " << src[i] << "\n";
}
return false;
}
float in_2d[16] = { 1, 1, 1, 0,
1, 1, 1, 1,
1, 1, 0, 1,
1, 1, 1, 1 };
float out_2d[16] = { 8, 4, 1, 0,
5, 2, 1, 1,
4, 1, 0, 1,
5, 2, 1, 2 };
C2DFImage src_img(C2DBounds(4,4));
transform(&in_2d[0], &in_2d[16],src_img.begin(), numeric_limits<float>::max() * _1 );
dt2d(src_img);
if (!equal(src_img.begin(), src_img.end(), &out_2d[0])) {
cverr() << "2D distance transform failed\n";
return false;
}
bool in_3d_1[9] = { 1, 0, 0, 0, 0, 0, 0, 0, 0};
bool in_3d_2[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0};
typedef unsigned short ushort;
unsigned short out_3d_1[9] = { 0, 16, 32,
16, ushort(16* sqrt(2.0f)), ushort(16* sqrt(5.0f)),
32, ushort(16* sqrt(5.0f)), ushort(16* sqrt(8.0f))};
unsigned short out_3d_2[9] = { 16, ushort(16* sqrt(2.0f)), ushort(16* sqrt(5.0f)),
ushort(16* sqrt(2.0f)),ushort(16* sqrt(3.0f)), ushort(16* sqrt(6.0f)),
ushort(16* sqrt(5.0f)),ushort(16* sqrt(6.0f)), 48};
C2DBounds Size2D(3,3);
C2DBitImage *src1 = new C2DBitImage(Size2D);
C2DBitImage *src2 = new C2DBitImage(Size2D);
copy(&in_3d_1[0], &in_3d_1[9], src1->begin());
copy(&in_3d_2[0], &in_3d_2[9], src2->begin());
C2DImageWrap s1(src1);
C2DImageWrap s2(src2);
initialise(s1);
if (!test_initialise(Size2D)) {
cverr() << "initialisation failed!\n";
return false;
}
read(s2, 1);
read(s1, 2);
read(s2, 3);
read(s1, 4);
C2DImageWrap r1 = get_slice(0);
C2DImageWrap r2 = get_slice(1);
C2DImageWrap r3 = get_slice(2);
C2DImageWrap r4 = get_slice(3);
C2DImageWrap r5 = get_slice(4);
if (r1.get_type() != it_ushort || r2.get_type() != it_ushort || r3.get_type() != it_ushort) {
cverr() << "3D image type wrong\n";
return false;
}
if (r1.get_size() != Size2D ||r2.get_size() != Size2D ||r3.get_size() != Size2D) {
cverr() << "3D image size wrong\n";
return false;
}
C2DUSImage *rimage = r1.getC2DUSImage();
if (!equal(rimage->begin(), rimage->end(), &out_3d_1[0])) {
cverr() << "slice 0 error\n";
print_diff(*rimage, out_3d_1);
return false;
}
rimage = r2.getC2DUSImage();
if (!equal(rimage->begin(), rimage->end(), &out_3d_2[0])) {
cverr() << "slice 1 error\n";
print_diff(*rimage, out_3d_2);
return false;
}
rimage = r3.getC2DUSImage();
if (!equal(rimage->begin(), rimage->end(), &out_3d_1[0])) {
cverr() << "slice 2 error\n";
print_diff(*rimage, out_3d_1);
return false;
}
rimage = r4.getC2DUSImage();
if (!equal(rimage->begin(), rimage->end(), &out_3d_2[0])) {
cverr() << "slice 1 error\n";
print_diff(*rimage, out_3d_2);
return false;
}
rimage = r5.getC2DUSImage();
if (!equal(rimage->begin(), rimage->end(), &out_3d_1[0])) {
cverr() << "slice 2 error\n";
print_diff(*rimage, out_3d_1);
return false;
}
return true;
}
int main( int argc, const char *argv[] )
{
string in_filename;
string out_filename;
string out_type;
vector<string> filter_chain;
C2DImageIOPluginHandler imageio;
bool self_test = false;
float scale;
#ifdef HAVE_OPENEXR
const char default_out_type[]="exr";
const char default_scale[]="1.0";
#else
const char default_out_type[]="png";
const char default_scale[]="256.0";
#endif
popt::COptions options;
options.push_back(popt::option( in_filename, "in-file", 'i', "input image(s) to be filtered", NULL));
options.push_back(popt::option( out_filename, "out-base", 'o', "output file name base", "filtered"));
options.push_back(popt::option( out_type, imageio.get_set(), "type", 't',"output file type" , default_out_type));
options.push_back(popt::option( scale, "scale", 's', "scaling factor for distance image", default_scale));
options.push_back(popt::option( self_test, "self-test", 0,"run a self test" ));
try {
popt::parse_options(argc, argv, options, filter_chain);
if (self_test) {
C3DDT dt(16.0, false);
return dt.test() ? EXIT_SUCCESS : EXIT_FAILURE;
}
cvdebug() << "IO supported types: " << imageio.plugin_names() << "\n";
if ( in_filename.empty() )
throw mona_runtime_error("'--in-file' ('i') option required");
if ( out_filename.empty() )
throw mona_runtime_error("'--out-base' ('o') option required");
bool use_src_format = out_type.empty();
if (!use_src_format && !imageio.find_plugin_by_name(out_type)) {
cvwarn() << "Output file format not supported, revert to input file format\n";
use_src_format = true;
}
CHistory::instance().append(argv[0], revision, options);
size_t start_filenum = 0;
size_t end_filenum = 0;
size_t format_width = 0;
std::string src_basename = get_filename_pattern_and_range(in_filename, start_filenum, end_filenum, format_width);
if (start_filenum >= end_filenum)
throw invalid_argument(string("no files match pattern ") + src_basename);
char new_line = cverb.show_debug() ? '\n' : '\r';
bool first = true;
auto_ptr<C3DDT> dt;
for (size_t i = start_filenum, k=0; i < end_filenum; ++i, ++k) {
string src_name = create_filename(src_basename.c_str(), i);
cvmsg() << new_line << "Read: " << i <<" out of "<< "[" << start_filenum<< "," << end_filenum << "]\r" ;
auto_ptr<C2DImageList> in_image_list(imageio.load(src_name));
if (in_image_list.get() && in_image_list->size()) {
if (use_src_format) {
out_type = in_image_list->get_sourceformat();
/* here should be a test, weather the output image format supports float valued images */
if (out_type != "exr" ) {
if (scale <= 1.0)
cvwarn() << "Output image format does not support float and the distance values will not be scaled\n";
}else
scale = 1.0;
}
if (first) {
dt.reset(new C3DDT(scale, out_type == "exr"));
dt->initialise(*in_image_list->begin());
first = false;
}else
dt->read(*in_image_list->begin(), k);
}
}
for (size_t i = start_filenum, k=0; i < end_filenum; ++i, ++k) {
cvmsg() << new_line << "Create: " << i <<" out of "<< "[" << start_filenum<< "," << end_filenum << "]\r" ;
stringstream ss;
ss << out_filename << setw(format_width) << setfill('0') << i << "." << out_type;
C2DImageList list;
list.push_back(dt->get_slice(k));
if ( !imageio.save(out_type, list, ss.str()) ){
string not_save = ("unable to save result to ") + ss.str();
throw mona_runtime_error(not_save);
}
}
cvmsg() << "\n";
return EXIT_SUCCESS;
}
catch (const mona_runtime_error& e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (const mona_fatal_error& e){
cerr << argv[0] << " fatal: " << e.what() << endl;
}
catch (const mona_exception& e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (const runtime_error &e){
cerr << argv[0] << " runtime: " << e.what() << endl;
}
catch (const invalid_argument &e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (const exception& e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (...){
cerr << argv[0] << " unknown exception" << endl;
}
return EXIT_FAILURE;
}