/* -*- 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: distance3d.cc,v 1.2 2006-08-31 08:09:55 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 <fstream>
#include <string>
#include <stdexcept>
#include <dlfcn.h>
#include <libmona/2DImageWrap.hh>
#include <libmona/monaException.hh>
#include <libmona/monaHistory.hh>
#include <libmona/3DVector.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();
void initialise(const C2DImageWrap& wimage);
void read( const C2DImageWrap& wimage, size_t q);
void get_slice(size_t s, const C2DImageWrap& image, vector<pair<C3DBounds, float> >& result) 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;
};
C3DDT::C3DDT()
{
}
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);
}
}
void C3DDT::get_slice(size_t s, const C2DImageWrap& image, vector<pair<C3DBounds, float> >& result) const
{
cvdebug() << "get slice " << s << "\n";
const C2DBitImage *src = image.getC2DBitImage();
if (!src) {
stringstream err;
err << "input image " << s << "not of type bit";
throw invalid_argument(err.str());
}
if (src->get_size() != _M_size) {
stringstream err;
err << "input image " << s << "has a dffernt size then reference";
throw invalid_argument(err.str());
}
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);
C2DBitImage::const_iterator isrc = src->begin();
C2DFImage::iterator iref = slice_tmp.begin();
for (size_t y = 0; y < _M_size.y; ++y)
for (size_t x = 0; x < _M_size.x; ++x, ++iref, ++isrc)
if (*isrc)
result.push_back(pair<C3DBounds, float>(C3DBounds(x,y,s), sqrt(*iref)));
}
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};
/*
float out_3d_1[9] = { 0, 1, 2,
1, sqrt(2.0f), sqrt(5.0f),
2, sqrt(5.0f), sqrt(8.0f)};
float out_3d_2[9] = { 1, sqrt(2.0f), sqrt(5.0f),
sqrt(2.0f), sqrt(3.0f), sqrt(6.0f),
0.0f, sqrt(6.0f), 3};
*/
vector<pair<C3DBounds, float> > ref_result;
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(0,0,0), 0));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(1,0,0), 1));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(2,0,0), 2));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(0,1,0), 1));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(1,1,0), sqrt(2.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(2,1,0), sqrt(5.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(0,2,0), 2));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(1,2,0), sqrt(5.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(2,2,0), sqrt(8.0f)));
//ref_result.push_back(pair<C3DBounds, float>(C3DBounds(0,0,1), 0));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(1,0,1), sqrt(2.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(2,0,1), sqrt(5.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(0,1,1), sqrt(2.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(1,1,1), sqrt(3.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(2,1,1), sqrt(6.0f)));
//ref_result.push_back(pair<C3DBounds, float>(C3DBounds(0,2,1), 0));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(1,2,1), sqrt(6.0f)));
ref_result.push_back(pair<C3DBounds, float>(C3DBounds(2,2,1), 3));
bool ref_3d_1[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1};
bool ref_3d_2[9] = { 0, 1, 1, 1, 1, 1, 0, 1, 1};
C2DBounds Size2D(3,3);
C2DBitImage *src1 = new C2DBitImage(Size2D);
C2DBitImage *src2 = new C2DBitImage(Size2D);
C2DBitImage *ref1 = new C2DBitImage(Size2D);
C2DBitImage *ref2 = 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());
copy(&ref_3d_1[0], &ref_3d_1[9], ref1->begin());
copy(&ref_3d_2[0], &ref_3d_2[9], ref2->begin());
C2DImageWrap s1(src1);
C2DImageWrap s2(src2);
C2DImageWrap wref1(ref1);
C2DImageWrap wref2(ref2);
initialise(s1);
if (!test_initialise(Size2D)) {
cverr() << "initialisation failed!\n";
return false;
}
read(s2, 1);
read(s1, 2);
read(s2, 3);
read(s1, 4);
// should implement another test
vector<pair<C3DBounds, float> > result;
get_slice(0, wref1, result);
get_slice(1, wref2, result);
if (!equal(result.begin(), result.end(), ref_result.begin())) {
cverr() << "part 3 failed\n";
return false;
}
return true;
}
int main( int argc, const char *argv[] )
{
string src_filename;
string ref_filename;
string out_filename;
string out_type;
vector<string> filter_chain;
C2DImageIOPluginHandler imageio;
bool self_test = false;
popt::COptions options;
options.push_back(popt::option( src_filename, "in-file", 'i', "input image(s) to be filtered", NULL));
options.push_back(popt::option( ref_filename, "ref-file", 'r', "reference to evaluate the distance from", NULL));
options.push_back(popt::option( out_filename, "out-file", 'o', "output file name", NULL));
// options.push_back(popt::option( out_type, imageio.get_set(), "type", 't',"output file type" , NULL));
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;
return dt.test() ? EXIT_SUCCESS : EXIT_FAILURE;
}
cvdebug() << "IO supported types: " << imageio.plugin_names() << "\n";
if ( src_filename.empty() )
throw mona_runtime_error("'--in-file' ('i') option required");
if ( ref_filename.empty() )
throw mona_runtime_error("'--ref-file' ('r') option required");
if ( out_filename.empty() )
throw mona_runtime_error("'--out-file' ('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 ref_basename = get_filename_pattern_and_range(ref_filename, start_filenum, end_filenum, format_width);
if (start_filenum >= end_filenum)
throw invalid_argument(string("no files match pattern ") + ref_basename);
char new_line = cverb.show_debug() ? '\n' : '\r';
// time_t start_time = time(NULL);
bool first = true;
C3DDT dt;
for (size_t i = start_filenum, k=0; i < end_filenum; ++i, ++k) {
string ref_name = create_filename(ref_basename.c_str(), i);
cvmsg() << new_line << "Read: " << i <<" out of "<< "[" << start_filenum<< "," << end_filenum << "]\r" ;
auto_ptr<C2DImageList> in_image_list(imageio.load(ref_name));
if (in_image_list.get() && in_image_list->size()) {
if (first) {
dt.initialise(*in_image_list->begin());
first = false;
}else
dt.read(*in_image_list->begin(), k);
}
}
size_t src_start_filenum = 0;
size_t src_end_filenum = 0;
std::string src_basename = get_filename_pattern_and_range(src_filename, src_start_filenum, src_end_filenum, format_width);
if (src_start_filenum >= src_end_filenum)
throw invalid_argument(string("no files match pattern ") + src_basename);
if (src_start_filenum != start_filenum || src_end_filenum != end_filenum)
throw invalid_argument(string("reference range and in-range are not equal"));
vector<pair<C3DBounds, float> > result;
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()) {
dt.get_slice(k, *in_image_list->begin(), result);
}
}
cvmsg() << "\n";
ofstream os( out_filename.c_str(), ios::out );
for (vector<pair<C3DBounds, float> >::const_iterator i = result.begin();
i != result.end(); ++i)
os << i->first.x << " " << i->first.y << " " << i->first.z << " " << i->second << '\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;
}