minimisation_d.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */


#ifndef AIMS_REGISTRATION_MINIMISATION_D_H
#define AIMS_REGISTRATION_MINIMISATION_D_H


#include <cartobase/type/string_conversion.h>
#include <cartobase/config/verbose.h>
#include <iostream>

#include <aims/io/motionW.h>
#include <aims/transform/least_square_estimation.h>
#include <aims/registration/minimisation.h>


//methode principale:

template <class T>
Motion Minimisation::quaternion(DisplacementField<T>& displacementField)
{
      
  //on charge l'ensemble des points a apparier
  std::vector<Point3df> ptstest = displacementField.getpointstest();
  std::vector<Point3df> ptsref = displacementField.getpointsref(); 
  int size = ptstest.size();
  unsigned sizekept = unsigned ( double(ptstest.size())*_Pkept );
  
  Motion m;

  if ((sizekept > 0) && (ptsref.size() > 0))
  {
    // Calcul de la transformation LS choisie
    switch(_transfo)
    {
      case 1:
      {
        // std::cout<<"Calcul transfo rigide..."<<std::endl;
        aims::RigidLeastSquareEstimation leastSquareEstimation(ptstest, ptsref);
        m = *(leastSquareEstimation.motion());
        break;
      }

      case 2:
      {
        // std::cout<<"Calcul transfo similitude..."<<std::endl;
        aims::SimiLeastSquareEstimation leastSquareEstimation(ptstest, ptsref);
        m = *(leastSquareEstimation.motion());
        break;
      }

      case 3:
      {
        // std::cout<<"Calcul transfo affine..."<<std::endl;
        aims::AffineLeastSquareEstimation leastSquareEstimation(ptstest, ptsref);
        m = *(leastSquareEstimation.motion());
        break;
      }

      default:
      {
        m.setToIdentity();
        throw std::runtime_error("The transformation chosen is not valid !") ;
        break;
      }
    }

    //ROBUSTESSE DEBUT
    
    // Le tri des residus pour ref et test  
    std::multimap<float, Point3df> mt;
    std::multimap<float, Point3df> mr;
    
    std::multimap<float, Point3df>::iterator itt;
    std::multimap<float, Point3df>::iterator itr;
    
    // Diverses variables utiles
    Point3df r; // point residu courant
    float r2; // norme carree de r
    unsigned indice; // indice de boucle multiple (variable locale boucle impossible)
    Motion rigide_kept;
    Motion similitude_kept;

    // La liste des points ref et test ayant les plus petits residus
    std::vector<Point3df> ptstestkept;
    std::vector<Point3df> ptsrefkept;
    
    
    for(int boucle=0; boucle < 3; boucle++)
    {

      //std::cout << "Minimisation - s : " << size << ", m : " << m << std::endl << std::flush;

      mt.clear(); mr.clear();
      ptstestkept.clear(); ptsrefkept.clear();

      for(int l = 0; l < size; l++)  
      {
        Point3df x = ptstest[l];
        Point3df y = ptsref[l];

        r = y - m.transform(x); //calcul du residu pour chaque point
        r2 = r.norm(); // sa norme carree
        mt.insert(std::pair<float, Point3df>(r2, x) ); //insertion du residu et du point test dans map  
        mr.insert(std::pair<float, Point3df>(r2, y) ); //insertion du residu et du point test dans map
      }

      // On remplit les listes avec une partie seulement des meilleurs residus      
      for (itt = mt.begin(),itr = mr.begin(), indice = 0 ;indice<sizekept;++itt,++itr,++indice)
      {
        ptstestkept.push_back( (*itt).second );
        ptsrefkept.push_back( (*itr).second );
      }

      // Calcul de la transformation LS choisie
      switch(_transfo)
      {
        case 1:
        {
          // cout<<"Calcul transfo rigide..."<<cout;
          aims::RigidLeastSquareEstimation leastSquareEstimation_kept(ptstestkept, ptsrefkept);
          m = *(leastSquareEstimation_kept.motion());
          break;
        }

        case 2:
        {
          // cout<<"Calcul transfo similitude..."<<cout;
          aims::SimiLeastSquareEstimation leastSquareEstimation_kept(ptstestkept, ptsrefkept);
          m = *(leastSquareEstimation_kept.motion());
          break;
        }

        case 3:
        {
          // cout<<"Calcul transfo affine..."<<cout;
          aims::AffineLeastSquareEstimation leastSquareEstimation_kept(ptstestkept, ptsrefkept);
          m = *(leastSquareEstimation_kept.motion());
          break;
        }

        default:
        {
          m.setToIdentity();
          throw std::runtime_error("The transformation chosen is not valid !") ;
          break;
        }
      }
    } // fin boucle robustesse          
  }
  else
  {
    cartoDbgMsg( 1, "Not enough information in reference or test image to estimate transformation" );
  }
  
  return (m);
}


#endif