Normalization.cpp

Go to the documentation of this file.

/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#include<ImLib3D/Normalization.hpp>
#include<ImLib3D/Arithmetic.hpp>
#include<ImLib3D/ImageStats.hpp>
#include<ImLib3D/Signal1D.hpp>
//  #include<ImLib3D/Display.hpp>

#include<ImLib3D/Display.hpp>
Signal1Df 
JointHistogram::ConfidenceWeights() const
{
    Signal1Df res(nbSrcBins,SrcVal(0),SrcVal(nbSrcBins));
    for(int i=0;i<nbSrcBins;i++){res[i]=GetSrcBinTotalProba(i);}
    return res;
}

Signal1Df 
JointHistogram::SmoothAndExtend(const Signal1Df &src,const Signal1Df &weights,double smoothSigma,double outlierSigma)
{
    Signal1Df res=src;
    // first compute reliable but imprecise curve
    // then refine curve, using  outlier downweighting 
    for(int iter=0;iter<2;iter++)
    {
        double sigma;
        if(iter==0){sigma=smoothSigma*5;}
        if(iter==1){sigma=smoothSigma;}
        
        // predicts values from neighboring tangents, using weights
        // and using gaussian smoothing
        for(int x=0;x<res.Size();x++)
        {
            double s=0;
            double totcoef=0;
            for(int x0=0;x0<res.Size()-1;x0++)
            {
                int x1=x0+1;
                if(x0<0 || x1>=res.Size()){continue;}
                // tangent
                double dx=(x-.5*(x0+x1));
                double dy=src[x1]-src[x0];
                if(dy<0){continue;}
                // confidence weights
                double w=min(weights[x0],weights[x1]);
                // gaussian smoothing
                double coef=exp(-dx*dx/(2*sigma*sigma));
                double predictedValue=dy*(float)(x-x0)/(float)(x1-x0) + src[x0];
                // outlier downweighting
                double reweight=1;
                if(iter>0)
                {
                    double error=predictedValue-res[x];
                    reweight=exp(-error*error/(2*outlierSigma*outlierSigma));
                }
                double finalCoef=w*reweight*coef;
                //              printf("*** x0:%d x1:%d coef:%f w:%f predictedValue:%f reweight:%f finalCoef:%f\n",x0,x1,coef,w,predictedValue,reweight,finalCoef);
                s+=finalCoef*predictedValue;
                totcoef+=finalCoef;
            }
            //          printf("*** x:%d s:%f totcoef:%f\n",x,s,totcoef);
            res[x]=(totcoef>0 ? s/totcoef : 0);
        }
        //      res.GnuPlot();
    }

    // impose increasing function, starting at max confidence
    float wmax=-1;
    int wxmax=-1;
    for(int x=0;x<weights.Size();x++){if(weights[x]>wmax){wmax=weights[x];wxmax=x;}}

    for(int x=wxmax+1;x<weights.Size();x++){res[x]=max(res[x],res[x-1]);}
    for(int x=wxmax-1;x>=0            ;x--){res[x]=min(res[x],res[x+1]);}

    // smooth result
    res.Smooth(2);
    res.Smooth(1);

    return res;
}

void
IP3D::NormalizeAverage(const Image3Df& src,const Image3Df& ref,Image3Df& res)
{
    double avRef,varRef;
    IP3D::AverageAndVariance(ref,avRef,varRef);
    double avSrc,varSrc;
    IP3D::AverageAndVariance(src,avSrc,varSrc);
    printf("src:%f\n",avSrc);
    res=src;
    res+=(avRef-avSrc);
}

void
IP3D::NormalizeAverageAndVariance(const Image3Df& src,const Image3Df& ref,Image3Df& res)
{
    double avRef,varRef;
    IP3D::AverageAndVariance(ref,avRef,varRef);
    double avSrc,varSrc;
//  double avRes,varRes;
    IP3D::AverageAndVariance(src,avSrc,varSrc);
    printf("src:%f %f\n",avSrc,varSrc);
    if(varSrc==0){fprintf(stderr,"WARNING: NormalizeAverageAndVariance: sourceimage variance is 0 (image is constant), just normalizing averge\n");varSrc=1;}
    res=src;
    res-=avSrc;
    res*=sqrt(varRef/varSrc);
    res+=avRef;

}

// **************************************************************************************
// **************************************************************************************
// **************************************************************************************

void 
JointHistogram::AddDisplayGrid()
{
    // add "grid" for visu
    if(!histogramImage.HasMask()){histogramImage.AddMask();}
    histogramImage.Mask().Fill(0);
    double ddSrc=pow(10.0,floor(log10(srcMax-srcMin)));
    double ddRef=pow(10.0,floor(log10(refMax-refMin)));
    for(Mask3D::iteratorXYZ p=histogramImage.Mask().begin();p!=histogramImage.Mask().end();++p)
    {
        double xx=srcMin+p.x*(srcMax-srcMin)/(double)nbSrcBins;
        double yy=refMin+p.y*(refMax-refMin)/(double)nbRefBins;
        double xx1=srcMin+(p.x+1)*(srcMax-srcMin)/(double)nbSrcBins;
        double yy1=refMin+(p.y+1)*(refMax-refMin)/(double)nbRefBins;
        if(floor(xx/ddSrc)!=floor(xx1/ddSrc))
        {
            if(p.y==0){mprintf("grid x:%f\n",xx);}
            *p=1;
        }
        if(floor(yy/ddRef)!=floor(yy1/ddRef))
        {
            if(p.x==0){mprintf("grid y:%f\n",yy);}
            *p=1;
        }
    }
}

float 
JointHistogram::GetSrcBinTotalProba(int x) const
{
    float res=0;
    for(int y=0;y<nbRefBins;y++){res+=histogramImage(x,y,0);}
    return res;
}
Signal1Df 
JointHistogram::MedianTransferFct() const 
{
    Signal1Df res(nbSrcBins,SrcVal(0),SrcVal(nbSrcBins));
    for(int x=0;x<nbSrcBins;x++)
    {
        double tot=0;
        int y;
        double avg=0;
        for(y=0;y<nbRefBins;y++)
        {
            tot+=histogramImage(x,y,0);
            avg+=y*histogramImage(x,y,0);
        }
        printf("tot %d:%f\n",x,tot);
        avg/=tot;
        double sum=0;
        int ly=0;
        for(y=0;y<nbRefBins;y++)
        {
            sum+=histogramImage(x,y,0);
            if(sum>tot/1.9999999)
            {
                if(ly==y-1)
                {
                    double ls=sum-histogramImage(x,y,0);
                    double ymiddle=(y-ly)*(0.5-ls/tot)/(sum/tot-ls/tot)+ly;                 
                    printf("tot:%f sum:%f y:%d ly:%d ls:%f ymiddle:%f(%f) avg:%f(%f)\n",tot,sum,y,ly,ls,ymiddle,RefVal(ymiddle),avg,RefVal(avg));
                    res[x]=RefVal(ymiddle);                 
                }
                else
                {
                    double ymiddle=(y+ly)/2.0;
                    printf("tot:%f sum:%f y:%d ly:%d ymiddle:%f\n",tot,sum,y,ly,ymiddle);
                    res[x]=RefVal(ymiddle);
                }
                break;
            }
            if(histogramImage(x,y,0)){ly=y;}
        }
    }
    return res;
}

Signal1Df 
JointHistogram::MaximumProbabilityTransferFct() const
{
    Signal1Df res(nbSrcBins,SrcVal(0),SrcVal(nbSrcBins));
    for(int x=0;x<nbSrcBins;x++)
    {
        Signal1Df proba0(nbRefBins,RefVal(0),RefVal(nbRefBins));
        for(int y=0;y<nbRefBins;y++)
        {
            proba0[y]=histogramImage(x,y,0);
        }
        int supersample=10;
        Signal1Df proba(nbRefBins*supersample,RefVal(0),RefVal(nbRefBins));
        for(int y=0;y<proba.Size();y++)
        {
            proba[y]=proba0(proba.ISampleToX(y));
        }
        proba.Smooth(2*supersample);
        proba.Smooth(2*supersample);
        proba.Smooth(2*supersample);
        double maxv=-1;
        int ymax=-1;
        for(int y=0;y<proba.Size();y++)
        {
            if(proba[y]>maxv){ymax=y;maxv=proba[y];}
        }
        // affine result with quadratic approx (not very usefull, just usefull if supersample is low)
//          if(ymax>0 && ymax<proba.Size()-1)
//          {
//              float dy0=proba[ymax-1]-proba[ymax];
//              float dy1=proba[ymax+1]-proba[ymax];
//              float qymax=.5*((dy0-dy1)/(dy0+dy1))+ymax;
//              float x0=proba.GetSupport().first;
//              float x1=proba.GetSupport().second;
//              res[x]=qymax*(x1-x0)/proba.Size()+x0;
//          }
//          else
        {
            res[x]=proba.ISampleToX(ymax);
        }
    }
    return res;
}


void 
JointHistogram::ComputeFromIndependentSamples(const Image3Df &src,const Image3Df &ref)
{
    // make joint histogram
    for(Image3Df::const_iteratorXYZ p=src.begin();p!=src.end();++p)
    {
        if(src.HasMask() && !src.Mask(p.Pos())){continue;}
        if(ref.HasMask() && !ref.Mask(p.Pos())){continue;}
        histogramImage(SSrcBin(*p),SRefBin(ref(p.Pos())),0)+=1;
    }
}
void 
JointHistogram::ComputeFromDirectionalLinearInterpolation(const Image3Df &src,const Image3Df &ref,int dir)
{
    // make joint histogram
    Vect3Di srcSize=src.SizeV();
    int d0=(dir+1)%3;
    int d1=(dir+2)%3;
    int d2=dir;
//  float minDSrc=.001*(srcMax-srcMin)/(double)nbSrcBins;
//  float minDRef=.001*(refMax-refMin)/(double)nbRefBins;
    Vect3Di dpos(0,0,0);dpos(d2)=1;
    Vect3Di pos(0,0,0);
    for(int x0=0;x0<srcSize(d0);x0++)
    {
        pos(d0)=x0;
        for(int x1=0;x1<srcSize(d1);x1++)
        {
            pos(d1)=x1;
            for(int x2=0;x2<srcSize(d2)-1;x2++)
            {
                pos(d2)=x2;
                Vect3Di pos0=pos;
                Vect3Di pos1=pos+dpos;
                // intenisty of src from pos0 to pos1 goes from src0 to src1
                // while intensity of ref goes from ref0 to ref1

                // skip if we're on unmasked
                if( (src.HasMask() && ( !src.Mask(pos0) || !src.Mask(pos1))) ||
                    (ref.HasMask() && ( !ref.Mask(pos0) || !ref.Mask(pos1)))   )  {continue;}
                float src0=src(pos0);
                float ref0=ref(pos0);
                float src1=src(pos1);
                float ref1=ref(pos1);
//              pos0.Show("pos0:","   ");
//              pos1.Show("pos1:","\n");
//              printf("src0:%f ref0:%f    src1:%f ref1:%f\n",src0,ref0,src1,ref1);
                // draw line of value 1/linesize in histogramImage between these two positions
                {
                    Vect3Df p0(SrcBinF(src0),RefBinF(ref0),0);
                    Vect3Df p1(SrcBinF(src1),RefBinF(ref1),0);
                    Vect3Df dP=p1-p0;
                    int dir=dP.MaxAbsCoord();
                    bool incr=dP(dir)>0;
                    dP/=fabs(dP(dir));
                    int linesize=0;
                    for(Vect3Df p=p0;(incr? p(dir)<p1(dir) : p(dir)>p1(dir));p+=dP)
                    {
                        linesize++;
                    }
                    for(Vect3Df p=p0;(incr? p(dir)<p1(dir) : p(dir)>p1(dir));p+=dP)
                    {
                        if(histogramImage.IsInside(rint(p))){histogramImage(rint(p))+=1.0/linesize;}
                    }
                }
            }
        }
    }
}

JointHistogram::JointHistogram(const string &filename):
    histogramImage(filename)
{
    nbSrcBins=histogramImage.Property<int>("nbSrcBins");
    nbRefBins=histogramImage.Property<int>("nbRefBins");
    srcMin=histogramImage.Property<double>("srcMin");
    srcMax=histogramImage.Property<double>("srcMax");
    refMin=histogramImage.Property<double>("refMin");
    refMax=histogramImage.Property<double>("refMax");
    printf("nbSrcBins:%d nbRefBins:%d srcMin:%f srcMax:%f refMin:%f refMax:%f\n",
           nbSrcBins, nbRefBins, srcMin, srcMax, refMin, refMax);
}
JointHistogram::JointHistogram(const Image3Df &src,const Image3Df &ref,int _nbSrcBins,int _nbRefBins):
    nbSrcBins(_nbSrcBins),nbRefBins(_nbRefBins),histogramImage(_nbSrcBins,_nbRefBins,1)     
{
    histogramImage.Fill(0);
    // find ranges/ bounds
    FindRanges(src,ref);

//      ComputeFromIndependentSamples(src,ref);
    ComputeFromDirectionalLinearInterpolation(src,ref,0);
    ComputeFromDirectionalLinearInterpolation(src,ref,1);
    ComputeFromDirectionalLinearInterpolation(src,ref,2);
    histogramImage*=1/3.0;
}

void 
JointHistogram::FindRanges(const Image3Df& src, const Image3Df &ref)
{
    vector<float> srcValues;
    vector<float> refValues;
    srcValues.reserve(src.GetNVoxels());
    refValues.reserve(ref.GetNVoxels());
    for(Image3Df::const_iteratorXYZ p=src.begin();p!=src.end();++p)
    {
        if(src.HasMask() && !src.Mask(p.Pos())){continue;}
        if(ref.HasMask() && !ref.Mask(p.Pos())){continue;}
        srcValues.push_back(*p);
        refValues.push_back(ref(p.Pos()));
    }

    // ProbabilityDistribution uses bins, which might not be reliable enough 

    // use nth_element to find histogram range (this is much faster than sorting)
    vector<float>::iterator pnth;
    pnth=srcValues.begin()+(size_t)floor(.005*srcValues.size());
    nth_element(srcValues.begin(),pnth,srcValues.end());
    srcMin=*pnth;
    pnth=srcValues.begin()+min(srcValues.size()-1,(size_t)ceil(.995*srcValues.size()));
    nth_element(srcValues.begin(),pnth,srcValues.end());
    srcMax=*pnth;
    float d=srcMax-srcMin;
    srcMin-=d*.3;
    srcMax+=d*.3;

    pnth=refValues.begin()+(size_t)floor(.005*refValues.size());
    nth_element(refValues.begin(),pnth,refValues.end());
    refMin=*pnth;
    pnth=refValues.begin()+min(refValues.size()-1,(size_t)ceil(.995*refValues.size()));
    nth_element(refValues.begin(),pnth,refValues.end());
    refMax=*pnth;
    d=refMax-refMin;
    refMin-=d*.6;
    refMax+=d*.6;

    histogramImage.AddProperty("srcMin",srcMin);
    histogramImage.AddProperty("srcMax",srcMax);
    histogramImage.AddProperty("refMin",refMin);
    histogramImage.AddProperty("refMax",refMax);
    histogramImage.AddProperty("nbSrcBins",nbSrcBins);
    histogramImage.AddProperty("nbRefBins",nbRefBins);

    histogramImage.AddSupport(Vect3Df(srcMin,refMin,0),Vect3Df(srcMax,refMax,1));
    mprintf("JointHistogram srcMin:%f srcMax:%f refMin:%f refMax:%f\n",srcMin,srcMax,refMin,refMax);
}

// **************************************************************************************
// **************************************************************************************
void
IP3D::NormalizeJointHistogram(const Image3Df& src0,const Image3Df& ref,Image3Df& res,int nbSrcBins,int nbRefBins)
{
    Image3Df src;
    IP3D::NormalizeAverageAndVariance(src0,ref,src);
    // compute joint histo
    JointHistogram joint(src,ref,nbSrcBins,nbRefBins);
    // compute transfer fct
    Signal1Df transferFct=JointHistogram::SmoothAndExtend(joint.MaximumProbabilityTransferFct(),
                                                          joint.ConfidenceWeights(),
                                                          joint.nbSrcBins/50.0,
                                                          (joint.refMax-joint.refMin)/25);

    // apply transfer fct
    // single transfer fct on src image only
    res=src;
    for(Image3Df::iteratorFast p=res.begin();p!=res.end();++p)
    {
        *p=transferFct(*p);
    }
}



// **************************************************************************************
void
IP3D::ComputeJointHistogramImage(const Image3Df& src0,const Image3Df& ref,Image3Df& res,int drawTransferFunction,int nbSrcBins,int nbRefBins)
{
    Image3Df src;
    IP3D::NormalizeAverageAndVariance(src0,ref,src);
    JointHistogram joint(src,ref,nbSrcBins,nbRefBins);

    Signal1Df transferFct;
    if(drawTransferFunction==1)
    {
        transferFct=joint.BestFitTransferFct();
    }
    if(drawTransferFunction==2)
    {
        transferFct=JointHistogram::SmoothAndExtend(joint.MaximumProbabilityTransferFct(),
                                                              joint.ConfidenceWeights(),
                                                              joint.nbSrcBins/50.0,
                                                              (joint.refMax-joint.refMin)/25);
    }
    if(drawTransferFunction>0)
    {
        for(int x=0;x<joint.histogramImage.Width();x++)
        {
            joint.histogramImage(x,joint.SRefBin(transferFct(joint.SrcVal(x))),0)=-10;
        }
    }
    res=joint.histogramImage;
}


// **************************************************************************************
// **************************************************************************************
// **************************************************************************************
// **************************************************************************************
// #include<gsl/gsl_version.h>
// #if GSL_VERSION="1.1.1"
// Signal1Df 
// JointHistogram::BestFitTransferFct(double smoothingCoef,bool allowDecreasing) const
// {
//  ThrowError("JointHistogram::BestFitTransferFct requires at least gsl 1.2 ... not compiled");
// }
// #else
#include<gsl/gsl_multimin.h>

class IntensityTransferFctFitter
{
    const JointHistogram &joint;
    double smoothingCoef,decreasingCoef;
    Signal1Df central;
    Signal1Df confidenceWeights;
    double totalWeight;
    double refImageVariance;
public:
    void SetDecreasingCoef(double _decreasingCoef){decreasingCoef=_decreasingCoef;}
    //  transfer function is linear by parts
    float TransferFctValue(const vector<double> &values,int xpos)
    {
        // just test with linear
        float pos=values.size()*xpos/(float)joint.nbSrcBins;
        int p0=min((int)values.size()-2,max(0,(int)floor(pos)));
        int p1=p0+1;
        return (pos-p0)*values[p1]+(p1-pos)*values[p0];
    }
    void Filter(const vector<double> &src,const vector<double> &filter,vector<double> &res)
    {
        if((filter.size()%2)==0){ThrowError("filter size must bue unpair");}
        res.resize(src.size());
        for(uint i=0;i<src.size();i++)
        {
            int center=(filter.size()-1)/2;
            // i+j-center>=0  => j>=center-i
            int jmin=max(center-(int)i,0);
            // i+j-center<=s-1  => j<s-1+center-i
            int jmax=min(((int)src.size())-1+center-(int)i,(int)filter.size()-1);               
            double total=0;
            double totalc=0;
            for(int j=jmin;j<=jmax;j++)
            {
                total+=filter[j]*src[i+j-center];
                totalc+=filter[j];
            }
            res[i]=total/totalc;
        }
    }
    // error of fit
    double TotalError(const vector<double> &values)
    {
        bool talk=false;//(rand()%000)==0;
        // precompute transfer fct values
        vector<double> splineYPos(joint.nbSrcBins);
        for(int x=0;x<joint.nbSrcBins;x++)
        {
            splineYPos[x]=TransferFctValue(values,x);
        }
        // compute error due to joint histogram
        double errorJointHistogram=0;
        vector<float> debugDisplayImageError;
        for(int x=0;x<joint.nbSrcBins;x++)
        {
            double d=joint.RefVal(splineYPos[x])-central[x];
            errorJointHistogram+=confidenceWeights[x]*(d*d);
            if(talk){debugDisplayImageError.push_back(confidenceWeights[x]*(d*d));}
        }
        if(talk){GnuPlot(debugDisplayImageError,"image error");}
        errorJointHistogram/=refImageVariance*totalWeight;
        // compute smoothness error 
        double errorSmoothness=0;
        {
            // trasnfer fct is a linear interpolated: (delta convoluted by 1st degree spline)
            // second derivative of 1st degree spline is
            vector<double> secondDeriv(values.size()-2);
            for(uint i=1;i<values.size()-1;i++)
            {
                secondDeriv[i-1]=values[i+1]-2*values[i]+values[i-1];
            }
            // smooth out secondDeriv
            vector<double> smoothed0;
            {
                double filter00[] =  {1,3,4,3,1};
                vector<double> filter0(filter00, filter00+sizeof(filter00)/sizeof(double));
                Filter(secondDeriv,filter0,smoothed0);
            }
            // smooth out secondDeriv
            vector<double> smoothed1;
            {
                double filter10[] = {1,2,3,4,5,6,5,4,3,2,1};
                vector<double> filter1(filter10, filter10+sizeof(filter10)/sizeof(double));
                Filter(secondDeriv,filter1,smoothed1);
            }

            if(talk)
            {
                GnuPlot(values,"transfer function values");
                GnuPlot(secondDeriv,"secondDeriv");
                GnuPlot(smoothed0,"smoothed0");
                GnuPlot(smoothed1,"smoothed1");
                printf("print\n");
            }
            // we dont want largeScaleSmoothnessCoef to be very strong, because
            // it introduces a bias on correct results
            // now compute error
            for(uint i=0;i<secondDeriv.size();i++)
            {
                errorSmoothness+=
                    //                  secondDeriv[i]*secondDeriv[i]*.0005+
                    smoothed0  [i]*smoothed0[i]+
                    smoothed1  [i]*smoothed1[i]*.2;
                    
            }
        }
        errorSmoothness*=40.0/(double)(joint.nbRefBins*joint.nbRefBins*values.size());
        // compute decreasing error 
        double errorDecreasing=0;
        for(uint i=0;i<values.size()-1;i++)
        {
            errorDecreasing+=max(0.0,-(values[i+1]-values[i]));
        }
        errorDecreasing*=.0001;
        return errorJointHistogram+smoothingCoef*errorSmoothness+decreasingCoef*errorDecreasing;
    }

    IntensityTransferFctFitter(const JointHistogram &_joint,double _smoothingCoef=1,double _decreasingCoef=0):
        joint(_joint),
        smoothingCoef(_smoothingCoef),
        decreasingCoef(_decreasingCoef)
    {
        central=joint.MaximumProbabilityTransferFct();
        confidenceWeights=joint.ConfidenceWeights();
        // debug : display central 
//          {
//              Image3Df tmp=joint.histogramImage;
//              for(int x=0;x<tmp.Width();x++)
//              {
//                  tmp(x,joint.SRefBin(central[x]),0)=-10;
//                  printf("%f:%f ",central.ISampleToX(x),central[x]);
//              }       
//              printf("\n");
//              tmp.WriteToFile("central.im3D");
//          }
        // ref image variance
        refImageVariance=0;
        totalWeight=0;
        double refImageAverage=0;
        for(Image3Df::const_iteratorXYZ p=joint.histogramImage.begin();p!=joint.histogramImage.end();++p)
        {
            totalWeight+=*p;
            refImageVariance+=(*p)*joint.RefVal(p.y)*joint.RefVal(p.y);
            refImageAverage+=(*p)*joint.RefVal(p.y);
        }
        refImageAverage/=totalWeight;
        refImageVariance=(refImageVariance/totalWeight)-refImageAverage;
        printf("refImageAverage:%f refImageVariance:%f stddev:%f\n",refImageAverage,refImageVariance,sqrt(refImageVariance));
    }

};
vector<double>
MakeSTLVector(const gsl_vector *v)
{
    vector<double> values(v->size);
    for(uint i=0;i<values.size();i++){values[i]=gsl_vector_get(v, i);}
    return values;
}
gsl_vector *
MakeGSLVector(const vector<double> &src)
{
    gsl_vector *res = gsl_vector_alloc (src.size());  
    for(uint i=0;i<src.size();i++)
    {
        gsl_vector_set (res, i, src[i]);
    }
    return res;
}

double
GSLMinimizationFunction (const gsl_vector *v, void *params)
{
    return ((IntensityTransferFctFitter *)params)->TotalError(MakeSTLVector(v));
}

int
FitIntensityTransferFctUsingGSLMinimization(IntensityTransferFctFitter &intensityTransferFctFitter,vector<double> &values0)
{
    gsl_vector *values=MakeGSLVector(values0);

    // Initial vertex size vector  (??)
    gsl_vector *ss;
    ss = gsl_vector_alloc (values->size);
    gsl_vector_set_all (ss, 1.0);


    // initialize function
    gsl_multimin_function minex_func;
    {
        minex_func.f = &GSLMinimizationFunction;
        minex_func.n = values->size;
        minex_func.params = (void *)&intensityTransferFctFitter;
    }

    // initialize minimizer
    gsl_multimin_fminimizer *minimizer = 
        minimizer = gsl_multimin_fminimizer_alloc (gsl_multimin_fminimizer_nmsimplex, values->size);
    gsl_multimin_fminimizer_set (minimizer, &minex_func, values, ss);


    // iterate minimizations
    int rval = GSL_CONTINUE;
    int status = GSL_SUCCESS;
    double ssval;
    size_t iter = 0;
    while (rval == GSL_CONTINUE)
    {
        iter++;
        status = gsl_multimin_fminimizer_iterate(minimizer);
      
        if (status) 
        break;

        rval = gsl_multimin_test_size (gsl_multimin_fminimizer_size (minimizer),
                                       1e-2);
        ssval = gsl_multimin_fminimizer_size (minimizer);

        if (rval == GSL_SUCCESS) 
        {
            printf ("converged to minimum at\n");

            printf("iter %4d : ", iter);
            for (uint i = 0; i < values->size; i++)
            {
            printf ("%5.2f ", gsl_vector_get (minimizer->x, i));
            }
            printf("f() = %-10.10f ssize = %.3f\n", minimizer->fval, ssval);
        }
    }
    // cleanup 
    values0=MakeSTLVector(minimizer->x);
    gsl_vector_free(values);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free (minimizer);
    return status;
}

Signal1Df 
JointHistogram::BestFitTransferFct(double smoothingCoef,bool allowDecreasing) const
{
    IntensityTransferFctFitter intensityTransferFctFitter(*this);

    // Starting point (identity transfer fct)
    vector<double> values(40);
    for(uint i=0;i<values.size();i++)
    {
        values[i]=(nbRefBins*i)/(float)values.size();
    }
    // if decreasing is not allow. iteratively fit
    // while progressively increasing decreasingCoef
    for(double decreasingCoef=.0001;decreasingCoef<1000000;decreasingCoef*=4)
    {
        if(!allowDecreasing){intensityTransferFctFitter.SetDecreasingCoef(decreasingCoef);}
        FitIntensityTransferFctUsingGSLMinimization(intensityTransferFctFitter,values);
        // check if decreasing
        bool decreasing=false;
        for(uint i=0;i<values.size()-1;i++)
        {
            if(values[i+1]<values[i]){decreasing=true;break;}
        }
        if(allowDecreasing || !decreasing){break;}
    }
    // compute final results from fitted values
    Signal1Df res(nbSrcBins,SrcVal(0),SrcVal(nbSrcBins));
    for(int i=0;i<res.Size();i++)
    {
        res[i]=RefVal(intensityTransferFctFitter.TransferFctValue(values,i));
    }
    return res;
}

---