MNE-CPP  0.1.9
A Framework for Electrophysiology
fiff_info.cpp
Go to the documentation of this file.
1 //=============================================================================================================
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "fiff_info.h"
43 #include "fiff_stream.h"
44 #include "fiff_file.h"
45 
46 #include <utils/ioutils.h>
47 #include <iostream>
48 
49 //=============================================================================================================
50 // USED NAMESPACES
51 //=============================================================================================================
52 
53 using namespace FIFFLIB;
54 using namespace UTILSLIB;
55 using namespace Eigen;
56 
57 //=============================================================================================================
58 // DEFINE MEMBER METHODS
59 //=============================================================================================================
60 
62 : FiffInfoBase()//nchan(-1)
63 , sfreq(-1.0)
64 , linefreq(-1.0)
65 , highpass(-1.0)
66 , lowpass(-1.0)
67 , gantry_angle(-1)
68 , acq_pars("")
69 , acq_stim("")
70 {
71  meas_date[0] = -1;
72 }
73 
74 //=============================================================================================================
75 
76 FiffInfo::FiffInfo(const FiffInfo& p_FiffInfo)
77 : FiffInfoBase(p_FiffInfo)
78 , file_id(p_FiffInfo.file_id)
79 , sfreq(p_FiffInfo.sfreq)
80 , linefreq(p_FiffInfo.linefreq)
81 , highpass(p_FiffInfo.highpass)
82 , lowpass(p_FiffInfo.lowpass)
83 , proj_id(p_FiffInfo.proj_id)
84 , proj_name(p_FiffInfo.proj_name)
85 , xplotter_layout(p_FiffInfo.xplotter_layout)
86 , experimenter(p_FiffInfo.experimenter)
87 , description(p_FiffInfo.description)
88 , utc_offset(p_FiffInfo.utc_offset)
89 , gantry_angle(p_FiffInfo.gantry_angle)
90 , dev_ctf_t(p_FiffInfo.dev_ctf_t)
91 , dig(p_FiffInfo.dig)
92 , dig_trans(p_FiffInfo.dig_trans)
93 , projs(p_FiffInfo.projs)
94 , comps(p_FiffInfo.comps)
95 , acq_pars(p_FiffInfo.acq_pars)
96 , acq_stim(p_FiffInfo.acq_stim)
97 {
98  meas_date[0] = p_FiffInfo.meas_date[0];
99  meas_date[1] = p_FiffInfo.meas_date[1];
100 }
101 
102 //=============================================================================================================
103 
105 {
106 }
107 
108 //=============================================================================================================
109 
111 {
113  file_id = FiffId();
114  meas_date[0] = -1;
115  sfreq = -1.0;
116  linefreq = -1.0;
117  highpass = -1.0;
118  lowpass = -1.0;
119  proj_id = -1;
120  proj_name = "";
121  xplotter_layout = "";
122  experimenter = "";
123  description = "";
124  utc_offset = "";
125  gantry_angle = -1;
126  dev_ctf_t.clear();
127  dig.clear();
128  dig_trans.clear();
129  projs.clear();
130  comps.clear();
131  acq_pars = "";
132  acq_stim = "";
133 }
134 
135 //=============================================================================================================
136 
138 {
139  qint32 comp = 0;
140  qint32 first_comp = -1;
141 
142  qint32 k = 0;
143  for (k = 0; k < this->nchan; ++k)
144  {
145  if (this->chs[k].kind == FIFFV_MEG_CH)
146  {
147  comp = this->chs[k].chpos.coil_type >> 16;
148  if (first_comp < 0)
149  first_comp = comp;
150  else if (comp != first_comp)
151  printf("Compensation is not set equally on all MEG channels");
152  }
153  }
154  return comp;
155 }
156 
157 //=============================================================================================================
158 
159 bool FiffInfo::make_compensator(fiff_int_t from, fiff_int_t to, FiffCtfComp& ctf_comp, bool exclude_comp_chs) const
160 {
161  MatrixXd C1, C2, comp_tmp;
162 
163 // if(ctf_comp.data)
164 // delete ctf_comp.data;
165  ctf_comp.data->clear();
166 
167  if (from == to)
168  {
169  ctf_comp.data->data = MatrixXd::Identity(this->nchan, this->nchan);
170  return false;
171  }
172 
173  if (from == 0)
174  C1 = MatrixXd::Zero(this->nchan,this->nchan);
175  else
176  {
177  if (!this->make_compensator(from, C1))
178  {
179  printf("Cannot create compensator C1\n");
180  printf("Desired compensation matrix (kind = %d) not found\n",from);
181  return false;
182  }
183  }
184 
185  if (to == 0)
186  C2 = MatrixXd::Zero(this->nchan,this->nchan);
187  else
188  {
189  if (!this->make_compensator(to, C2))
190  {
191  printf("Cannot create compensator C2\n");
192  printf("Desired compensation matrix (kind = %d) not found\n",to);
193  return false;
194  }
195  }
196  //
197  // s_orig = s_from + C1*s_from = (I + C1)*s_from
198  // s_to = s_orig - C2*s_orig = (I - C2)*s_orig
199  // s_to = (I - C2)*(I + C1)*s_from = (I + C1 - C2 - C2*C1)*s_from
200  //
201  comp_tmp = MatrixXd::Identity(this->nchan,this->nchan) + C1 - C2 - C2*C1;
202 
203  qint32 k;
204  if (exclude_comp_chs)
205  {
206  VectorXi pick = MatrixXi::Zero(1,this->nchan);
207  qint32 npick = 0;
208  for (k = 0; k < this->nchan; ++k)
209  {
210  if (this->chs[k].kind != FIFFV_REF_MEG_CH)
211  {
212  pick(npick) = k;
213  ++npick;
214  }
215  }
216  if (npick == 0)
217  {
218  printf("Nothing remains after excluding the compensation channels\n");
219  return false;
220  }
221 
222  ctf_comp.data->data.resize(npick,this->nchan);
223  for (k = 0; k < npick; ++k)
224  ctf_comp.data->data.row(k) = comp_tmp.block(pick(k), 0, 1, this->nchan);
225  } else {
226  ctf_comp.data->data = comp_tmp;
227  }
228 
229  return true;
230 }
231 
232 //=============================================================================================================
233 
234 bool FiffInfo::make_compensator(fiff_int_t kind, MatrixXd& this_comp) const//private method
235 {
236  FiffNamedMatrix::SDPtr this_data;
237  MatrixXd presel, postsel;
238  qint32 k, col, c, ch=0, row, row_ch=0, channelAvailable;
239 
240  for (k = 0; k < this->comps.size(); ++k)
241  {
242  if (this->comps[k].kind == kind)
243  {
244  this_data = this->comps[k].data;
245 
246  //
247  // Create the preselector
248  //
249  presel = MatrixXd::Zero(this_data->ncol,this->nchan);
250 
251  for(col = 0; col < this_data->ncol; ++col)
252  {
253  channelAvailable = 0;
254  for (c = 0; c < this->ch_names.size(); ++c)
255  {
256  if (QString::compare(this_data->col_names.at(col),this->ch_names.at(c)) == 0)
257  {
258  ++channelAvailable;
259  ch = c;
260  }
261  }
262  if (channelAvailable == 0)
263  {
264  printf("Channel %s is not available in data\n",this_data->col_names.at(col).toUtf8().constData());
265  return false;
266  }
267  else if (channelAvailable > 1)
268  {
269  printf("Ambiguous channel %s",this_data->col_names.at(col).toUtf8().constData());
270  return false;
271  }
272  presel(col,ch) = 1.0;
273  }
274  //
275  // Create the postselector
276  //
277  postsel = MatrixXd::Zero(this->nchan,this_data->nrow);
278 
279  for (c = 0; c < this->nchan; ++c)
280  {
281  channelAvailable = 0;
282  for (row = 0; row < this_data->row_names.size(); ++row)
283  {
284  if (QString::compare(this->ch_names.at(c),this_data->row_names.at(row)) == 0)
285  {
286  ++channelAvailable;
287  row_ch = row;
288  }
289  }
290  if (channelAvailable > 1)
291  {
292  printf("Ambiguous channel %s", this->ch_names.at(c).toUtf8().constData());
293  return false;
294  }
295  else if (channelAvailable == 1)
296  {
297  postsel(c,row_ch) = 1.0;
298  }
299  }
300  this_comp = postsel*this_data->data*presel;
301  return true;
302  }
303  }
304  this_comp = defaultMatrixXd;
305  return false;
306 }
307 
308 //=============================================================================================================
309 
310 FiffInfo FiffInfo::pick_info(const RowVectorXi &sel) const
311 {
312  FiffInfo res = *this;//new FiffInfo(this);
313  if (sel.size() == 0)
314  return res;
315 
316  //ToDo when pointer List do delation
317  res.chs.clear();
318  res.ch_names.clear();
319 
320  qint32 idx;
321  for(qint32 i = 0; i < sel.size(); ++i)
322  {
323  idx = sel[i];
324  res.chs.append(this->chs[idx]);
325  res.ch_names.append(this->ch_names[idx]);
326  }
327  res.nchan = static_cast<int>(sel.size());
328 
329  return res;
330 }
331 
332 //=============================================================================================================
333 
334 QList<FiffChInfo> FiffInfo::set_current_comp(QList<FiffChInfo>& listFiffChInfo, fiff_int_t value)
335 {
336  QList<FiffChInfo> newList;
337  qint32 k;
338  fiff_int_t coil_type;
339 
340  for(k = 0; k < listFiffChInfo.size(); ++k)
341  newList.append(listFiffChInfo[k]);
342 
343  qint32 lower_half = 65535;// hex2dec('FFFF');
344  for (k = 0; k < listFiffChInfo.size(); ++k)
345  {
346  if (listFiffChInfo[k].kind == FIFFV_MEG_CH)
347  {
348  coil_type = listFiffChInfo[k].chpos.coil_type & lower_half;
349  newList[k].chpos.coil_type = (coil_type | (value << 16));
350  }
351  }
352  return newList;
353 }
354 
355 //=============================================================================================================
356 
357 void FiffInfo::writeToStream(FiffStream* p_pStream) const
358 {
359  //
360  // We will always write floats
361  //
362  fiff_int_t data_type = 4;
363  qint32 k;
364  QList<FiffChInfo> chs;
365 
366  for(k = 0; k < this->nchan; ++k)
367  chs << this->chs[k];
368 
369  fiff_int_t nchan = chs.size();
370 
371  //
372  // write the essentials
373  //
374  p_pStream->start_block(FIFFB_MEAS);//4
375  p_pStream->write_id(FIFF_BLOCK_ID);//5
376  if(this->meas_id.version != -1)
377  {
378  p_pStream->write_id(FIFF_PARENT_BLOCK_ID,this->meas_id);//6
379  }
380  //
381  // Measurement info
382  //
383  p_pStream->start_block(FIFFB_MEAS_INFO);//7
384 
385  //
386  // Blocks from the original -> skip this
387  //
388 // QList<fiff_int_t> blocks;
389 // blocks << FIFFB_SUBJECT << FIFFB_HPI_MEAS << FIFFB_HPI_RESULT << FIFFB_ISOTRAK << FIFFB_PROCESSING_HISTORY;
390  bool have_hpi_result = false;
391  bool have_isotrak = false;
392  //
393  // megacq parameters
394  //
395  if (!this->acq_pars.isEmpty() || !this->acq_stim.isEmpty())
396  {
397  p_pStream->start_block(FIFFB_DACQ_PARS);
398  if (!this->acq_pars.isEmpty())
399  p_pStream->write_string(FIFF_DACQ_PARS, this->acq_pars);
400 
401  if (!this->acq_stim.isEmpty())
402  p_pStream->write_string(FIFF_DACQ_STIM, this->acq_stim);
403 
404  p_pStream->end_block(FIFFB_DACQ_PARS);
405  }
406  //
407  // Coordinate transformations if the HPI result block was not there
408  //
409  if (!have_hpi_result)
410  {
411  if (!this->dev_head_t.isEmpty())
412  p_pStream->write_coord_trans(this->dev_head_t);
413 
414  if (!this->ctf_head_t.isEmpty())
415  p_pStream->write_coord_trans(this->ctf_head_t);
416  }
417  //
418  // Polhemus data
419  //
420  if (this->dig.size() > 0 && !have_isotrak)
421  {
422  p_pStream->start_block(FIFFB_ISOTRAK);
423  for (qint32 k = 0; k < this->dig.size(); ++k)
424  p_pStream->write_dig_point(this->dig[k]);
425 
426  p_pStream->end_block(FIFFB_ISOTRAK);
427  }
428  //
429  // Projectors
430  //
431  p_pStream->write_proj(this->projs);
432  //
433  // CTF compensation info
434  //
435  p_pStream->write_ctf_comp(this->comps);
436  //
437  // Bad channels
438  //
439  if (this->bads.size() > 0)
440  {
441  p_pStream->start_block(FIFFB_MNE_BAD_CHANNELS);
442  p_pStream->write_name_list(FIFF_MNE_CH_NAME_LIST,this->bads);
443  p_pStream->end_block(FIFFB_MNE_BAD_CHANNELS);
444  }
445  //
446  // General
447  //
448  p_pStream->write_float(FIFF_SFREQ,&this->sfreq);
449  p_pStream->write_float(FIFF_LINE_FREQ,&this->linefreq);
450  p_pStream->write_float(FIFF_HIGHPASS,&this->highpass);
451  p_pStream->write_float(FIFF_LOWPASS,&this->lowpass);
452  p_pStream->write_string(FIFF_EXPERIMENTER,this->experimenter);
453  p_pStream->write_string(FIFF_DESCRIPTION,this->description);
454  p_pStream->write_string(FIFF_UNIT_C,this->utc_offset);
455  p_pStream->write_string(FIFF_PROJ_NAME,this->proj_name);
456  p_pStream->write_int(FIFF_PROJ_ID,&this->proj_id);
457  p_pStream->write_int(FIFF_GANTRY_ANGLE,&this->gantry_angle);
458  p_pStream->write_int(FIFF_NCHAN,&nchan);
459  p_pStream->write_int(FIFF_DATA_PACK,&data_type);
460  if (this->meas_date[0] != -1)
461  p_pStream->write_int(FIFF_MEAS_DATE,this->meas_date, 2);
462  //
463  // Channel info
464  //
465  MatrixXd cals(1,nchan);
466 
467  for(k = 0; k < nchan; ++k)
468  {
469  //
470  // Scan numbers may have been messed up
471  //
472  chs[k].scanNo = k+1;//+1 because
473 // chs[k].range = 1.0f;//Why? -> cause its already calibrated through reading
474  cals(0,k) = static_cast<double>(chs[k].cal); //ToDo whats going on with cals?
475  p_pStream->write_ch_info(chs[k]);
476  }
477  //
478  //
479  p_pStream->end_block(FIFFB_MEAS_INFO);
480 }
481 
482 //=============================================================================================================
483 
484 void FiffInfo::print() const
485 {
486  std::cout << "Sample frequency: " << sfreq << "\n";
487  std::cout << "LineFreq: " << linefreq << " | Highpass: " << highpass << " | Lowpass: " << lowpass << "\n";
488  std::cout << "Number of digitizer points: " << dig.size() << "\n";
489  for (auto& point : dig){
490  if (point.kind == FIFFV_POINT_HPI){
491  std::cout << "HPI Point " << point.ident << " - " << point.r[0] << ", " << point.r[1] << ", " << point.r[2] << "\n";
492  }
493  }
494 }
#define FIFF_NCHAN
Definition: fiff_file.h:453
QList< FiffChInfo > chs
fiff_long_t write_ctf_comp(const QList< FiffCtfComp > &comps)
#define FIFF_DESCRIPTION
Definition: fiff_file.h:486
FiffInfo class declaration.
#define FIFF_DATA_PACK
Definition: fiff_file.h:455
QString acq_stim
Definition: fiff_info.h:271
FiffCoordTrans ctf_head_t
FiffStream class declaration.
fiff_int_t gantry_angle
Definition: fiff_info.h:264
FiffNamedMatrix::SDPtr data
CTF software compensation data.
Definition: fiff_ctf_comp.h:73
IOUtils class declaration.
fiff_long_t write_name_list(fiff_int_t kind, const QStringList &data)
#define FIFF_PROJ_ID
Definition: fiff_file.h:577
fiff_long_t write_dig_point(const FiffDigPoint &dig)
Header file describing the numerical values used in fif files.
#define FIFF_EXPERIMENTER
Definition: fiff_file.h:465
#define FIFF_HIGHPASS
Definition: fiff_file.h:476
QString proj_name
Definition: fiff_info.h:259
fiff_long_t write_proj(const QList< FiffProj > &projs)
fiff_long_t write_coord_trans(const FiffCoordTrans &trans)
#define FIFFV_REF_MEG_CH
FiffCoordTrans dev_head_t
fiff_long_t end_block(fiff_int_t kind, fiff_int_t next=FIFFV_NEXT_SEQ)
bool make_compensator(fiff_int_t from, fiff_int_t to, FiffCtfComp &ctf_comp, bool exclude_comp_chs=false) const
Definition: fiff_info.cpp:159
QList< FiffCtfComp > comps
Definition: fiff_info.h:269
qint32 get_current_comp()
Definition: fiff_info.cpp:137
QSharedDataPointer< FiffNamedMatrix > SDPtr
fiff_long_t start_block(fiff_int_t kind)
fiff_long_t write_float(fiff_int_t kind, const float *data, fiff_int_t nel=1)
FiffCoordTrans dev_ctf_t
Definition: fiff_info.h:265
#define FIFF_LINE_FREQ
Definition: fiff_file.h:489
fiff_long_t write_id(fiff_int_t kind, const FiffId &id=FiffId::getDefault())
fiff_long_t write_int(fiff_int_t kind, const fiff_int_t *data, fiff_int_t nel=1, fiff_int_t next=FIFFV_NEXT_SEQ)
QString xplotter_layout
Definition: fiff_info.h:260
void set_current_comp(fiff_int_t value)
Definition: fiff_info.h:292
fiff_int_t meas_date[2]
Definition: fiff_info.h:253
FiffInfo pick_info(const Eigen::RowVectorXi &sel=defaultVectorXi) const
Definition: fiff_info.cpp:310
#define FIFF_MEAS_DATE
Definition: fiff_file.h:457
QString utc_offset
Definition: fiff_info.h:263
FiffCoordTrans dig_trans
Definition: fiff_info.h:267
QList< FiffProj > projs
Definition: fiff_info.h:268
FIFF measurement file information.
Definition: fiff_info.h:84
void print() const
Definition: fiff_info.cpp:484
#define FIFF_GANTRY_ANGLE
Definition: fiff_file.h:557
#define FIFF_PROJ_NAME
Definition: fiff_file.h:578
QString acq_pars
Definition: fiff_info.h:270
#define FIFF_SFREQ
Definition: fiff_file.h:454
fiff_int_t version
Definition: fiff_id.h:152
QList< FiffDigPoint > dig
Definition: fiff_info.h:266
Universially unique identifier.
Definition: fiff_id.h:68
#define FIFF_LOWPASS
Definition: fiff_file.h:472
fiff_long_t write_ch_info(const FiffChInfo &ch)
QString description
Definition: fiff_info.h:262
FIFF File I/O routines.
Definition: fiff_stream.h:104
void writeToStream(FiffStream *p_pStream) const
Definition: fiff_info.cpp:357
QString experimenter
Definition: fiff_info.h:261
fiff_long_t write_string(fiff_int_t kind, const QString &data)
light measurement info