MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_epoch_data_list.cpp
Go to the documentation of this file.
1 //=============================================================================================================
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "mne_epoch_data_list.h"
43 
44 #include <utils/mnemath.h>
45 
46 //=============================================================================================================
47 // QT INCLUDES
48 //=============================================================================================================
49 
50 #include <QPointer>
51 #include <QtConcurrent>
52 #include <QDebug>
53 
54 //=============================================================================================================
55 // USED NAMESPACES
56 //=============================================================================================================
57 
58 using namespace FIFFLIB;
59 using namespace MNELIB;
60 using namespace UTILSLIB;
61 using namespace Eigen;
62 
63 //=============================================================================================================
64 // DEFINE MEMBER METHODS
65 //=============================================================================================================
66 
67 MNEEpochDataList::MNEEpochDataList()
68 {
69 }
70 
71 //=============================================================================================================
72 
73 MNEEpochDataList::~MNEEpochDataList()
74 {
75 // MNEEpochDataList::iterator i;
76 // for( i = this->begin(); i!=this->end(); ++i) {
77 // if (*i)
78 // delete (*i);
79 // }
80 }
81 
82 //=============================================================================================================
83 
84 MNEEpochDataList MNEEpochDataList::readEpochs(const FiffRawData& raw,
85  const MatrixXi& events,
86  float tmin,
87  float tmax,
88  qint32 event,
89  const QMap<QString,double>& mapReject,
90  const QStringList& lExcludeChs,
91  const RowVectorXi& picks)
92 {
93  MNEEpochDataList data;
94 
95  // Select the desired events
96  qint32 count = 0;
97  qint32 p;
98  MatrixXi selected = MatrixXi::Zero(1, events.rows());
99  for (p = 0; p < events.rows(); ++p)
100  {
101  if (events(p,1) == 0 && events(p,2) == event)
102  {
103  selected(0,count) = p;
104  ++count;
105  }
106  }
107  selected.conservativeResize(1, count);
108  if (count > 0) {
109  qInfo("[MNEEpochDataList::readEpochs] %d matching events found",count);
110  } else {
111  qWarning("[MNEEpochDataList::readEpochs] No desired events found.");
112  return MNEEpochDataList();
113  }
114 
115  // If picks are empty, pick all
116  RowVectorXi picksNew = picks;
117  if(picks.cols() <= 0) {
118  picksNew.resize(raw.info.chs.size());
119  for(int i = 0; i < raw.info.chs.size(); ++i) {
120  picksNew(i) = i;
121  }
122  }
123 
124  fiff_int_t event_samp, from, to;
125  fiff_int_t dropCount = 0;
126  MatrixXd timesDummy;
127  MatrixXd times;
128 
129  QScopedPointer<MNEEpochData> epoch(Q_NULLPTR);
130 
131  for (p = 0; p < count; ++p) {
132  // Read a data segment
133  event_samp = events(selected(p),0);
134  from = event_samp + tmin*raw.info.sfreq;
135  to = event_samp + floor(tmax*raw.info.sfreq + 0.5);
136 
137  epoch.reset(new MNEEpochData());
138 
139  if(raw.read_raw_segment(epoch->epoch, timesDummy, from, to, picksNew)) {
140  if (p == 0) {
141  times.resize(1, to-from+1);
142  for (qint32 i = 0; i < times.cols(); ++i)
143  times(0, i) = ((float)(from-event_samp+i)) / raw.info.sfreq;
144  }
145 
146  epoch->event = event;
147  epoch->tmin = tmin;
148  epoch->tmax = tmax;
149 
150  epoch->bReject = checkForArtifact(epoch->epoch,
151  raw.info,
152  mapReject,
153  lExcludeChs);
154 
155  if (epoch->bReject) {
156  dropCount++;
157  }
158 
159  //Check if data block has the same size as the previous one
160  if(!data.isEmpty()) {
161  if(epoch->epoch.size() == data.last()->epoch.size()) {
162  data.append(MNEEpochData::SPtr(epoch.take()));//List takes ownwership of the pointer - no delete need
163  }
164  } else {
165  data.append(MNEEpochData::SPtr(epoch.take()));//List takes ownwership of the pointer - no delete need
166  }
167  } else {
168  qWarning("[MNEEpochDataList::readEpochs] Can't read the event data segments.");
169  }
170  }
171 
172  qInfo().noquote() << "[MNEEpochDataList::readEpochs] Read a total of"<< data.size() <<"epochs of type" << event << "and marked"<< dropCount <<"for rejection.";
173 
174  return data;
175 }
176 
177 //=============================================================================================================
178 
179 FiffEvoked MNEEpochDataList::average(const FiffInfo& info,
180  fiff_int_t first,
181  fiff_int_t last,
182  VectorXi sel,
183  bool proj)
184 {
185  FiffEvoked p_evoked;
186 
187  qInfo("[MNEEpochDataList::average] Calculate evoked. ");
188 
189  MatrixXd matAverage;
190 
191  if(this->size() > 0) {
192  matAverage = MatrixXd::Zero(this->at(0)->epoch.rows(), this->at(0)->epoch.cols());
193  } else {
195  return p_evoked;
196  }
197 
198  if(sel.size() > 0) {
199  p_evoked.nave = sel.size();
200 
201  for(qint32 i = 0; i < sel.size(); ++i) {
202  matAverage.array() += this->at(sel(i))->epoch.array();
203  }
204  } else {
205  p_evoked.nave = this->size();
206 
207  for(qint32 i = 0; i < this->size(); ++i) {
208  matAverage.array() += this->at(i)->epoch.array();
209  }
210  }
211  matAverage.array() /= p_evoked.nave;
212 
213  qInfo("[MNEEpochDataList::average] %d averages used [done]", p_evoked.nave);
214 
215  p_evoked.setInfo(info, proj);
216 
218 
219  p_evoked.first = first;
220  p_evoked.last = last;
221 
222  p_evoked.times = RowVectorXf::LinSpaced(this->first()->epoch.cols(), this->first()->tmin, this->first()->tmax);
223 
224  p_evoked.times[static_cast<int>(this->first()->tmin * -1 * info.sfreq)] = 0;
225 
226  p_evoked.comment = QString::number(this->at(0)->event);
227 
228  if(p_evoked.proj.rows() > 0) {
229  matAverage = p_evoked.proj * matAverage;
230  qInfo("[MNEEpochDataList::average] SSP projectors applied to the evoked data");
231  }
232 
233  p_evoked.data = matAverage;
234 
235  return p_evoked;
236 }
237 
238 //=============================================================================================================
239 
240 void MNEEpochDataList::applyBaselineCorrection(const QPair<float, float> &baseline)
241 {
242  // Run baseline correction
243  QMutableListIterator<MNEEpochData::SPtr> i(*this);
244  while (i.hasNext()) {
245  i.next()->applyBaselineCorrection(baseline);
246  }
247 }
248 
249 //=============================================================================================================
250 
251 void MNEEpochDataList::dropRejected()
252 {
253  QMutableListIterator<MNEEpochData::SPtr> i(*this);
254  while (i.hasNext()) {
255  if (i.next()->bReject) {
256  i.remove();
257  }
258  }
259 }
260 
261 //=============================================================================================================
262 
263 void MNEEpochDataList::pick_channels(const RowVectorXi& sel)
264 {
265  QMutableListIterator<MNEEpochData::SPtr> i(*this);
266  while (i.hasNext()) {
267  i.next()->pick_channels(sel);
268  }
269 }
270 
271 //=============================================================================================================
272 
273 bool MNEEpochDataList::checkForArtifact(const MatrixXd& data,
274  const FiffInfo& pFiffInfo,
275  const QMap<QString,double>& mapReject,
276  const QStringList& lExcludeChs)
277 {
278  //qDebug() << "MNEEpochDataList::checkForArtifact - Doing artifact reduction for" << mapReject;
279 
280  bool bReject = false;
281 
282  //Prepare concurrent data handling
283  QList<ArtifactRejectionData> lchData;
284  QList<int> lChTypes;
285 
286  if(mapReject.contains("grad") ||
287  mapReject.contains("mag") ) {
288  lChTypes << FIFFV_MEG_CH;
289  }
290 
291  if(mapReject.contains("eeg")) {
292  lChTypes << FIFFV_EEG_CH;
293  }
294 
295  if(mapReject.contains("eog")) {
296  lChTypes << FIFFV_EOG_CH;
297  }
298 
299  if(lChTypes.isEmpty()) {
300  return bReject;
301  }
302 
303  for(int i = 0; i < pFiffInfo.chs.size(); ++i) {
304  if(lChTypes.contains(pFiffInfo.chs.at(i).kind)
305  && !lExcludeChs.contains(pFiffInfo.chs.at(i).ch_name)
306  && !pFiffInfo.bads.contains(pFiffInfo.chs.at(i).ch_name)
307  && pFiffInfo.chs.at(i).chpos.coil_type != FIFFV_COIL_BABY_REF_MAG
308  && pFiffInfo.chs.at(i).chpos.coil_type != FIFFV_COIL_BABY_REF_MAG2) {
309  ArtifactRejectionData tempData;
310  tempData.data = data.row(i);
311 
312  switch (pFiffInfo.chs.at(i).kind) {
313  case FIFFV_MEG_CH:
314  if(pFiffInfo.chs.at(i).unit == FIFF_UNIT_T) {
315  tempData.dThreshold = mapReject["mag"];
316  } else if(pFiffInfo.chs.at(i).unit == FIFF_UNIT_T_M) {
317  tempData.dThreshold = mapReject["grad"];
318  }
319  break;
320 
321  case FIFFV_EEG_CH:
322  tempData.dThreshold = mapReject["eeg"];
323  break;
324 
325  case FIFFV_EOG_CH:
326  tempData.dThreshold = mapReject["eog"];
327  break;
328  }
329 
330  tempData.sChName = pFiffInfo.chs.at(i).ch_name;
331  lchData.append(tempData);
332  }
333  }
334 
335  if(lchData.isEmpty()) {
336  qWarning() << "[MNEEpochDataList::checkForArtifact] No channels found to scan for artifacts. Do not reject. Returning.";
337 
338  return bReject;
339  }
340 
341  //qDebug() << "MNEEpochDataList::checkForArtifact - lchData.size()" << lchData.size();
342 
343  //Start the concurrent processing
344  QFuture<void> future = QtConcurrent::map(lchData, checkChThreshold);
345  future.waitForFinished();
346 
347  for(int i = 0; i < lchData.size(); ++i) {
348  if(lchData.at(i).bRejected) {
349  bReject = true;
350  qInfo().noquote() << "[MNEEpochDataList::checkForArtifact] Reject trial because of channel"<<lchData.at(i).sChName;
351  break;
352  }
353  }
354 
355  return bReject;
356 }
357 
358 //=============================================================================================================
359 
360 void MNEEpochDataList::checkChThreshold(ArtifactRejectionData& inputData)
361 {
362  RowVectorXd temp = inputData.data;
363 
364  // Remove offset
365  //temp = temp.array() - temp(0);
366 
367  double min = temp.minCoeff();
368  double max = temp.maxCoeff();
369 
370  // Peak to Peak
371  double pp = max - min;
372 
373  if(std::fabs(pp) > inputData.dThreshold) {
374  inputData.bRejected = true;
375  } else {
376  inputData.bRejected = false;
377  }
378 
379 // qDebug() << "MNEEpochDataList::checkChThreshold - min" << min;
380 // qDebug() << "MNEEpochDataList::checkChThreshold - max" << max;
381 // qDebug() << "MNEEpochDataList::checkChThreshold - pp" << pp;
382 // qDebug() << "MNEEpochDataList::checkChThreshold - inputData.dThreshold" << inputData.dThreshold;
383 
384 // //If absolute vaue of min or max if bigger than threshold -> reject
385 // if((std::fabs(min) > inputData.dThreshold) || (std::fabs(max) > inputData.dThreshold)) {
386 // inputData.bRejected = true;
387 // } else {
388 // inputData.bRejected = false;
389 // }
390 }
Eigen::RowVectorXf times
Definition: fiff_evoked.h:230
FIFF measurement file information.
Definition: fiff_info.h:84
#define FIFFV_ASPECT_STD_ERR
Definition: fiff_file.h:439
FIFF raw measurement data.
Definition: fiff_raw_data.h:78
fiff_int_t aspect_kind
Definition: fiff_evoked.h:226
bool read_raw_segment(Eigen::MatrixXd &data, Eigen::MatrixXd &times, fiff_int_t from=-1, fiff_int_t to=-1, const Eigen::RowVectorXi &sel=defaultRowVectorXi, bool do_debug=false) const
evoked data
Definition: fiff_evoked.h:77
#define FIFFV_COIL_BABY_REF_MAG2
MNEMath class declaration.
MNEEpochDataList class declaration.
#define FIFFV_ASPECT_AVERAGE
Definition: fiff_file.h:438
#define FIFFV_COIL_BABY_REF_MAG
QSharedPointer< MNEEpochData > SPtr
void setInfo(const FiffInfo &p_info, bool proj=true)
Eigen::MatrixXd proj
Definition: fiff_evoked.h:232
Eigen::MatrixXd data
Definition: fiff_evoked.h:231
QList< FiffChInfo > chs