Algorithme de Levenberg Marquardt

**TopCao** · 06/01/2010, 22h26

Bonjour à tous,

Je tiens tout d'abord à préciser que je suis débutant et totalement autodidact en C++.

Alors je dispose d'une série de données (48 valeurs) dont je sais qu'elles correspondent à une fonction du type :

S(b) = Sa.exp(-b.ADCa) + Sb.exp(-b.ADCb) + Sc.exp(-b.ADCc) + bruit

S étant le signal et b ma variable.

Je voudrais donc à l'aide d'un algorithme de Levenberg Marquardt retrouver les coefficients Sa, Sb, Sc, ADCa, ADCb et ADCc.

Je voulais mettre mes valeurs dans un tableau de pointeurs et générer mes 6 coeff mais je n'y arrive pas.

Je suis donc preneur de toute aide et je peux envoyer les bouts de code que j'ai commencé à taper.

Précisions : je prog en c++ avec Xcode (sous Mac) pour faire des progs en ligne de commande.

Merci à vous

alex

**Klaim** · 07/01/2010, 02h13

Il faudrait que tu montres le code que tu as déjà pour que tout ça soit bien clair stp. Comme ça on peut te guider.

**TopCao** · 07/01/2010, 18h17

Bonjour,

Voici mon code
3 fichiers :

itkCompex.h : mon header
itkComputeCompex.cxx : mon prog
itkCompex.txx : le template qui contient l'algo

Pour bien comprendre :
J'écris ce programme dans le but de l'exécuter sur une série d'images (d'ou la structure ITK rencontrée partout) mais je souhaite tester la fiabilité de l'algo de Levenberg sur une série de chiffres (en fait la valeur d'un voxel au cours du temps)

J'espère que vous n'allez pas être pris de convulsions en regardant mon code qui je l'admets est vraiment cracra et pas compilable en l'état !!

Merci pour vos conseils

Alex

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
/*
 *  itkCompex.h
 *  compex
 *
 *  Created by Alexandre Ladoux on 07/01/10.
 *  Copyright 2010 __LIMEC__. All rights reserved.
 *
 */
 
 
#ifndef _itk_Compex_h_
#define _itk_Compex_h_
 
#include <itkImageToImageFilter.h>
 
namespace itk
{
 
	template <class TInputImage, class TOutputImage>
    class ITK_EXPORT Compex :
	public ImageToImageFilter<TInputImage, TOutputImage>
	{
 
	public:
 
		typedef Compex			Self;
		typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
		typedef SmartPointer<Self> Pointer;
		typedef SmartPointer<const Self> ConstPointer;
 
		itkNewMacro  (Self);
		itkTypeMacro (Compex, ImageToImagFilter);
 
		typedef TInputImage  InputImageType;
		typedef TOutputImage OutputImageType;
 
		typedef typename InputImageType::PixelType   InputPixelType;
		typedef typename OutputImageType::PixelType  OutputPixelType;
		typedef typename OutputImageType::RegionType OutputImageRegionType;
		typedef OutputPixelType                      ScalarType;
 
		// Dé-commenter lors de l'application du process à des images
		// itkSetObjectMacro (DataSetImages, InputImageType);
		// itkGetObjectMacro (DataSetImages, InputImageType);
 
 
		itkSetMacro (v01, double);
		itkGetMacro (v01, double);
		itkSetMacro (v02, double);
		itkGetMacro (v02, double);
		itkSetMacro (v03, double);
		itkGetMacro (v03, double);
		itkSetMacro (v04, double);
		itkGetMacro (v04, double);
		itkSetMacro (v05, double);
		itkGetMacro (v05, double);
		itkSetMacro (v06, double);
		itkGetMacro (v06, double);
		itkSetMacro (v07, double);
		itkGetMacro (v07, double);
		itkSetMacro (v08, double);
		itkGetMacro (v08, double);
		itkSetMacro (v09, double);
		itkGetMacro (v09, double);
		itkSetMacro (v10, double);
		itkGetMacro (v10, double);
		itkSetMacro (v11, double);
		itkGetMacro (v11, double);
		itkSetMacro (v12, double);
		itkGetMacro (v12, double);
		itkSetMacro (v13, double);
		itkGetMacro (v13, double);
		itkSetMacro (v14, double);
		itkGetMacro (v14, double);		
		itkSetMacro (v15, double);
		itkGetMacro (v15, double);
		itkSetMacro (v16, double);
		itkGetMacro (v16, double);
		itkSetMacro (v17, double);
		itkGetMacro (v17, double);
		itkSetMacro (v18, double);
		itkGetMacro (v18, double);
		itkSetMacro (v19, double);
		itkGetMacro (v19, double);
		itkSetMacro (v20, double);
		itkGetMacro (v20, double);
		itkSetMacro (v21, double);
		itkGetMacro (v21, double);
		itkSetMacro (v22, double);
		itkGetMacro (v22, double);
		itkSetMacro (v23, double);
		itkGetMacro (v23, double);
		itkSetMacro (v24, double);
		itkGetMacro (v24, double);
		itkSetMacro (v25, double);
		itkGetMacro (v25, double);
		itkSetMacro (v26, double);
		itkGetMacro (v26, double);
		itkSetMacro (v27, double);
		itkGetMacro (v27, double);
		itkSetMacro (v28, double);
		itkGetMacro (v28, double);				
		itkSetMacro (v29, double);
		itkGetMacro (v29, double);
		itkSetMacro (v30, double);
		itkGetMacro (v30, double);
		itkSetMacro (v31, double);
		itkGetMacro (v31, double);
		itkSetMacro (v32, double);
		itkGetMacro (v32, double);
		itkSetMacro (v33, double);
		itkGetMacro (v33, double);
		itkSetMacro (v34, double);
		itkGetMacro (v34, double);
		itkSetMacro (v35, double);
		itkGetMacro (v35, double);
		itkSetMacro (v36, double);
		itkGetMacro (v36, double);
		itkSetMacro (v37, double);
		itkGetMacro (v37, double);
		itkSetMacro (v38, double);
		itkGetMacro (v38, double);
		itkSetMacro (v39, double);
		itkGetMacro (v39, double);
		itkSetMacro (v40, double);
		itkGetMacro (v40, double);
		itkSetMacro (v41, double);
		itkGetMacro (v41, double);
		itkSetMacro (v42, double);
		itkGetMacro (v42, double);		
		itkSetMacro (v43, double);
		itkGetMacro (v43, double);
		itkSetMacro (v44, double);
		itkGetMacro (v44, double);
		itkSetMacro (v45, double);
		itkGetMacro (v45, double);
		itkSetMacro (v46, double);
		itkGetMacro (v46, double);
		itkSetMacro (v47, double);
		itkGetMacro (v47, double);
		itkSetMacro (v48, double);
		itkGetMacro (v48, double);
 
 
	protected:
		Compex();
		~Compex(){};
 
 
		void BeforeThreadedGenerateData(void);
		void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, int threadId);
		void PrintSelf(std::ostream& os, Indent indent) const
		{
			Superclass::PrintSelf(os,indent);
		}
 
 
 
	private:
		CompexFilter (const Self&);
		void operator=(const Self&);
 
		double                  m_v01;
		double					m_v02;
		double                  m_v03;
		double					m_v04;
		double                  m_v05;
		double					m_v06;
		double                  m_v07;
		double					m_v08;
		double                  m_v09;
		double					m_v10;
		double                  m_v11;
		double					m_v12;
		double                  m_v13;
		double					m_v14;
		double                  m_v15;
		double					m_v16;
		double                  m_v17;
		double					m_v18;
		double                  m_v19;
		double					m_v20;
		double                  m_v21;
		double					m_v22;
		double                  m_v23;
		double					m_v24;
		double                  m_v25;
		double					m_v26;
		double                  m_v27;
		double					m_v28;
		double                  m_v29;
		double					m_v30;
		double                  m_v31;
		double					m_v32;
		double                  m_v33;
		double					m_v34;
		double                  m_v35;
		double					m_v36;
		double                  m_v37;
		double					m_v38;
		double                  m_v39;
		double					m_v40;
		double                  m_v41;
		double					m_v42;
		double                  m_v43;
		double					m_v44;
		double                  m_v45;
		double					m_v46;
		double                  m_v47;
		double					m_v48;
 
 
		// Dé-commenter lors de l'application du process à des images
		//typename InputImageType::Pointer m_B0Image;
 
	};
 
} // end of namespace
 
#include "itkCompex.txx"
 
#endif

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
// ATTENTION : COMMENTAIRE POUR MOI MEME
//
// Copie du template itkCompex.txx créé uniquement pour avoir l'indentation xcode
// non compilé par Xcode
// Copier le corps du fichier dans itkCompex.txx
// après toute modification et avant toute compilation
 
 
#ifndef _itk_Compex_txx_
#define _itk_Compex_txx_
 
#include <itkCompex.h>
 
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
 
namespace itk
{
 
	template<class TInputImage, class TOutputImage>
	CompexFilter<TInputImage,TOutputImage>
	::CompexFilter()
	{
		m_v01  = 0.0;
		m_v02  = 0.0;
		m_v03  = 0.0;
		m_v04  = 0.0;
		m_v05  = 0.0;
		m_v06  = 0.0;
		m_v07  = 0.0;
		m_v08  = 0.0;
		m_v09  = 0.0;
		m_v10  = 0.0;
		m_v11  = 0.0;
		m_v12  = 0.0;
		m_v13  = 0.0;
		m_v14  = 0.0;
		m_v15  = 0.0;
		m_v16  = 0.0;
		m_v17  = 0.0;
		m_v18  = 0.0;
		m_v19  = 0.0;
		m_v20  = 0.0;
		m_v21  = 0.0;
		m_v22  = 0.0;
		m_v23  = 0.0;
		m_v24  = 0.0;
		m_v25  = 0.0;
		m_v26  = 0.0;
		m_v27  = 0.0;
		m_v28  = 0.0;
		m_v29  = 0.0;
		m_v30  = 0.0;
		m_v31  = 0.0;
		m_v32  = 0.0;
		m_v33  = 0.0;
		m_v34  = 0.0;
		m_v35  = 0.0;
		m_v36  = 0.0;
		m_v37  = 0.0;
		m_v38  = 0.0;
		m_v39  = 0.0;
		m_v40  = 0.0;
		m_v41  = 0.0;
		m_v42  = 0.0;
		m_v43  = 0.0;
		m_v44  = 0.0;
		m_v45  = 0.0;
		m_v46  = 0.0;
		m_v47  = 0.0;
		m_v48  = 0.0;
 
		// Dé-commenter lors de l'application du process à des images
		// m_B0Image = 0;
	}
 
 
	template<class TInputImage, class TOutputImage>
	void
	CompexFilter<TInputImage,TOutputImage>
	::BeforeThreadedGenerateData()
	{
		if( this->GetNumberOfInputs()==0 )
		{
			throw itk::ExceptionObject(__FILE__,__LINE__,"Error: Input is not set.");
		}
 
		if( m_B0Image.IsNull() )
		{
			throw itk::ExceptionObject(__FILE__,__LINE__,"Error: B0 image is not set.");
		}    
	}
 
 
	template<class TInputImage, class TOutputImage>
	void
	CompexFilter<TInputImage,TOutputImage>
	::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, int threadId)
	{
 
		typedef ImageRegionIterator<OutputImageType>      IteratorOutputType;
		typedef ImageRegionConstIterator<InputImageType>  IteratorInputType;
 
 
		IteratorOutputType  itOut (this->GetOutput(),  outputRegionForThread);
		IteratorInputType   itIn  (this->GetInput(),   outputRegionForThread);
		//IteratorInputType   itB0  (this->GetB0Image(), outputRegionForThread);
 
 
		double v01 = this->GetBValue();
		double v02 = this->GetBValue();
		double v03 = this->GetBValue();
		double v04 = this->GetBValue();
		double v05 = this->GetBValue();
		double v06 = this->GetBValue();
		double v07 = this->GetBValue();
		double v08 = this->GetBValue();
		double v09 = this->GetBValue();
		double v10 = this->GetBValue();
		double v11 = this->GetBValue();
		double v12 = this->GetBValue();
		double v13 = this->GetBValue();
		double v14 = this->GetBValue();
		double v15 = this->GetBValue();
		double v16 = this->GetBValue();
		double v17 = this->GetBValue();
		double v18 = this->GetBValue();
		double v19 = this->GetBValue();
		double v20 = this->GetBValue();
		double v21 = this->GetBValue();
		double v22 = this->GetBValue();
		double v23 = this->GetBValue();
		double v24 = this->GetBValue();
		double v25 = this->GetBValue();
		double v26 = this->GetBValue();
		double v27 = this->GetBValue();
		double v28 = this->GetBValue();
		double v29 = this->GetBValue();
		double v30 = this->GetBValue();
		double v31 = this->GetBValue();
		double v32 = this->GetBValue();
		double v33 = this->GetBValue();
		double v34 = this->GetBValue();
		double v35 = this->GetBValue();
		double v36 = this->GetBValue();
		double v37 = this->GetBValue();
		double v38 = this->GetBValue();
		double v39 = this->GetBValue();
		double v40 = this->GetBValue();
		double v41 = this->GetBValue();
		double v42 = this->GetBValue();
		double v43 = this->GetBValue();
		double v44 = this->GetBValue();
		double v45 = this->GetBValue();
		double v46 = this->GetBValue();
		double v47 = this->GetBValue();
		double v48 = this->GetBValue();
 
 
 
// BOUCLE DE L'ALGORITHME DE LEVENBERG MARQUARDT		
 
 
		// Paramètres des données
		// int xRange = 128     ;	//Résolution en x quand image
		// int yRange = 128     ;	//Résolution en y quand image
		// int nSlices = 8      ;	//Résolution en z quand imageS
		int nbbValues = 48   ;  //Résolution en z
		int tableau[bValues] ;  //Tableau contenant les valeurs de b pour un voxel donné
 
 
		// Paramètres d'initialisation
		double Sa0  = 125	  ;
		double Sb0  = 99,001  ;
		double Sc0  = 77,28	  ;
		double ADCa = 0,0001  ;
		double ADCb = 0,00002 ;
		double ADCc = 0,00003 ;
 
		// Paramètres estimés
		double estSa0  = Sa0  ;
		double estSb0  = Sb0  ;
		double estSc0  = Sc0  ;
		double estADCa = ADCa ;
		double estADCb = ADCb ;
		double estADCc = ADCc ;
 
		// Fonction d'initialisation
		double estSignal = Sa0 * exp(-ADCa * bValue) + Sb0 * exp(-ADCb * bValue) +Sc0 * exp(-ADCc * bValue);
		int nbIter = 100;				// Nombre d'irérations
		double lambda = 0,001;          // Valeur initiale du facteur d'amortissement
		int nbParams = 6;				// Nombre de paramètres à estimer
 
 
		// Initialisation
		for (iter = 1; iter = nbIter; iter++) 
		{
 
			{
 
				tableau[0] = v01 ;				// tableau contenant à la suite les n valeurs de b
				...
				tableau[47] = v48 ;
 
 
				// ????? Evaluer la matrice jacobienne pour les paramètres estimés.
				double J = ???;
 
				// Faire l'erreur de distance pour les paramètres estimés
				estSignal = Sa0 * exp(-ADCa * bValue) + Sb0 * exp(-ADCb * bValue) +Sc0 * exp(-ADCc * bValue);
 
				// Calculer la matrice d'Hessian approximée avec J' la trasposée de J
				H = J' * J ;
 
				// appliquer le factor d'amortissement à la matrice hessienne
				double matriceLM = H + (lambda * nbParams);
 
				//incrémenter tous les paramètres
				double inc = - inv(matriceLM) * (J'*k);
 
 
 
			}		
 
		}
 
 
#endif

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#include "GetPot.h"
 
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkCompex.h>
 
//#include "itkGISImageIOFactory.h"
 
void PrintHelp (const char* exec)
{
	std::cout << "Usage: " << std::endl;
	std::cout << exec << "	<-v01 value01> 
	<-v02 value02> 
	<-v03 value03> 
	<-v04 value04> 
	<-v05 value06> 
	<-v06 value06> 
	<-v07 value07>
	<-v08 value08> 
	<-v09 value09> 
	<-v10 value10> 
	<-v11 value11> 
	<-v12 value12> 
	<-v13 value13>
	<-v14 value14> 
	<-v15 value15> 
	<-v16 value16> 
	<-v17 value17> 
	<-v18 value18>	
	<-v19 value19>
	<-v20 value20> 
	<-v21 value21> 
	<-v22 value22> 
	<-v23 value23> 
	<-v24 value24> 
	<-v25 value25>	
	<-v26 value26> 
	<-v27 value27> 
	<-v28 value28>
	<-v29 value29> 
	<-v30 value30> 
	<-v31 value31> 
	<-v32 value32> 
	<-v33 value33>	
	<-v34 value34>
	<-v35 value35> 
	<-v36 value36> 
	<-v37 value37> 
	<-v38 value38> 
	<-v39 value39> 
	<-v40 value40>	
	<-v41 value41> 
	<-v42 value42> 
	<-v43 value43> 
	<-v44 value44>
	<-v45 value45> 
	<-v46 value46> 
	<-v47 value47> 
	<-v48 value48> 
 
 
 
	" << std::endl;
}
 
 
int main (int narg, char *arg[])
{
 
	// parse arguments
	GetPot cl (narg, arg);
	if( cl.size() == 1 || cl.search (2,"--help","-h") ) {
		PrintHelp (cl[0]);
		return -1;
	}
 
 
	// scanner les parametres
	double bvalue = cl.follow (1000.0, 2, "-b", "-B");
 
	const char* input   = cl.follow ("", 2, "-i", "-I");
	const char* b0image = cl.follow ("", 2, "-b0", "-B0");
	const char* output  = cl.follow ("", 2, "-o", "-O");
 
	//itk::GISImageIOFactory::RegisterOneFactory();
 
 
	// lire l'image donnee en parametre
	typedef itk::Image<double, 3>           ImageType;
	typedef itk::ImageFileReader<ImageType> ReaderType;
 
	ImageType::Pointer InputImage = 0;
	{
		ReaderType::Pointer reader = ReaderType::New();
		reader->SetFileName ( input );
		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cerr << e;
			return -1;
		}
 
		InputImage = reader->GetOutput(); 
	}
 
 
 
	ImageType::Pointer B0Image = 0;
	{
		ReaderType::Pointer reader = ReaderType::New();
		reader->SetFileName ( b0image );
		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cerr << e;
			return -1;
		}
 
		B0Image = reader->GetOutput();
	}
 
 
 
	typedef itk::DiffusionWeightedImageToADCImageFilter<ImageType, ImageType> ADCCalculatorType;
 
	ADCCalculatorType::Pointer myFilter =  ADCCalculatorType::New();
	myFilter->SetBValue  ( bvalue );
	myFilter->SetB0Image ( B0Image );
	myFilter->SetInput   ( InputImage );
 
	// execute le filtre
	try
	{
		myFilter->Update();
	}
	catch (itk::ExceptionObject &e)
	{
		std::cerr << e;
		return -1;
	}
 
	InputImage = myFilter->GetOutput();
 
 
 
	// ecrire resultat
	typedef itk::ImageFileWriter<ImageType> WriterType;
	WriterType::Pointer writer = WriterType::New();
	writer->SetInput ( InputImage );
	writer->SetFileName ( output );
	try
	{
		writer->Update();
	}
	catch (itk::ExceptionObject &e)
	{
		std::cerr << e;
		return -1;
	}
 
	return 0;
 
}

**Davidbrcz** · 07/01/2010, 19h43

OMG, tu ne connais ni les boucles ni les tableaux ? Car déclarer à coup de copier/coller 48 variables, il faut être courageux (et un peu fou sur les bords)