/* -*- 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: defmod.cc 938 2006-07-11 11:57:01Z write1 $ */
#if HAVE_CONFIG_H
#include <config.h>
#endif
#include <string>
// MONA specific
#include <mona.hh>
#include <set>
using namespace std;
using namespace mona;
class CImage2Float: public TUnaryImageFilter<C3DFImage> {
public:
template <class Data3D>
typename CImage2Float::result_type operator () (Data3D& data) {
C3DFImage result(data.get_size());
copy(data.begin(), data.end(), result.begin());
return result;
}
};
typedef set<const CTriangleMesh::vertex_type*> CNeighbors;
struct SLocation {
CTriangleMesh::normal_type *normal;
CNeighbors vertices;
};
typedef map<CTriangleMesh::vertex_type*, SLocation> CModel;
void mesh_deform(CTriangleMesh& mesh, const C3DFImage& ref_image, const C3DFVectorfield& gradient,
float iso_value, float w1, float w2, float w3, float kappa, int max_iter)
{
// build model
CModel model;
mesh.evaluate_normals();
CTriangleMesh::vertex_iterator vb = mesh.vertices_begin();
CTriangleMesh::vertex_iterator ve = mesh.vertices_end();
CTriangleMesh::normal_iterator nb = mesh.normals_begin();
while (vb != ve) {
SLocation loc = {&*nb};
model[&*vb] = loc;
++vb;
++nb;
}
CTriangleMesh::triangle_iterator tb = mesh.triangles_begin();
CTriangleMesh::triangle_iterator te = mesh.triangles_end();
while (tb != te) {
const CTriangleMesh::vertex_type& a = mesh.vertex_at(tb->x);
const CTriangleMesh::vertex_type& b = mesh.vertex_at(tb->y);
const CTriangleMesh::vertex_type& c = mesh.vertex_at(tb->z);
SLocation& na = model[const_cast<CTriangleMesh::vertex_type*>(&a)];
na.vertices.insert(&a);
na.vertices.insert(&b);
na.vertices.insert(&c);
SLocation& nb = model[const_cast<CTriangleMesh::vertex_type*>(&b)];
nb.vertices.insert(&a);
nb.vertices.insert(&b);
nb.vertices.insert(&c);
SLocation& nc = model[const_cast<CTriangleMesh::vertex_type*>(&c)];
nc.vertices.insert(&a);
nc.vertices.insert(&b);
nc.vertices.insert(&c);
++tb;
}
CModel::iterator model_end = model.end();
double old_sum_dist = numeric_limits<double>::max();
double sum_dist = 0.0;
double max_shift = numeric_limits<double>::max();
while (max_iter-- && max_shift > 0.001) {
old_sum_dist = sum_dist;
sum_dist = 0.0;
max_shift = 0.0;
CModel::iterator model_i = model.begin();
while (model_i != model_end) {
// values for external forces
float iso_delta = (iso_value - ref_image(*model_i->first))/kappa;
float grad_scale = (gradient(*model_i->first) * *model_i->second.normal);
float f3 = tanh( iso_delta );
float f2 = grad_scale * f3 / 100.0;
float f1 = w2 * f2 + w3 * f3;
C3DFVector center(0,0,0);
if (model_i->second.vertices.size() > 0) {
// values for internal force
CNeighbors::const_iterator ni = model_i->second.vertices.begin();
CNeighbors::const_iterator ne = model_i->second.vertices.end();
while (ni != ne) {
center += **ni;
++ni;
}
center /= model_i->second.vertices.size() ;
}else
center = *model_i->first;
C3DFVector shift = w1 * (center - *model_i->first) + f1 * *model_i->second.normal;
sum_dist += iso_delta > 0 ? iso_delta : -iso_delta;
*model_i->first += shift;
float snorm = shift.norm();
if (max_shift < snorm)
max_shift = snorm;
++model_i;
}
mesh.evaluate_normals();
cvmsg() << "iterations to go " << max_iter << "; distance = " << sum_dist << "; max shift = " << max_shift << " \r";
}
cvmsg() << endl;
}
void scale_gradient (C3DFVectorfield & field)
{
C3DFVectorfield::iterator i = field.begin();
C3DFVectorfield::iterator e = field.end();
while ( i != e ) {
float norm2 = i->norm2();
if (norm2 > 0.0)
*i /= norm2;
++i;
}
}
int main(int argc, const char *argv[])
{
try {
string in_filename;
string out_filename;
string ref_filename;
float w1, w2, w3;
float kappa;
float iso_value;
int maxiter;
int gauss;
popt::COptions opts;
opts.push_back(popt::option( in_filename, "in-mesh", 'i', "input mesh", NULL ));
opts.push_back(popt::option( out_filename, "out-mesh", 'o', "defirmed mesh", "-" ));
opts.push_back(popt::option( ref_filename, "ref-image", 'r', "reference image", NULL));
opts.push_back(popt::option( iso_value, "iso-value", 's', "iso-value for mesh adaption", "64"));
opts.push_back(popt::option( w1, "w1", 'a', "weight for mesh smoothing force", "0.04"));
opts.push_back(popt::option( w2, "w2", 'b', "weight for gradient force", "0.04"));
opts.push_back(popt::option( w3, "w3", 'c', "weight for intensity adaptinf force", "0.02"));
opts.push_back(popt::option( kappa, "kappa", 'k', "scaling factor", "1.0"));
opts.push_back(popt::option( gauss, "gauss", 'g', "gauss filter put on the image", "2"));
opts.push_back(popt::option( maxiter,"maxiter", 'm', "maximum number of iterations", "200"));
vector<string> non_options;
popt::parse_options(argc, argv, opts, non_options);
if (in_filename.empty()) {
throw mona_runtime_error("'--in-mesh' ('i') option required");
}
if (ref_filename.empty()) {
throw mona_runtime_error("'--ref-image' ('r') option required");
}
if (non_options.size() > 0) {
cverr() << "ignoring unknown options" << endl;
}
CIOMeshPluginHandler meshio;
auto_ptr<CTriangleMesh> mesh_ptr(meshio.load(in_filename));
if ( !mesh_ptr->get_available_data() ) {
string not_found = ("No supported data found in ") + in_filename;
throw mona_runtime_error(not_found);
}
C3DImageIOPluginHandler imageio;
auto_ptr<C3DImageList> images(imageio.load(ref_filename));
if (images.get() && images->size() > 0) {
CImage2Float img2float;
C3DFImage ref_image = image_filter(img2float, *images->begin());
if (gauss > 0) {
cvmsg() << "gauss filtering image" << endl;
gauss_3d(&ref_image, gauss, gauss, gauss);
}
cvmsg() << "evaluate gradient" << endl;
C3DFVectorfield grad_field= ref_image.get_gradient_field<C3DFVectorfield>();
scale_gradient( grad_field);
cvmsg() << "start deform ..." << endl;
mesh_deform(*mesh_ptr, ref_image, grad_field, iso_value, w1, w2, w3, kappa, maxiter);
cvmsg() << "finished deform." << endl;
}else {
throw mona_runtime_error(string("No image found in ")+ref_filename );
}
if ( !meshio.save(mesh_ptr->get_sourceformat() , *mesh_ptr, out_filename) ) {
string not_found = ("Unable to wite mesh to ") + out_filename;
throw mona_runtime_error(not_found);
}
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 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;
}
/*
$Log$
Revision 1.2 2005/06/29 13:32:02 wollny
move to libmona-version 0.7
Revision 1.2 2005/06/02 13:33:23 gerddie
adapt code to new plugin handling
Revision 1.1.1.1 2005/03/17 13:48:32 gerddie
initial checkin
Revision 1.1 2005/01/20 17:39:37 wollny
first checkin
*/