MNE-CPP  0.1.9
A Framework for Electrophysiology
dipole_fit.cpp
1 
2 
3 #include "dipole_fit.h"
4 #include "../c/mne_meas_data_set.h"
5 #include "guess_data.h"
6 
7 #include <string.h>
8 #include <QScopedPointer>
9 
10 using namespace INVERSELIB;
11 using namespace MNELIB;
12 using namespace FWDLIB;
13 
14 #if defined(_WIN32) || defined(_WIN64)
15 #define snprintf _snprintf
16 #define vsnprintf _vsnprintf
17 #define strcasecmp _stricmp
18 #define strncasecmp _strnicmp
19 #endif
20 
21 //=============================================================================================================
22 // STATIC DEFINITIONS
23 //=============================================================================================================
24 
25 #ifndef PROGRAM_VERSION
26 #define PROGRAM_VERSION "1.00"
27 #endif
28 
29 #ifndef FAIL
30 #define FAIL -1
31 #endif
32 
33 #ifndef OK
34 #define OK 0
35 #endif
36 
37 #ifndef TRUE
38 #define TRUE 1
39 #endif
40 
41 #ifndef FALSE
42 #define FALSE 0
43 #endif
44 
45 #define BIG_TIME 1e6
46 
47 #define SEG_LEN 10.0
48 
49 #define EPS_VALUES 0.05
50 
51 //=============================================================================================================
52 // STATIC DEFINITIONS ToDo make members
53 //=============================================================================================================
54 
55 #define MALLOC(x,t) (t *)malloc((x)*sizeof(t))
56 #define REALLOC(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
57 #define FREE(x) if ((char *)(x) != NULL) free((char *)(x))
58 
59 #define ALLOC_CMATRIX(x,y) mne_cmatrix((x),(y))
60 
61 #define FREE_CMATRIX(m) mne_free_cmatrix((m))
62 
63 static void matrix_error(int kind, int nr, int nc)
64 
65 {
66  if (kind == 1)
67  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
68  else if (kind == 2)
69  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
70  else
71  printf("Allocation error for a %d x %d matrix\n",nr,nc);
72  if (sizeof(void *) == 4) {
73  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
74  printf("Please consider moving to a 64-bit platform.");
75  }
76  printf("Cannot continue. Sorry.\n");
77  exit(1);
78 }
79 
80 float **mne_cmatrix (int nr,int nc)
81 
82 {
83  int i;
84  float **m;
85  float *whole;
86 
87  m = MALLOC(nr,float *);
88  if (!m) matrix_error(1,nr,nc);
89  whole = MALLOC(nr*nc,float);
90  if (!whole) matrix_error(2,nr,nc);
91 
92  for(i=0;i<nr;i++)
93  m[i] = whole + i*nc;
94  return m;
95 }
96 
97 void mne_free_cmatrix (float **m)
98 {
99  if (m) {
100  FREE(*m);
101  FREE(m);
102  }
103 }
104 
105 //char *mne_compose_mne_name(const char *path, const char *filename)
107 // * Compose a filename under the "$MNE_ROOT" directory
108 // */
109 //{
110 // char *res;
111 // char *mne_root;
112 
113 // if (filename == NULL) {
114 // qCritical("No file name specified to mne_compose_mne_name");
115 // return NULL;
116 // }
117 // mne_root = getenv(MNE_ENV_ROOT);
118 // if (mne_root == NULL || strlen(mne_root) == 0) {
119 // qCritical("Environment variable MNE_ROOT not set");
120 // return NULL;
121 // }
122 // if (path == NULL || strlen(path) == 0) {
123 // res = MALLOC(strlen(mne_root)+strlen(filename)+2,char);
124 // strcpy(res,mne_root);
125 // strcat(res,"/");
126 // strcat(res,filename);
127 // }
128 // else {
129 // res = MALLOC(strlen(mne_root)+strlen(filename)+strlen(path)+3,char);
130 // strcpy(res,mne_root);
131 // strcat(res,"/");
132 // strcat(res,path);
133 // strcat(res,"/");
134 // strcat(res,filename);
135 // }
136 // return res;
137 //}
138 
139 //============================= misc_util.c =============================
140 
141 //============================= mne_named_matrix.c =============================
142 
143 void mne_free_name_list(char **list, int nlist)
144 /*
145  * Free a name list array
146  */
147 {
148  int k;
149  if (list == NULL || nlist == 0)
150  return;
151  for (k = 0; k < nlist; k++) {
152 #ifdef FOO
153  printf("%d %s\n",k,list[k]);
154 #endif
155  FREE(list[k]);
156  }
157  FREE(list);
158  return;
159 }
160 
161 //============================= mne_ch_selections.c =============================
162 
163 /*
164  * Mandatory allocation functions
165  */
166 static mneChSelection new_ch_selection()
167 
168 {
169  mneChSelection newsel = MALLOC(1,mneChSelectionRec);
170 
171  newsel->name.clear();
172  newsel->chdef.clear();
173  newsel->chspick.clear();
174  newsel->chspick_nospace.clear();
175  newsel->pick = NULL;
176  newsel->pick_deriv = NULL;
177  newsel->ch_kind = NULL;
178  newsel->ndef = 0;
179  newsel->nchan = 0;
180  newsel->kind = MNE_CH_SELECTION_UNKNOWN;
181  return newsel;
182 }
183 
184 mneChSelection mne_ch_selection_these(const QString& selname, const QStringList& names, int nch)
185 /*
186  * Give an explicit list of interesting channels
187  */
188 {
189  int c;
190  mneChSelection sel;
191 
192  sel = new_ch_selection();
193  sel->name = selname;
194 // sel->chdef;
195  sel->ndef = nch;
196  sel->kind = MNE_CH_SELECTION_USER;
197 
198  for (c = 0; c < nch; c++)
199  sel->chdef.append(names[c]);
200 
201  return sel;
202 }
203 
204 static void omit_spaces(QStringList names, int nnames)
205 
206 {
207  char *c,*cc;
208  for (int k = 0; k < names.size(); k++) {
209  while( names[k].startsWith(" ") ) {
210  names[k].remove(0,1);
211  }
212  }
213  return;
214 }
215 
216 int mne_ch_selection_assign_chs(mneChSelection sel,
217  MneRawData* data)
218 /*
219  * Make the channel picking real easy
220  */
221 {
222  int c,rc,d;
223  MneRawInfo* info;
224  int nch;
225  QString dash;
226 
227  if (!sel || !data)
228  return 0;
229 
230  info = data->info;
231  sel->chspick.clear();
232  sel->chspick_nospace.clear();
233  /*
234  * Expansion of possible regular expressions must be added eventually
235  */
236  sel->chspick = sel->chdef;
237  sel->chspick_nospace = sel->chdef;
238  omit_spaces(sel->chspick_nospace,sel->ndef);
239  sel->nchan = sel->ndef;
240 
241  sel->pick = REALLOC(sel->pick,sel->nchan,int); /* Just in case */
242  sel->pick_deriv = REALLOC(sel->pick_deriv,sel->nchan,int);
243  sel->ch_kind = REALLOC(sel->ch_kind,sel->nchan,int);
244 
245  for (c = 0; c < sel->nchan; c++) {
246  sel->pick[c] = -1;
247  sel->pick_deriv[c] = -1;
248  sel->ch_kind[c] = -1;
249  for (rc = 0; rc < info->nchan; rc++) {
250  if (QString::compare(sel->chspick[c],info->chInfo[rc].ch_name,Qt::CaseInsensitive) == 0 ||
251  QString::compare(sel->chspick_nospace[c],info->chInfo[rc].ch_name,Qt::CaseInsensitive) == 0) {
252  sel->pick[c] = rc;
253  sel->ch_kind[c] = info->chInfo[rc].kind;
254  break;
255  }
256  }
257  }
258  /*
259  * Maybe the derivations will help
260  */
261  sel->nderiv = 0;
262  if (data->deriv_matched) {
263  QStringList deriv_names = data->deriv_matched->deriv_data->rowlist;
264  int nderiv = data->deriv_matched->deriv_data->nrow;
265 
266  for (c = 0; c < sel->nchan; c++) {
267  if (sel->pick[c] == -1) {
268  for (d = 0; d < nderiv; d++) {
269  if (QString::compare(sel->chspick[c],deriv_names[d],Qt::CaseInsensitive) == 0 &&
270  data->deriv_matched->valid && data->deriv_matched->valid[d]) {
271  sel->pick_deriv[c] = d;
272  sel->ch_kind[c] = data->deriv_matched->chs[d].kind;
273  sel->nderiv++;
274  break;
275  }
276  }
277  }
278  }
279  }
280  /*
281  * Try simple channels again without the part after dashes
282  */
283  for (c = 0; c < sel->nchan; c++) {
284  if (sel->pick[c] == -1 && sel->pick_deriv[c] == -1) {
285  for (rc = 0; rc < info->nchan; rc++) {
286  dash = QString(info->chInfo[rc].ch_name).mid(QString(info->chInfo[rc].ch_name).indexOf("-")+1);// strchr(info->chInfo[rc].ch_name,'-');
287  if (!dash.isNull()) {
288 // *dash = '\0';
289  if (QString::compare(sel->chspick[c],info->chInfo[rc].ch_name,Qt::CaseInsensitive) == 0 ||
290  QString::compare(sel->chspick_nospace[c],info->chInfo[rc].ch_name,Qt::CaseInsensitive) == 0) {
291 // *dash = '-';
292  sel->pick[c] = rc;
293  sel->ch_kind[c] = info->chInfo[rc].kind;
294  break;
295  }
296 // *dash = '-';
297  }
298  }
299  }
300  }
301  for (c = 0, nch = 0; c < sel->nchan; c++) {
302  if (sel->pick[c] >= 0)
303  nch++;
304  }
305  if (sel->nderiv > 0)
306  printf("Selection %c%s%c has %d matched derived channels.\n",'"',sel->name.toUtf8().constData(),'"',sel->nderiv);
307  return nch;
308 }
309 
310 //============================= mne_get_values.c =============================
311 
312 int mne_get_values_from_data (float time, /* Interesting time point */
313  float integ, /* Time integration */
314  float **data, /* The data values (time by time) */
315  int nsamp, /* How many time points? */
316  int nch, /* How many channels */
317  float tmin, /* Time of first sample */
318  float sfreq, /* Sampling frequency */
319  int use_abs, /* Use absolute values */
320  float *value) /* The picked values */
321 /*
322  * Pick a signal value using linear interpolation
323  */
324 {
325  int n1,n2,k;
326  float s1,s2;
327  float f1,f2;
328  float sum;
329  float width;
330  int ch;
331 
332  for (ch = 0; ch < nch; ch++) {
333  /*
334  * Find out the correct samples
335  */
336  if (std::fabs(sfreq*integ) < EPS_VALUES) { /* This is the single-sample case */
337  s1 = sfreq*(time - tmin);
338  n1 = floor(s1);
339  f1 = 1.0 + n1 - s1;
340  if (n1 < 0 || n1 > nsamp-1) {
341  printf("Sample value out of range %d (0..%d)",n1,nsamp-1);
342  return(-1);
343  }
344  /*
345  * Avoid rounding error
346  */
347  if (n1 == nsamp-1) {
348  if (std::fabs(f1-1.0) < 1e-3)
349  f1 = 1.0;
350  }
351  if (f1 < 1.0 && n1 > nsamp-2) {
352  printf("Sample value out of range %d (0..%d) %.4f",n1,nsamp-1,f1);
353  return(-1);
354  }
355  if (f1 < 1.0) {
356  if (use_abs)
357  sum = f1*std::fabs(data[n1][ch]) + (1.0-f1)*std::fabs(data[n1+1][ch]);
358  else
359  sum = f1*data[n1][ch] + (1.0-f1)*data[n1+1][ch];
360  }
361  else {
362  if (use_abs)
363  sum = std::fabs(data[n1][ch]);
364  else
365  sum = data[n1][ch];
366  }
367  }
368  else { /* Multiple samples */
369  s1 = sfreq*(time - 0.5*integ - tmin);
370  s2 = sfreq*(time + 0.5*integ - tmin);
371  n1 = ceil(s1); n2 = floor(s2);
372  if (n2 < n1) { /* We are within one sample interval */
373  n1 = floor(s1);
374  if (n1 < 0 || n1 > nsamp-2)
375  return (-1);
376  f1 = s1 - n1;
377  f2 = s2 - n1;
378  if (use_abs)
379  sum = 0.5*((f1+f2)*std::fabs(data[n1+1][ch]) + (2.0-f1-f2)*std::fabs(data[n1][ch]));
380  else
381  sum = 0.5*((f1+f2)*data[n1+1][ch] + (2.0-f1-f2)*data[n1][ch]);
382  }
383  else {
384  f1 = n1 - s1;
385  f2 = s2 - n2;
386  if (n1 < 0 || n1 > nsamp-1) {
387  printf("Sample value out of range %d (0..%d)",n1,nsamp-1);
388  return(-1);
389  }
390  if (n2 < 0 || n2 > nsamp-1) {
391  printf("Sample value out of range %d (0..%d)",n2,nsamp-1);
392  return(-1);
393  }
394  if (f1 != 0.0 && n1 < 1)
395  return(-1);
396  if (f2 != 0.0 && n2 > nsamp-2)
397  return(-1);
398  sum = 0.0;
399  width = 0.0;
400  if (n2 > n1) { /* Do the whole intervals */
401  if (use_abs) {
402  sum = 0.5 * std::fabs(data[n1][ch]);
403  for (k = n1+1; k < n2; k++)
404  sum = sum + std::fabs(data[k][ch]);
405  sum = sum + 0.5 * std::fabs(data[n2][ch]);
406  }
407  else {
408  sum = 0.5*data[n1][ch];
409  for (k = n1+1; k < n2; k++)
410  sum = sum + data[k][ch];
411  sum = sum + 0.5*data[n2][ch];
412  }
413  width = n2 - n1;
414  }
415  /*
416  * Add tails
417  */
418  if (use_abs) {
419  if (f1 != 0.0)
420  sum = sum + 0.5 * f1 * (f1 * std::fabs(data[n1-1][ch]) + (2.0 - f1) * std::fabs(data[n1][ch]));
421  if (f2 != 0.0)
422  sum = sum + 0.5 * f2 * (f2 * std::fabs(data[n2+1][ch]) + (2.0 - f2) * std::fabs(data[n2][ch]));
423  }
424  else {
425  if (f1 != 0.0)
426  sum = sum + 0.5*f1*(f1*data[n1-1][ch] + (2.0-f1)*data[n1][ch]);
427  if (f2 != 0.0)
428  sum = sum + 0.5*f2*(f2*data[n2+1][ch] + (2.0-f2)*data[n2][ch]);
429  }
430  width = width + f1 + f2;
431  sum = sum/width;
432  }
433  }
434  value[ch] = sum;
435  }
436  return (0);
437 }
438 
439 int mne_get_values_from_data_ch (float time, /* Interesting time point */
440  float integ, /* Time integration */
441  float **data, /* The data values (channel by channel) */
442  int nsamp, /* How many time points? */
443  int nch, /* How many channels */
444  float tmin, /* Time of first sample */
445  float sfreq, /* Sampling frequency */
446  int use_abs, /* Use absolute values */
447  float *value) /* The picked values */
448 /*
449  * Pick a signal value using linear interpolation
450  */
451 {
452  int n1,n2,k;
453  float s1,s2;
454  float f1,f2;
455  float sum;
456  float width;
457  int ch;
458 
459  for (ch = 0; ch < nch; ch++) {
460  /*
461  * Find out the correct samples
462  */
463  if (std::fabs(sfreq * integ) < EPS_VALUES) { /* This is the single-sample case */
464  s1 = sfreq*(time - tmin);
465  n1 = floor(s1);
466  f1 = 1.0 + n1 - s1;
467  if (n1 < 0 || n1 > nsamp-1)
468  return(-1);
469  if (f1 < 1.0 && n1 > nsamp-2)
470  return(-1);
471  if (f1 < 1.0) {
472  if (use_abs)
473  sum = f1 * std::fabs(data[ch][n1]) + (1.0 - f1) * std::fabs(data[ch][n1+1]);
474  else
475  sum = f1*data[ch][n1] + (1.0-f1)*data[ch][n1+1];
476  }
477  else {
478  if (use_abs)
479  sum = std::fabs(data[ch][n1]);
480  else
481  sum = data[ch][n1];
482  }
483  }
484  else { /* Multiple samples */
485  s1 = sfreq*(time - 0.5*integ - tmin);
486  s2 = sfreq*(time + 0.5*integ - tmin);
487  n1 = ceil(s1); n2 = floor(s2);
488  if (n2 < n1) { /* We are within one sample interval */
489  n1 = floor(s1);
490  if (n1 < 0 || n1 > nsamp-2)
491  return (-1);
492  f1 = s1 - n1;
493  f2 = s2 - n1;
494  if (use_abs)
495  sum = 0.5*((f1+f2)*std::fabs(data[ch][n1+1]) + (2.0-f1-f2)*std::fabs(data[ch][n1]));
496  else
497  sum = 0.5*((f1+f2)*data[ch][n1+1] + (2.0-f1-f2)*data[ch][n1]);
498  }
499  else {
500  f1 = n1 - s1;
501  f2 = s2 - n2;
502  if (n1 < 0 || n1 > nsamp-1 || n2 < 0 || n2 > nsamp-1)
503  return(-1);
504  if (f1 != 0.0 && n1 < 1)
505  return(-1);
506  if (f2 != 0.0 && n2 > nsamp-2)
507  return(-1);
508  sum = 0.0;
509  width = 0.0;
510  if (n2 > n1) { /* Do the whole intervals */
511  if (use_abs) {
512  sum = 0.5 * std::fabs(data[ch][n1]);
513  for (k = n1+1; k < n2; k++)
514  sum = sum + std::fabs(data[ch][k]);
515  sum = sum + 0.5 * std::fabs(data[ch][n2]);
516  }
517  else {
518  sum = 0.5*data[ch][n1];
519  for (k = n1+1; k < n2; k++)
520  sum = sum + data[ch][k];
521  sum = sum + 0.5*data[ch][n2];
522  }
523  width = n2 - n1;
524  }
525  /*
526  * Add tails
527  */
528  if (use_abs) {
529  if (f1 != 0.0)
530  sum = sum + 0.5 * f1 * (f1 * std::fabs(data[ch][n1-1]) + (2.0 - f1) * std::fabs(data[ch][n1]));
531  if (f2 != 0.0)
532  sum = sum + 0.5 * f2 * (f2 * std::fabs(data[ch][n2+1]) + (2.0 - f2) * std::fabs(data[ch][n2]));
533  }
534  else {
535  if (f1 != 0.0)
536  sum = sum + 0.5*f1*(f1*data[ch][n1-1]+ (2.0-f1)*data[ch][n1]);
537  if (f2 != 0.0)
538  sum = sum + 0.5*f2*(f2*data[ch][n2+1] + (2.0-f2)*data[ch][n2]);
539  }
540  width = width + f1 + f2;
541  sum = sum/width;
542  }
543  }
544  value[ch] = sum;
545  }
546  return (0);
547 }
548 
549 //=============================================================================================================
550 // DEFINE MEMBER METHODS
551 //=============================================================================================================
552 
554 : settings(p_settings)
555 {
556 }
557 
558 //=============================================================================================================
559 //todo split in initFit where the settings are handed over and the actual fit
560 ECDSet DipoleFit::calculateFit() const
561 {
562  QScopedPointer<GuessData>guess (Q_NULLPTR);
563  ECDSet set;
564  FwdEegSphereModel* eeg_model = NULL;
565  DipoleFitData* fit_data = NULL;
566  MneMeasData* data = NULL;
567  MneRawData* raw = NULL;
568  mneChSelection sel = NULL;
569 
570  printf("---- Setting up...\n\n");
571  if (settings->include_eeg) {
572  if ((eeg_model = FwdEegSphereModel::setup_eeg_sphere_model(settings->eeg_model_file,settings->eeg_model_name,settings->eeg_sphere_rad)) == NULL)
573  goto out;
574  }
575 
576  if ((fit_data = DipoleFitData::setup_dipole_fit_data( settings->mriname,
577  settings->measname,
578  settings->bemname.isEmpty() ? NULL : settings->bemname.toUtf8().data(),
579  &settings->r0,
580  eeg_model,
581  settings->accurate,
582  settings->badname,
583  settings->noisename,
584  settings->grad_std,
585  settings->mag_std,
586  settings->eeg_std,
587  settings->mag_reg,
588  settings->grad_reg,
589  settings->eeg_reg,
590  settings->diagnoise,
591  settings->projnames,
592  settings->include_meg,
593  settings->include_eeg)) == NULL)
594  goto out;
595 
596  fit_data->fit_mag_dipoles = settings->fit_mag_dipoles;
597  if (settings->is_raw) {
598  int c;
599  float t1,t2;
600 
601  printf("\n---- Opening a raw data file...\n\n");
602  if ((raw = MneRawData::mne_raw_open_file(settings->measname.isEmpty() ? NULL : settings->measname.toUtf8().data(),TRUE,FALSE,&(settings->filter))) == NULL)
603  goto out;
604  /*
605  * A channel selection is needed to access the data
606  */
607  sel = mne_ch_selection_these("fit",fit_data->ch_names,fit_data->nmeg+fit_data->neeg);
608  mne_ch_selection_assign_chs(sel,raw);
609  for (c = 0; c < sel->nchan; c++)
610  if (sel->pick[c] < 0) {
611  qCritical ("All desired channels were not available");
612  goto out;
613  }
614  printf("\tChannel selection created.\n");
615  /*
616  * Let's be a little generous here
617  */
618  t1 = raw->first_samp/raw->info->sfreq;
619  t2 = (raw->first_samp+raw->nsamp-1)/raw->info->sfreq;
620  if (settings->tmin < t1 + settings->integ)
621  settings->tmin = t1 + settings->integ;
622  if (settings->tmax > t2 - settings->integ)
623  settings->tmax = t2 - settings->integ;
624  if (settings->tstep < 0)
625  settings->tstep = 1.0/raw->info->sfreq;
626 
627  printf("\tOpened raw data file %s : %d MEG and %d EEG \n",
628  settings->measname.toUtf8().data(),fit_data->nmeg,fit_data->neeg);
629  }
630  else {
631  printf("\n---- Reading data...\n\n");
632  if ((data = MneMeasData::mne_read_meas_data(settings->measname,
633  settings->setno,
634  NULL,
635  NULL,
636  fit_data->ch_names,
637  fit_data->nmeg+fit_data->neeg)) == NULL)
638  goto out;
639  if (settings->do_baseline)
640  data->adjust_baselines(settings->bmin,settings->bmax);
641  else
642  printf("\tNo baseline setting in effect.\n");
643  if (settings->tmin < data->current->tmin + settings->integ/2.0)
644  settings->tmin = data->current->tmin + settings->integ/2.0;
645  if (settings->tmax > data->current->tmin + (data->current->np-1)*data->current->tstep - settings->integ/2.0)
646  settings->tmax = data->current->tmin + (data->current->np-1)*data->current->tstep - settings->integ/2.0;
647  if (settings->tstep < 0)
648  settings->tstep = data->current->tstep;
649 
650  printf("\tRead data set %d from %s : %d MEG and %d EEG \n",
651  settings->setno,settings->measname.toUtf8().data(),fit_data->nmeg,fit_data->neeg);
652  if (!settings->noisename.isEmpty()) {
653  printf("\nScaling the noise covariance...\n");
654  if (DipoleFitData::scale_noise_cov(fit_data,data->current->nave) == FAIL)
655  goto out;
656  }
657  }
658 
659  /*
660  * Proceed to computing the fits
661  */
662  printf("\n---- Computing the forward solution for the guesses...\n\n");
663  guess.reset(new GuessData( settings->guessname,
664  settings->guess_surfname,
665  settings->guess_mindist, settings->guess_exclude, settings->guess_grid, fit_data));
666  if (guess.isNull())
667  goto out;
668 
669  fprintf (stderr,"\n---- Fitting : %7.1f ... %7.1f ms (step: %6.1f ms integ: %6.1f ms)\n\n",
670  1000*settings->tmin,1000*settings->tmax,1000*settings->tstep,1000*settings->integ);
671 
672  if (raw) {
673  if (fit_dipoles_raw(settings->measname,raw,sel,fit_data,guess.take(),settings->tmin,settings->tmax,settings->tstep,settings->integ,settings->verbose) == FAIL)
674  goto out;
675  }
676  else {
677  if (fit_dipoles(settings->measname,data,fit_data,guess.take(),settings->tmin,settings->tmax,settings->tstep,settings->integ,settings->verbose,set) == FAIL)
678  goto out;
679  }
680  printf("%d dipoles fitted\n",set.size());
681 
682 out : {
683  return set;
684  }
685 }
686 
687 //=============================================================================================================
688 
689 int DipoleFit::fit_dipoles( const QString& dataname, MneMeasData* data, DipoleFitData* fit, GuessData* guess, float tmin, float tmax, float tstep, float integ, int verbose, ECDSet& p_set)
690 {
691  float *one = MALLOC(data->nchan,float);
692  float time;
693  ECDSet set;
694  ECD dip;
695  int s;
696  int report_interval = 10;
697 
698  set.dataname = dataname;
699 
700  printf("Fitting...%c",verbose ? '\n' : '\0');
701  for (s = 0, time = tmin; time < tmax; s++, time = tmin + s*tstep) {
702  /*
703  * Pick the data point
704  */
705  if (mne_get_values_from_data(time,integ,data->current->data,data->current->np,data->nchan,data->current->tmin,
706  1.0/data->current->tstep,FALSE,one) == FAIL) {
707  printf("Cannot pick time: %7.1f ms\n",1000*time);
708  continue;
709  }
710 
711  if (!DipoleFitData::fit_one(fit,guess,time,one,verbose,dip))
712  printf("t = %7.1f ms : %s\n",1000*time,"error (tbd: catch)");
713  else {
714  set.addEcd(dip);
715  if (verbose)
716  dip.print(stdout);
717  else {
718  if (set.size() % report_interval == 0)
719  printf("%d..",set.size());
720  }
721  }
722  }
723  if (!verbose)
724  printf("[done]\n");
725  FREE(one);
726  p_set = set;
727  return OK;
728 }
729 
730 //=============================================================================================================
731 
732 int DipoleFit::fit_dipoles_raw(const QString& dataname, MneRawData* raw, mneChSelection sel, DipoleFitData* fit, GuessData* guess, float tmin, float tmax, float tstep, float integ, int verbose, ECDSet& p_set)
733 {
734  float *one = MALLOC(sel->nchan,float);
735  float sfreq = raw->info->sfreq;
736  float myinteg = integ > 0.0 ? 2*integ : 0.1;
737  int overlap = ceil(myinteg*sfreq);
738  int length = SEG_LEN*sfreq;
739  int step = length - overlap;
740  int stepo = step + overlap/2;
741  int start = raw->first_samp;
742  int s,picks;
743  float time,stime;
744  float **data = ALLOC_CMATRIX(sel->nchan,length);
745  ECD dip;
746  ECDSet set;
747  int report_interval = 10;
748 
749  set.dataname = dataname;
750 
751  /*
752  * Load the initial data segment
753  */
754  stime = start/sfreq;
755  if (MneRawData::mne_raw_pick_data_filt(raw,sel,start,length,data) == FAIL)
756  goto bad;
757  printf("Fitting...%c",verbose ? '\n' : '\0');
758  for (s = 0, time = tmin; time < tmax; s++, time = tmin + s*tstep) {
759  picks = time*sfreq - start;
760  if (picks > stepo) { /* Need a new data segment? */
761  start = start + step;
762  if (MneRawData::mne_raw_pick_data_filt(raw,sel,start,length,data) == FAIL)
763  goto bad;
764  picks = time*sfreq - start;
765  stime = start/sfreq;
766  }
767  /*
768  * Get the values
769  */
770  if (mne_get_values_from_data_ch (time,integ,data,length,sel->nchan,stime,sfreq,FALSE,one) == FAIL) {
771  printf("Cannot pick time: %8.3f s\n",time);
772  continue;
773  }
774  /*
775  * Fit
776  */
777  if (!DipoleFitData::fit_one(fit,guess,time,one,verbose,dip))
778  qWarning() << "Error";
779  else {
780  set.addEcd(dip);
781  if (verbose)
782  dip.print(stdout);
783  else {
784  if (set.size() % report_interval == 0)
785  printf("%d..",set.size());
786  }
787  }
788  }
789  if (!verbose)
790  printf("[done]\n");
791  FREE_CMATRIX(data);
792  FREE(one);
793  p_set = set;
794  return OK;
795 
796 bad : {
797  FREE_CMATRIX(data);
798  FREE(one);
799  return FAIL;
800  }
801 }
802 
803 //=============================================================================================================
804 
805 int DipoleFit::fit_dipoles_raw(const QString& dataname, MneRawData* raw, mneChSelection sel, DipoleFitData* fit, GuessData* guess, float tmin, float tmax, float tstep, float integ, int verbose)
806 {
807  ECDSet set;
808  return fit_dipoles_raw(dataname, raw, sel, fit, guess, tmin, tmax, tstep, integ, verbose, set);
809 }
GuessData class declaration.
static bool fit_one(DipoleFitData *fit, GuessData *guess, float time, float *B, int verbose, ECD &res)
Dipole Fit setting implementation.
A comprehensive raw data structure.
Definition: mne_raw_data.h:86
static DipoleFitData * setup_dipole_fit_data(const QString &mriname, const QString &measname, const QString &bemname, Eigen::Vector3f *r0, FWDLIB::FwdEegSphereModel *eeg_model, int accurate_coils, const QString &badname, const QString &noisename, float grad_std, float mag_std, float eeg_std, float mag_reg, float grad_reg, float eeg_reg, int diagnoise, const QList< QString > &projnames, int include_meg, int include_eeg)
GuessData description.
Definition: guess_data.h:78
QList< FIFFLIB::FiffChInfo > chInfo
Definition: mne_raw_info.h:130
static int fit_dipoles_raw(const QString &dataname, MNELIB::MneRawData *raw, MNELIB::mneChSelection sel, DipoleFitData *fit, GuessData *guess, float tmin, float tmax, float tstep, float integ, int verbose, ECDSet &p_set)
Definition: dipole_fit.cpp:732
Electric Current Dipole description.
Definition: ecd.h:72
Dipole Fit Data implementation.
QString dataname
Definition: ecd_set.h:176
DipoleFit(DipoleFitSettings *p_settings)
Definition: dipole_fit.cpp:553
Holds a set of Electric Current Dipoles.
Definition: ecd_set.h:80
Electric Current Dipole description.
static int fit_dipoles(const QString &dataname, MneMeasData *data, DipoleFitData *fit, GuessData *guess, float tmin, float tmax, float tstep, float integ, int verbose, ECDSet &p_set)
Definition: dipole_fit.cpp:689
Dipole Fit class declaration.
easurement data representation in MNE calculations
Definition: mne_meas_data.h:92
Information about raw data in fiff file.
Definition: mne_raw_info.h:80
void adjust_baselines(float bmin, float bmax)
void print(FILE *f) const
Definition: ecd.cpp:88