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