MNE-CPP  0.1.9
A Framework for Electrophysiology
warp.cpp
Go to the documentation of this file.
1 //=============================================================================================================
35 //=============================================================================================================
36 // INCLUDES
37 //=============================================================================================================
38 
39 #include "warp.h"
40 
41 #include <iostream>
42 #include <fstream>
43 
44 //=============================================================================================================
45 // EIGEN INCLUDES
46 //=============================================================================================================
47 
48 #include <Eigen/LU>
49 
50 //=============================================================================================================
51 // QT INCLUDES
52 //=============================================================================================================
53 
54 #include <QDebug>
55 #include <QFile>
56 #include <QList>
57 
58 //=============================================================================================================
59 // USED NAMESPACES
60 //=============================================================================================================
61 
62 using namespace UTILSLIB;
63 using namespace Eigen;
64 
65 //=============================================================================================================
66 // DEFINE MEMBER METHODS
67 //=============================================================================================================
68 
69 MatrixXf Warp::calculate(const MatrixXf &sLm, const MatrixXf &dLm, const MatrixXf &sVert)
70 {
71  MatrixXf warpWeight, polWeight;
72  calcWeighting(sLm, dLm, warpWeight, polWeight);
73  MatrixXf wVert = warpVertices(sVert, sLm, warpWeight, polWeight);
74  return wVert;
75 }
76 
77 //=============================================================================================================
78 
79 void Warp::calculate(const MatrixXf & sLm, const MatrixXf &dLm, QList<MatrixXf> & vertList)
80 {
81  MatrixXf warpWeight, polWeight;
82  calcWeighting(sLm, dLm, warpWeight, polWeight);
83 
84  for (int i=0; i<vertList.size(); i++)
85  {
86  vertList.replace(i,warpVertices(vertList.at(i), sLm, warpWeight, polWeight));
87  }
88  return;
89 }
90 
91 //=============================================================================================================
92 
93 bool Warp::calcWeighting(const MatrixXf &sLm, const MatrixXf &dLm, MatrixXf& warpWeight, MatrixXf& polWeight)
94 {
95  MatrixXf K = MatrixXf::Zero(sLm.rows(),sLm.rows()); //K(i,j)=||sLm(i)-sLm(j)||
96  for (int i=0; i<sLm.rows(); i++)
97  K.col(i)=((sLm.rowwise()-sLm.row(i)).rowwise().norm());
98 
99 // std::cout << "Here is the matrix K:" << std::endl << K << std::endl;
100 
101  MatrixXf P (sLm.rows(),4); //P=[ones,sLm]
102  P << MatrixXf::Ones(sLm.rows(),1),sLm;
103 // std::cout << "Here is the matrix P:" << std::endl << P << std::endl;
104 
105  MatrixXf L ((sLm.rows()+4),(sLm.rows()+4)); //L=Full Matrix of the linear eq.
106  L << K,P,
107  P.transpose(),MatrixXf::Zero(4,4);
108 // std::cout << "Here is the matrix L:" << std::endl << L << std::endl;
109 
110  MatrixXf Y ((dLm.rows()+4),3); //Y=[dLm,Zero]
111  Y << dLm,
112  MatrixXf::Zero(4,3);
113 // std::cout << "Here is the matrix Y:" << std::endl << Y << std::endl;
114 
115  //
116  // calculate the weighting matrix (Y=L*W)
117  //
118  MatrixXf W ((dLm.rows()+4),3); //W=[warpWeight,polWeight]
119  Eigen::FullPivLU <MatrixXf> Lu(L); //LU decomposition is one method to solve lin. eq.
120  W=Lu.solve(Y);
121 // std::cout << "Here is the matrix W:" << std::endl << W << std::endl;
122 
123  warpWeight = W.topRows(sLm.rows());
124  polWeight = W.bottomRows(4);
125 
126  return true;
127 }
128 
129 //=============================================================================================================
130 
131 MatrixXf Warp::warpVertices(const MatrixXf &sVert, const MatrixXf & sLm, const MatrixXf& warpWeight, const MatrixXf& polWeight)
132 {
133  MatrixXf wVert = sVert * polWeight.bottomRows(3); //Pol. Warp
134  wVert.rowwise() += polWeight.row(0); //Translation
135 
136  //
137  // TPS Warp
138  //
139  MatrixXf K = MatrixXf::Zero(sVert.rows(),sLm.rows()); //K(i,j)=||sLm(i)-sLm(j)||
140  for (int i=0; i<sVert.rows(); i++)
141  K.row(i)=((sLm.rowwise()-sVert.row(i)).rowwise().norm().transpose());
142 // std::cout << "Here is the matrix K:" << std::endl << K << std::endl;
143 
144  wVert += K*warpWeight;
145 // std::cout << "Here is the matrix wVert:" << std::endl << wVert << std::endl;
146  return wVert;
147 }
148 
149 //=============================================================================================================
150 
151 MatrixXf Warp::readsLm(const QString &electrodeFileName)
152 {
153  MatrixXf electrodes;
154  QFile file(electrodeFileName);
155 
156  if(!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
157  qDebug()<<"Error opening file";
158 // return false;
159  }
160 
161  //Start reading from file
162  double numberElectrodes;
163  QTextStream in(&file);
164  int i=0;
165 
166  while(!in.atEnd())
167  {
168  QString line = in.readLine();
169  QStringList fields = line.split(QRegExp("\\s+"));
170 
171  //Delete last element if it is a blank character
172  if(fields.at(fields.size()-1) == "")
173  fields.removeLast();
174 
175  //Read number of electrodes
176  if(i == 0){
177  numberElectrodes = fields.at(fields.size()-1).toDouble();
178  electrodes = MatrixXf::Zero(numberElectrodes, 3);
179  }
180  //Read actual electrode positions
181  else{
182  Vector3f x;
183  x << fields.at(fields.size()-3).toFloat(),fields.at(fields.size()-2).toFloat(),fields.at(fields.size()-1).toFloat();
184  electrodes.row(i-1)=x.transpose();
185  }
186  i++;
187  }
188  return electrodes;
189 }
190 
191 //=============================================================================================================
192 
193 MatrixXf Warp::readsLm(const std::string &electrodeFileName)
194 {
195  MatrixXf electrodes;
196  std::ifstream inFile(electrodeFileName);
197 
198  if(!inFile.is_open()) {
199  qDebug()<<"Error opening file";
200  //Why are we not returning?
201 // return false;
202  }
203 
204  //Start reading from file
205  double numberElectrodes;
206  int i = 0;
207 
208  std::string line;
209  while(std::getline(inFile, line)){
210  std::vector<std::string> fields;
211  std::stringstream stream{line};
212  std::string element;
213 
214  stream >> std::ws;
215  while(stream >> element){
216  fields.push_back(std::move(element));
217  stream >> std::ws;
218  }
219 
220  //Read number of electrodes
221  if(i == 0){
222  numberElectrodes = std::stod(fields.at(fields.size()-1));
223  electrodes = MatrixXf::Zero(numberElectrodes, 3);
224  }
225 
226  //Read actual electrode positions
227  else{
228  Vector3f x;
229  x << std::stof(fields.at(fields.size()-3)), std::stof(fields.at(fields.size()-2)), std::stof(fields.at(fields.size()-1));
230  electrodes.row(i-1)=x.transpose();
231  }
232  i++;
233  }
234 
235  return electrodes;
236 }
Eigen::MatrixXf readsLm(const QString &electrodeFileName)
Definition: warp.cpp:151
Eigen::MatrixXf calculate(const Eigen::MatrixXf &sLm, const Eigen::MatrixXf &dLm, const Eigen::MatrixXf &sVert)
Warp class declaration.