peronaMalikSmoother_d.h

Go to the documentation of this file.

/* This software and supporting documentation are distributed by
 *     Institut Federatif de Recherche 49
 *     CEA/NeuroSpin, Batiment 145,
 *     91191 Gif-sur-Yvette cedex
 *     France
 *
 * This software is governed by the CeCILL-B license under
 * French law and abiding by the rules of distribution of free software.
 * You can  use, modify and/or redistribute the software under the
 * terms of the CeCILL-B license as circulated by CEA, CNRS
 * and INRIA at the following URL "http://www.cecill.info".
 *
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 *
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the
 * same conditions as regards security.
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-B license and that you accept its terms.
 */


#ifndef AIMS_PRIMALSKETCH_PERONAMALIKSMOOTHER_D_H
#define AIMS_PRIMALSKETCH_PERONAMALIKSMOOTHER_D_H

#include <cstdlib>
#include <aims/primalsketch/peronaMalikSmoother.h>
#include <aims/signalfilter/ggradient.h>
#include <aims/signalfilter/gjacobian.h>

namespace aims
{

  template<class T> AimsData<T> PeronaMalikSmoother<T>::doSmoothing(const AimsData<T> & ima, int maxiter, bool /*verbose*/)
  {
    if (maxiter >= 0)
      {
      AimsData<float> kernel(3, 3, 3, 1);
          AimsData<float> *tmp1, *tmp2, *swap, grad;
      int i;
      int x,y,z, count=0;
      int sx=ima.dimX(), sy=ima.dimY(), sz=ima.dimZ();
      std::map<float, int> valueSet;
      std::map<float,int>::iterator itVal;
      float nbThresh;
      float gxm, gxp, gym, gyp, gzm, gzp, cxm, cxp, cym, cyp, czm, czp;

      GaussianGradient<float>  gaussGrad(_sigma, _sigma, _sigma);
      GaussianJacobian<float> gaussJac(_sigma, _sigma, _sigma);
      AimsVector< AimsData< float >, 3 > jaco;

          carto::Converter< AimsData<T>, AimsData<float> > conv;
          AimsData< float > imaF( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() ),
                                    ima2( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
          conv.convert( ima, imaF );


      std::cout << "Evaluating gradient bound : " << std::endl;

      grad=gaussGrad.doit(imaF);
      Writer<AimsData<float> > writerG( "gradient.ima" );
      writerG.write( grad );

      // Evaluation of the gradient parameter from K

      std::cout << "Ordering values" << std::endl;

      nbThresh=sx*sz*sy*_K;

      for (z=0; z<sz; z++)
        for (y=0; y<sy; y++)
          for (x=0; x<sx; x++)
          {
            if (valueSet.find(grad(x,y,z)) != valueSet.end())
              valueSet[grad(x,y,z)]++;
            else
              valueSet.insert(std::pair<float,int>(grad(x,y,z), 1));
          }
      itVal=valueSet.begin();

      std::cout << "Looking for bound with K= " << _K << " at value number " << nbThresh << std::endl;
      for (; itVal!=valueSet.end(); ++itVal)
      {
        count+=(*itVal).second;
        if ((count>=nbThresh) && ((count-1)<nbThresh))
          _gradK=(*itVal).first;
      }
      std::cout << "Gradient bound has been set to " << _gradK << std::endl;

          tmp1=&imaF; tmp2=&ima2;
          std::cout << "Starting " << maxiter << " iterations of diffusion process" << std::endl;

          for (i=0; i<maxiter; i++)
      {
        std::cout << "+" << std::flush;
        grad=gaussGrad.doit((*tmp1));

        for (z=0; z<sz; z++)
          for (y=0; y<sy; y++)
            for (x=0; x<sx; x++)
            {
              if (x==0)
              {
                gxm=0; gxp=(*tmp1)(x+1,y,z) - (*tmp1)(x,y,z);
                cxm=0; cxp=conductance( (grad(x+1,y,z)+grad(x,y,z))/2.0);
              }
              else if (x==sx-1)
              {
                gxm=(*tmp1)(x-1,y,z) - (*tmp1)(x,y,z); gxp=0;
                cxm=conductance( (grad(x-1,y,z)+grad(x,y,z))/2.0); cxp=0;
              }
              else
              {
                gxm=(*tmp1)(x-1,y,z) - (*tmp1)(x,y,z); gxp=(*tmp1)(x+1,y,z) - (*tmp1)(x,y,z);
                cxm=conductance( (grad(x-1,y,z)+grad(x,y,z))/2.0); cxp=conductance( (grad(x+1,y,z)+grad(x,y,z))/2.0);
              }

              if (y==0)
              {
                gym=0; gyp=(*tmp1)(x,y+1,z) - (*tmp1)(x,y,z);
                cym=0; cyp=conductance( (grad(x,y+1,z)+grad(x,y,z))/2.0);
              }
              else if (y==sy-1)
              {
                gym=(*tmp1)(x,y-1,z) - (*tmp1)(x,y,z); gyp=0;
                cym=conductance( (grad(x,y-1,z)+grad(x,y,z))/2.0); cyp=0;
              }
              else
              {
                gym=(*tmp1)(x,y-1,z) - (*tmp1)(x,y,z); gyp=(*tmp1)(x,y+1,z) - (*tmp1)(x,y,z);
                cym=conductance( (grad(x,y-1,z)+grad(x,y,z))/2.0); cyp=conductance( (grad(x,y+1,z)+grad(x,y,z))/2.0);
              }

              if (z==0)
              {
                gzm=0; gzp=(*tmp1)(x,y,z+1) - (*tmp1)(x,y,z);
                czm=0; czp=conductance( (grad(x,y,z+1)+grad(x,y,z))/2.0);
              }
              else if (z==sz-1)
              {
                gzm=(*tmp1)(x,y,z-1) - (*tmp1)(x,y,z); gzp=0;
                czm=conductance( (grad(x,y,z-1)+grad(x,y,z))/2.0); czp=0;
              }
              else
              {
                gzm=(*tmp1)(x,y,z-1) - (*tmp1)(x,y,z); gzp=(*tmp1)(x,y,z+1) - (*tmp1)(x,y,z);
                czm=conductance( (grad(x,y,z-1)+grad(x,y,z))/2.0); czp=conductance( (grad(x,y,z+1)+grad(x,y,z))/2.0);
              }

              (*tmp2)(x,y,z)=(*tmp1)(x,y,z) + _dt *
                                      ( cxm*gxm + cxp*gxp
                                      +cym*gym + cyp*gyp
                                      +czm*gzm + czp*gzp );
            }
        swap=tmp1;
          tmp1=tmp2;
          tmp2=swap;
      }
          std::cout << "Finished" << std::endl;
          carto::Converter< AimsData<float>, AimsData<T> > conv2;
          AimsData<T>  ima3( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
          conv2.convert( (*tmp1), ima3);

          return ima3;
    }
      else
      {
        std::cerr << "PeronaMalik smoother  must have tIn < tOut" << std::endl;
    exit(EXIT_FAILURE);
      }
  }

}

#endif