MNE-CPP  0.1.9
A Framework for Electrophysiology
ecd_set.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "ecd_set.h"
42 #include <fiff/fiff_types.h>
43 
44 //ToDo don't use access and unlink -> use QT stuff instead -> QFile
45 #if defined(WIN32) || defined(_WIN32) || defined(__WIN32) && !defined(__CYGWIN__)
46 #include <io.h>
47 #else
48 #include <unistd.h>
49 #endif
50 
51 //=============================================================================================================
52 // QT INCLUDES
53 //=============================================================================================================
54 
55 #include <QDebug>
56 #include <QFile>
57 #include <QString>
58 
59 //=============================================================================================================
60 // EIGEN INCLUDES
61 //=============================================================================================================
62 
63 //=============================================================================================================
64 // USED NAMESPACES
65 //=============================================================================================================
66 
67 using namespace INVERSELIB;
68 using namespace FIFFLIB;
69 
70 //=============================================================================================================
71 // DEFINES
72 //=============================================================================================================
73 
74 #define X 0
75 #define Y 1
76 #define Z 2
77 
78 //=============================================================================================================
79 // DEFINE STATIC METHODS
80 //=============================================================================================================
81 
82 //TODO OBSOLETE - USE ALREADY DEFINED ONE
83 static fiff_int_t swap_int (fiff_int_t source)
84 {
85  unsigned char *csource = (unsigned char *)(&source);
86  fiff_int_t result;
87  unsigned char *cresult = (unsigned char *)(&result);
88 
89  cresult[0] = csource[3];
90  cresult[1] = csource[2];
91  cresult[2] = csource[1];
92  cresult[3] = csource[0];
93  return (result);
94 }
95 
96 //=============================================================================================================
97 //TODO OBSOLETE - USE ALREADY DEFINED ONE
98 static float swap_float (float source)
99 {
100  unsigned char *csource = (unsigned char *)(&source);
101  float result;
102  unsigned char *cresult = (unsigned char *)(&result);
103 
104  cresult[0] = csource[3];
105  cresult[1] = csource[2];
106  cresult[2] = csource[1];
107  cresult[3] = csource[0];
108  return result;
109 }
110 
111 //=============================================================================================================
112 
113 namespace INVERSELIB {
114 
115 typedef struct {
116  int dipole; /* Which dipole in a multi-dipole set */
117  float begin,end; /* Fitting time range */
118  float r0[3]; /* Sphere model origin */
119  float rd[3]; /* Dipole location */
120  float Q[3]; /* Dipole amplitude */
121  float goodness; /* Goodness-of-fit */
122  int errors_computed; /* Have we computed the errors */
123  float noise_level; /* Noise level used for error computations */
124  float single_errors[5]; /* Single parameter error limits */
125  float error_matrix[5][5]; /* This fully describes the conf. ellipsoid */
126  float conf_vol; /* The xyz confidence volume */
127  float khi2; /* The khi^2 value */
128  float prob; /* Probability to exceed khi^2 by chance */
129  float noise_est; /* Total noise estimate */
131 
132 } // Namespace
133 
134 //=============================================================================================================
135 // DEFINE MEMBER METHODS
136 //=============================================================================================================
137 
139 {
140 }
141 
142 //=============================================================================================================
143 
144 ECDSet::ECDSet(const ECDSet &p_ECDSet)
145 : dataname(p_ECDSet.dataname)
146 , m_qListDips(p_ECDSet.m_qListDips)
147 {
148 }
149 
150 //=============================================================================================================
151 
153 {
154 }
155 
156 //=============================================================================================================
157 
158 void ECDSet::addEcd(const ECD& p_ecd)
159 {
160  m_qListDips.append(p_ecd);
161 }
162 
163 //=============================================================================================================
164 
165 ECDSet ECDSet::read_dipoles_dip(const QString& fileName)
166 {
167  ECDSet set;
168 
169  QFile inputFile(fileName);
170  if (inputFile.open(QIODevice::ReadOnly|QIODevice::Text))
171  {
172  QTextStream in(&inputFile);
173  while (!in.atEnd())
174  {
175  QString line = in.readLine();
176  QStringList list = line.split(QRegExp("\\s+"));
177 
178  if(list[0].contains("#") || list.size() != 11) {
179  continue;
180  }
181  else {
182  ECD one;
183  one.valid = true;
184  one.time = list[1].toFloat() / 1000.0f;
185  one.rd[X] = list[3].toFloat() / 1000.0f;
186  one.rd[Y] = list[4].toFloat() / 1000.0f;
187  one.rd[Z] = list[5].toFloat() / 1000.0f;
188  one.Q[X] = list[7].toFloat() / 1e9f;
189  one.Q[Y] = list[8].toFloat() / 1e9f;
190  one.Q[Z] = list[9].toFloat() / 1e9f;
191  one.good = list[10].toFloat() / 100.0f;
192  set << one;
193  }
194  }
195  inputFile.close();
196 
197  printf("Read %d dipoles in dip format from %s\n",set.size(),fileName.toUtf8().data());
198  }
199  else {
200  printf("Not able to read from: %s\n", fileName.toUtf8().data());
201  }
202 
203  return set;
204 }
205 
206 //=============================================================================================================
207 
208 bool ECDSet::save_dipoles_bdip(const QString& fileName)
209 /*
210  * Save dipoles in the bdip format employed by xfit
211  */
212 {
213  FILE *out = NULL;
214  bdipEcdRec one_out;
215  ECD one;
216  int k,p;
217  int nsave;
218 
219  if (fileName.isEmpty() || this->size() == 0)
220  return true;
221 
222  if ((out = fopen(fileName.toUtf8().data(),"w")) == NULL) {
223  printf(fileName.toUtf8().data());
224  return false;
225  }
226 
227  for (k = 0, nsave = 0; k < this->size(); k++) {
228  one = m_qListDips[k];
229  if (one.valid) {
230  one_out.dipole = swap_int(1);
231  one_out.begin = swap_float(one.time);
232  for (p = 0; p < 3; p++) {
233  one_out.r0[p] = swap_float(0.0);
234  one_out.rd[p] = swap_float(one.rd[p]);
235  one_out.Q[p] = swap_float(one.Q[p]);
236  }
237  one_out.goodness = swap_float(one.good);
238  one_out.errors_computed = swap_int(0);
239  one_out.khi2 = swap_float(one.khi2);
240  if (fwrite(&one_out,sizeof(bdipEcdRec),1,out) != 1) {
241  printf("Failed to write a dipole");
242  goto bad;
243  }
244  nsave++;
245  }
246  }
247  if (fclose(out) != 0) {
248  out = NULL;
249  printf(fileName.toUtf8().data());
250  return false;
251  }
252  printf("Save %d dipoles in bdip format to %s\n",nsave,fileName.toUtf8().data());
253  return true;
254 
255 bad : {
256  if (out) {
257  fclose(out);
258  unlink(fileName.toUtf8().data());
259  }
260  return false;
261  }
262 }
263 
264 //=============================================================================================================
265 
266 bool ECDSet::save_dipoles_dip(const QString& fileName) const
267 {
268  FILE *out = NULL;
269  int k,nsave;
270  ECD one;
271 
272  if (fileName.isEmpty() || this->size() == 0)
273  return true;
274  if ((out = fopen(fileName.toUtf8().data(),"w")) == NULL) {
275  printf(fileName.toUtf8().data());
276  return false;
277  }
278  fprintf(out,"# CoordinateSystem \"Head\"\n");
279  fprintf (out,"# %7s %7s %8s %8s %8s %8s %8s %8s %8s %6s\n",
280  "begin","end","X (mm)","Y (mm)","Z (mm)","Q(nAm)","Qx(nAm)","Qy(nAm)","Qz(nAm)","g/%");
281  for (k = 0, nsave = 0; k < this->size(); k++) {
282  one = this->m_qListDips[k];
283  if (one.valid) {
284  fprintf(out," %7.1f %7.1f %8.2f %8.2f %8.2f %8.3f %8.3f %8.3f %8.3f %6.1f\n",
285  1000*one.time,1000*one.time,
286  1000*one.rd[X],1000*one.rd[Y],1000*one.rd[Z],
287  1e9*one.Q.norm(),1e9*one.Q[X],1e9*one.Q[Y],1e9*one.Q[Z],100.0*one.good);
288  nsave++;
289  }
290  }
291  fprintf(out,"## Name \"%s dipoles\" Style \"Dipoles\"\n","ALL");
292  if (fclose(out) != 0) {
293  out = NULL;
294  printf(fileName.toUtf8().data());
295  goto bad;
296  }
297  printf("Save %d dipoles in dip format to %s\n",nsave,fileName.toUtf8().data());
298  return true;
299 
300 bad : {
301  if (out) {
302  fclose(out);
303  unlink(fileName.toUtf8().data());
304  }
305  return false;
306  }
307 }
308 
309 //=============================================================================================================
310 
311 const ECD& ECDSet::operator[] (int idx) const
312 {
313  if (idx>=m_qListDips.length())
314  {
315  qWarning("Warning: Required ECD doesn't exist! Returning ECD '0'.");
316  idx=0;
317  }
318  return m_qListDips[idx];
319 }
320 
321 //=============================================================================================================
322 
324 {
325  if (idx >= m_qListDips.length())
326  {
327  qWarning("Warning: Required ECD doesn't exist! Returning ECD '0'.");
328  idx = 0;
329  }
330  return m_qListDips[idx];
331 }
332 
333 //=============================================================================================================
334 
336 {
337  this->m_qListDips.append(p_ecd);
338  return *this;
339 }
INVERSELIB::ECD
Electric Current Dipole description.
Definition: ecd.h:72
INVERSELIB::ECD::valid
bool valid
Definition: ecd.h:107
INVERSELIB::ECDSet::addEcd
void addEcd(const ECD &p_ecd)
Definition: ecd_set.cpp:158
INVERSELIB::ECD::rd
Eigen::Vector3f rd
Definition: ecd.h:109
ecd_set.h
FiffDigPointSet class declaration.
INVERSELIB::ECD::khi2
float khi2
Definition: ecd.h:112
INVERSELIB::ECDSet::operator<<
ECDSet & operator<<(const ECD &p_ecd)
Definition: ecd_set.cpp:335
INVERSELIB::ECDSet::ECDSet
ECDSet()
Definition: ecd_set.cpp:138
k
int k
Definition: fiff_tag.cpp:322
INVERSELIB::ECD::good
float good
Definition: ecd.h:111
INVERSELIB::ECDSet::operator[]
const ECD & operator[](int idx) const
Definition: ecd_set.cpp:311
INVERSELIB::ECD::Q
Eigen::Vector3f Q
Definition: ecd.h:110
INVERSELIB::ECDSet::save_dipoles_bdip
bool save_dipoles_bdip(const QString &fileName)
Definition: ecd_set.cpp:208
INVERSELIB::ECDSet
Holds a set of Electric Current Dipoles.
Definition: ecd_set.h:80
INVERSELIB::ECDSet::read_dipoles_dip
static ECDSet read_dipoles_dip(const QString &fileName)
Definition: ecd_set.cpp:165
fiff_types.h
Definitions for describing the objects in a FIFF file.
INVERSELIB::ECD::time
float time
Definition: ecd.h:108
INVERSELIB::bdipEcd
Definition: ecd_set.cpp:115
INVERSELIB::ECDSet::save_dipoles_dip
bool save_dipoles_dip(const QString &fileName) const
Definition: ecd_set.cpp:266
INVERSELIB::ECDSet::~ECDSet
~ECDSet()
Definition: ecd_set.cpp:152
INVERSELIB::ECDSet::size
qint32 size() const
Definition: ecd_set.h:193