MNE-CPP  beta 1.0
layoutmaker.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //*************************************************************************************************************
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "layoutmaker.h"
43 
44 
45 //*************************************************************************************************************
46 //=============================================================================================================
47 // Qt INCLUDES
48 //=============================================================================================================
49 
50 #include <QVector>
51 
52 
53 //*************************************************************************************************************
54 //=============================================================================================================
55 // USED NAMESPACES
56 //=============================================================================================================
57 
58 using namespace UTILSLIB;
59 
60 
61 //*************************************************************************************************************
62 //=============================================================================================================
63 // DEFINE MEMBER METHODS
64 //=============================================================================================================
65 
67 {
68 }
69 
70 
71 //*************************************************************************************************************
72 
73 bool LayoutMaker::makeLayout(const QList<QVector<double> > &inputPoints,
74  QList<QVector<double> > &outputPoints,
75  const QStringList &names,
76  QFile &outFile,
77  bool do_fit,
78  float prad,
79  float w,
80  float h,
81  bool writeFile)
82 {
83  /*
84  * Automatically make a layout according to the
85  * channel locations in inputPoints
86  */
87  int res = FAIL;
88  VectorXf r0(3);
89  VectorXf rr(3);
90  float rad,th,phi;
91 
92  float xmin,xmax,ymin,ymax;
93  int k;
94  int nchan = inputPoints.size();
95 
96  MatrixXf rrs(nchan,3);
97  VectorXf xx(nchan);
98  VectorXf yy(nchan);
99 
100  if (nchan <= 0) {
101  std::cout<<"No input points to lay out.";
102  return false;
103  }
104 
105  //Fill matrix with 3D points
106  for(k = 0; k<nchan; k++) {
107  rrs(k,0) = inputPoints.at(k)[0]; //x
108  rrs(k,1) = inputPoints.at(k)[1]; //y
109  rrs(k,2) = inputPoints.at(k)[2]; //z
110  }
111 
112  std::cout<<"Channels found for layout: "<<nchan<<std::endl;
113 
114  //Fit to sphere if wanted by the user
115  if (!do_fit)
116  std::cout<<"Using default origin:"<<r0[0]<<r0[1]<<r0[2]<<std::endl;
117  else {
118  if(fit_sphere_to_points(rrs,nchan,0.05,r0,rad) == FAIL) {
119  std::cout<<"Using default origin:"<<r0[0]<<r0[1]<<r0[2]<<std::endl;
120  }
121  else{
122  std::cout<<"best fitting sphere:"<<std::endl;
123  std::cout<<"torigin: "<<r0[0]<<r0[1]<<r0[2]<<rad<<std::endl<<"tradius: "<<rad<<std::endl;
124  }
125  }
126 
127  /*
128  * Do the azimuthal equidistant projection
129  */
130  for (k = 0; k < nchan; k++) {
131  rr = r0 - static_cast<VectorXf>(rrs.row(k));
132  sphere_coord(rr[0],rr[1],rr[2],&rad,&th,&phi);
133  xx[k] = prad*(2.0*th/M_PI)*cos(phi);
134  yy[k] = prad*(2.0*th/M_PI)*sin(phi);
135  }
136 
137  /*
138  * Find suitable range of viewports
139  */
140  xmin = xmax = xx[0];
141  ymin = ymax = yy[0];
142 
143  for(k = 1; k < nchan; k++) {
144  if (xx[k] > xmax)
145  xmax = xx[k];
146  else if (xx[k] < xmin)
147  xmin = xx[k];
148  if (yy[k] > ymax)
149  ymax = yy[k];
150  else if (yy[k] < ymin)
151  ymin = yy[k];
152  }
153 
154  if(xmin == xmax || ymin == ymax) {
155  std::cout<<"Cannot make a layout. All positions are identical"<<std::endl;
156  return res;
157  }
158 
159  xmax = xmax + 0.6*w;
160  xmin = xmin - 0.6*w;
161 
162  ymax = ymax + 0.6*h;
163  ymin = ymin - 0.6*h;
164 
165  /*
166  * Compose the viewports
167  */
168  QVector<double> point;
169  QTextStream out;
170 
171  if(writeFile) {
172  if (!outFile.open(QIODevice::WriteOnly)) {
173  std::cout<<"could not open output file";
174  qDebug()<<"could not open output file";
175  return FAIL;
176  }
177 
178  out.setDevice(&outFile);
179  }
180 
181  for(k = 0; k < nchan; k++) {
182  point.clear();
183  point.append(xx[k]-0.5*w);
184  point.append(-(yy[k]-0.5*h)); //mirror y axis
185  outputPoints.append(point);
186 
187  if(writeFile)
188  out << k+1 << " " << point[0] << " " << point[1] << " " << w << " " << h << " " << names.at(k)<<endl;
189  }
190 
191  if(writeFile) {
192  std::cout<<"success while wrtiting to output file";
193 
194  outFile.close();
195  }
196 
197  res = OK;
198 
199  return res;
200 }
201 
202 
203 //*************************************************************************************************************
204 
205 void LayoutMaker::sphere_coord (float x,
206  float y,
207  float z,
208  float *r,
209  float *theta,
210  float *phi)
211 {
212  /* Rectangular to spherical coordinates */
213  float rxy = sqrt(x*x+y*y);
214  if (rxy < EPS) { /* Let's hope this is reasonable */
215  *r = z;
216  *theta = 0.0;
217  *phi = 0.0;
218  }
219  else {
220  *r = sqrt(x*x+y*y+z*z);
221  *theta = acos(z/(*r));
222  *phi = atan2 (y,x);
223  if (*phi < 0.0)
224  *phi = *phi + 2.0*M_PI;
225  }
226 }
227 
228 //*************************************************************************************************************
229 
230 
231 int LayoutMaker::report_func(int loop,
232  const VectorXf &fitpar,
233  int npar,
234  double fval)
235 {
236  /*
237  * Report periodically
238  */
239  VectorXf r0 = fitpar;
240 
241  std::cout<<"loop: "<<loop<<"r0: "<<1000*r0[0]<<1000*r0[1]<<1000*r0[2]<<"fval: "<<fval<<std::endl;
242 
243  return OK;
244 }
245 
246 
247 //*************************************************************************************************************
248 
249 float LayoutMaker::fit_eval(const VectorXf &fitpar,
250  int npar,
251  void *user_data)
252 {
253  /*
254  * Calculate the cost function value
255  * Optimize for the radius inside here
256  */
257  fitUser user = (fitUser)user_data;
258  VectorXf r0 = fitpar;
259  VectorXf diff(3);
260  int k;
261  float sum,sum2,one,F;
262 
263  for (k = 0, sum = sum2 = 0.0; k < user->np; k++) {
264  diff = r0 - static_cast<VectorXf>(user->rr.row(k));
265  one = sqrt(pow(diff(0),2) + pow(diff(1),2) + pow(diff(2),2));
266  sum += one;
267  sum2 += one*one;
268  }
269  F = sum2 - sum*sum/user->np;
270 
271  if(user->report)
272  std::cout<<"r0: "<<1000*r0[0]<<1000*r0[1]<<1000*r0[2]<<"R: "<<1000*sum/user->np<<"fval: "<<F<<std::endl;
273 
274  return F;
275 }
276 
277 
278 //*************************************************************************************************************
279 
280 float LayoutMaker::opt_rad(VectorXf &r0,fitUser user)
281 {
282  float sum, one;
283  VectorXf diff(3);
284  int k;
285 
286  for (k = 0, sum = 0.0; k < user->np; k++) {
287  diff = r0 - static_cast<VectorXf>(user->rr.row(k));
288  one = sqrt(pow(diff(0),2) + pow(diff(1),2) + pow(diff(2),2));
289  sum += one;
290  }
291 
292  return sum/user->np;
293 }
294 
295 
296 //*************************************************************************************************************
297 
298 void LayoutMaker::calculate_cm_ave_dist(MatrixXf &rr,
299  int np,
300  VectorXf &cm,
301  float &avep)
302 {
303  int k,q;
304  float ave;
305  VectorXf diff(3);
306 
307  for (q = 0; q < 3; q++)
308  cm[q] = 0.0;
309 
310  for (k = 0; k < np; k++)
311  for (q = 0; q < 3; q++)
312  cm[q] += rr(k,q);
313 
314  if (np > 0) {
315  for (q = 0; q < 3; q++)
316  cm[q] = cm[q]/np;
317 
318  for (k = 0, ave = 0.0; k < np; k++) {
319  for (q = 0; q < 3; q++)
320  diff[q] = rr(k,q) - cm[q];
321  ave += sqrt(pow(diff(0),2) + pow(diff(1),2) + pow(diff(2),2));
322  }
323  avep = ave/np;
324  }
325 }
326 
327 
328 //*************************************************************************************************************
329 
330 MatrixXf LayoutMaker::make_initial_simplex(VectorXf &pars,
331  int npar,
332  float size)
333 {
334  /*
335  * Make the initial tetrahedron
336  */
337  MatrixXf simplex(npar+1,npar);
338  int k;
339 
340  for (k = 0; k < npar+1; k++)
341  simplex.row(k) = pars;
342 
343  for (k = 1; k < npar+1; k++)
344  simplex(k,k-1) = simplex(k,k-1) + size;
345 
346  return simplex;
347 }
348 
349 
350 //*************************************************************************************************************
351 
352 int LayoutMaker::fit_sphere_to_points(MatrixXf &rr,
353  int np,
354  float simplex_size,
355  VectorXf &r0,
356  float &R)
357 {
358  /*
359  * Find the optimal sphere origin
360  */
361  fitUserRec user;
362  float ftol = 1e-3;
363  int max_eval = 5000;
364  int report_interval = -1;
365  int neval;
366  MatrixXf init_simplex;
367  VectorXf init_vals(4);
368 
369  VectorXf cm(3);
370  float R0;
371  int k;
372 
373  int res = FAIL;
374 
375  user.rr = rr;
376  user.np = np;
377 
378  R0 = 0.1;
379  calculate_cm_ave_dist(rr,np,cm,R0);
380 
381  init_simplex = make_initial_simplex(cm,3,simplex_size);
382 
383  std::cout << "sphere origin calcuated" << cm[0] << " " << cm[1] << " " << cm[2] << std::endl;
384 
385  user.report = FALSE;
386 
387  for (k = 0; k < 4; k++)
388  init_vals[k] = fit_eval(static_cast<VectorXf>(init_simplex.row(k)),3,&user);
389 
390  user.report = FALSE;
391 
392  //Start the minimization
393  if(MinimizerSimplex::mne_simplex_minimize(init_simplex, /* The initial simplex */
394  init_vals, /* Function values at the vertices */
395  3, /* Number of variables */
396  ftol, /* Relative convergence tolerance */
397  fit_eval, /* The function to be evaluated */
398  &user, /* Data to be passed to the above function in each evaluation */
399  max_eval, /* Maximum number of function evaluations */
400  neval, /* Number of function evaluations */
401  report_interval, /* How often to report (-1 = no_reporting) */
402  report_func) != OK) /* The function to be called when reporting */
403  return FALSE;
404 
405  r0[0] = init_simplex(0,0);
406  r0[1] = init_simplex(0,1);
407  r0[2] = init_simplex(0,2);
408  R = opt_rad(r0,&user);
409 
410  res = OK;
411 
412  return res;
413 }
LayoutLoader class declaration.
static bool makeLayout(const QList< QVector< double > > &inputPoints, QList< QVector< double > > &outputPoints, const QStringList &names, QFile &outFile, bool do_fit, float prad, float w, float h, bool writeFile=false)
Definition: layoutmaker.cpp:73