/* -*- 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: costFkt.cc 670 2005-06-29 13:22:23Z wollny $
/*! \brief implementation of costFkt- classes
\author G. Wollny, tittge@cbs.mpg.de, 2004
*/
#pragma implementation "costFkt.hh"
#include <iostream>
#include <libmona/costFkt.hh>
namespace mona {
CCostFunction::CCostFunction()
{
}
CCostFunction::~CCostFunction()
{
}
double CCostFunction::get_cost() const
{
if (!value.is_valid())
value.set(evaluate_cost());
return value.get();
}
TCostEvaluator::~TCostEvaluator()
{
}
CSSDCostEvaluator::CSSDCostEvaluator():_M_cost(0)
{
}
void CSSDCostEvaluator::clear()
{
_M_cost = 0;
}
void CSSDCostEvaluator::accumulate(double src_val, double ref_val)
{
double h = src_val - ref_val;
_M_cost += h*h;
}
double CSSDCostEvaluator::get_value() const {
return _M_cost;
}
#ifndef NDEBUG
bool CSSDCostEvaluator::test()
{
CSSDCostEvaluator tester;
tester.clear();
for (int i = 0; i < 10; i++)
tester.accumulate(i,i);
if (tester.get_value() != 0.0)
return false;
tester.clear();
for (int i = 0; i < 3; i++)
for (int j = 0; j < 5; j++)
tester.accumulate(i,j);
if (tester.get_value() != 55.0)
return false;
return true;
}
#endif
CCCCostEvaluator::CCCCostEvaluator():
_M_sumS(0),
_M_sumR(0),
_M_sumS2(0),
_M_sumR2(0),
_M_sumDP(0),
_M_n(0)
{
}
void CCCCostEvaluator::clear()
{
_M_sumS = 0.0;
_M_sumR = 0.0;
_M_sumS2 = 0.0;
_M_sumR2 = 0.0;
_M_sumDP = 0.0;
_M_n = 0;
}
void CCCCostEvaluator::accumulate(double src_val, double ref_val)
{
_M_sumS += src_val;
_M_sumR += ref_val;
_M_sumS2 += src_val * src_val;
_M_sumR2 += ref_val * ref_val;
_M_sumDP += src_val * ref_val;
++_M_n;
}
double CCCCostEvaluator::get_value() const
{
double p1 = (_M_n * _M_sumS2 - _M_sumS * _M_sumS) *
(_M_n * _M_sumR2 - _M_sumR * _M_sumR);
// if p1 == 0, then one image has constant intensity and the registration
// can not get better
if (p1 == 0.0)
return 0.0;
double p2 = _M_n * _M_sumDP - _M_sumR * _M_sumS;
return 1.0 - p2 * p2 / p1;
}
#ifndef NDEBUG
bool CCCCostEvaluator::test()
{
CCCCostEvaluator tester;
tester.clear();
for (int i = 0; i < 10; i++)
tester.accumulate(i,i);
if ( tester.get_value() != 0.0)
return false;
tester.clear();
for (int i = 0; i < 10; i++)
for (int j = 0; j < 10; j++)
tester.accumulate(i,j);
if (tester.get_value() != 1.0)
return false;
CCCCostEvaluator tester2;
tester.clear();
for (int i = 0; i < 10; i++)
for (int j = 0; j < 5; j++){
tester.accumulate(i,j+10);
tester2.accumulate(j+10,i);
}
if (tester.get_value() != tester2.get_value())
return false;
if (tester.get_value() != tester.get_value())
return false;
return true;
}
#endif
CNMICostEvaluator::CNMICostEvaluator(double maxi, int hsize):
_M_n(0),
_M_hsize(hsize),
_M_factor(hsize/maxi),
_M_shisto(hsize),
_M_rhisto(hsize),
_M_xhisto(hsize * hsize)
{
}
CNMICostEvaluator::CNMICostEvaluator(const CNMICostEvaluator& org):
_M_n(0),
_M_hsize(org._M_hsize),
_M_factor(org._M_factor),
_M_shisto(org._M_hsize),
_M_rhisto(org._M_hsize),
_M_xhisto(org._M_hsize * org._M_hsize)
{
}
void CNMICostEvaluator::clear()
{
_M_n = 0;
fill(_M_rhisto.begin(),_M_rhisto.end(),0);
fill(_M_shisto.begin(),_M_shisto.end(),0);
fill(_M_xhisto.begin(),_M_xhisto.end(),0);
}
void CNMICostEvaluator::accumulate(double src_val, double ref_val)
{
unsigned int xc = (unsigned int) (src_val * _M_factor);
unsigned int yc = (unsigned int) (ref_val * _M_factor);
assert(xc < _M_hsize && yc < _M_hsize);
++_M_shisto[xc];
++_M_rhisto[yc];
++_M_xhisto[xc + _M_hsize * yc];
++_M_n;
}
double CNMICostEvaluator::get_value() const
{
double hrs = 0.0;
double hr = 0.0;
double hs = 0.0;
CHistogram::const_iterator ih = _M_xhisto.begin();
CHistogram::const_iterator eh = _M_xhisto.end();
while (ih != eh) {
double p = double(*ih) / _M_n;
if (p > 0.0)
hrs -= p * log(p);
++ih;
}
ih = _M_shisto.begin();
eh = _M_shisto.end();
while (ih != eh) {
double p = double(*ih) / _M_n;
if (p > 0.0)
hs -= p * log(p);
++ih;
}
ih = _M_rhisto.begin();
eh = _M_rhisto.end();
while (ih != eh) {
double p = double(*ih) / _M_n;
if (p > 0.0)
hr -= p * log(p);
++ih;
}
// both are imges have only one intensity value
// therefore, they match perfectly
if (hr + hs == 0.0)
return 0.0;
return 2.0 * hrs / (hr + hs) - 1.0;
}
#ifndef NDEBUG
bool CNMICostEvaluator::test()
{
CNMICostEvaluator tester(1024, 64);
// check ideal match
for (int i = 0; i < 1024; ++i)
tester.accumulate(i,i);
if (tester.get_value() != 0.0)
return false;
// check worst match
tester.clear();
for (int i = 0; i < 1024; ++i)
for (int j = 0; j < 1024; ++j)
tester.accumulate(i,j);
//
if ( tester.get_value() - 1.0 > 1e-10)
return false;
CNMICostEvaluator tester2(4096, 64);
// check symmetry
for (int i = 0; i < 1024; ++i)
for (int j = 0; j < 512; ++j) {
tester.accumulate(i,j+5);
tester2.accumulate(j+5,i);
}
if (tester2.get_value() - tester.get_value() > 1e-10)
return false;
// Checjk whether the values are always the same
long double r1 = tester.get_value();
long double r2 = tester.get_value();
if ( r1 != r2)
return false;
return true;
}
} // namespace mona
#endif
/* CVS LOG
$Log$
Revision 1.3 2005/06/29 13:22:22 wollny
switch to version 0.7
Revision 1.2 2005/05/31 19:21:12 gerddie
dont use namespace std, write it long
Revision 1.1.1.1 2005/03/17 13:44:20 gerddie
initial import
Revision 1.2 2004/06/03 09:57:32 wollny
Changed (hopefully) all instancable class names to Cxxxxx
Revision 1.1 2004/02/12 14:00:25 tittge
Major addaptions from libmia0
*/