Lemma is an Electromagnetics API
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634
  1. #ifndef __MATPLOT_VTK
  2. #define __MATPLOT_VTK
  3. /*
  4. * MatPlot_VTK
  5. *
  6. * Simple plotting of vectors and matrices based on
  7. * the VTK libraries.
  8. *
  9. * Four classes:
  10. * Plot2D_VTK - 2d plotting akin to plot(x,y)
  11. * Surf_VTK - Surface plotting akin to surf(x,y,z)
  12. * Contour_VTK - Contour plotting
  13. * Quiver_VTK - Vector field plot
  14. *
  15. * See examples.cpp for usage instructions.
  16. *
  17. * These classes are released "as is", without any implied
  18. * warranty or fitness for any particular purpose.
  19. *
  20. * Dag Lindbo, dag@csc.kth.se, 2007-12-09
  21. */
  22. // matrix size functions
  23. // #include "ublas_dims.h"
  24. //#include "mtl4_dims.h"
  25. // system includes
  26. #include <assert.h>
  27. #include <string>
  28. #include <cmath>
  29. // vtk includes used by templates
  30. #include "vtkFloatArray.h"
  31. #include "vtkPointData.h"
  32. #include "vtkPoints.h"
  33. #include "vtkRectilinearGrid.h"
  34. #include "vtkRenderer.h"
  35. #include "vtkStructuredGrid.h"
  36. #include "vtkXYPlotActor.h"
  37. namespace matplot {
  38. enum SCALE {LINEAR, LOG10};
  39. void render_interactive(vtkRenderer *rend, int xPix, int yPix);
  40. void render_interactive_camera(vtkRenderer *rend, int xPix, int yPix,
  41. vtkCamera* mycam);
  42. void render_interactive_cam(vtkRenderer *rend, int xPix, int yPix,
  43. double cam[3], double focal[3]);
  44. void render_to_png(vtkRenderer *rend, int xPix, int yPix,std::string fname);
  45. void render_to_png_cam(vtkRenderer *rend, int xPix, int yPix,
  46. std::string fname, double cam[3], double focal[3]);
  47. /** Plot vectors x versus y.
  48. * Inspired by "plot" in Matlab.
  49. *
  50. * Notes:
  51. *
  52. * Call plot(x,y,<...>) multiple times before show() to
  53. * get multiple curves in one figure.
  54. *
  55. * \author Dag Lindbo
  56. */
  57. class Plot2D_VTK {
  58. public:
  59. Plot2D_VTK(std::string x_label="x", std::string y_label="y",
  60. int xpix = 800, int ypix = 600, bool semilogx=false );
  61. ~Plot2D_VTK();
  62. /** Insert curve x vs.\ y into plot.
  63. * - Default color is blue, {0 , 0, 1}.
  64. * - Default line style is solid line, "-".
  65. */
  66. template<typename Vec_t>
  67. void plot(const Vec_t &x, const Vec_t &y) {
  68. double color[3] = { 0, 0, 1.0 };
  69. plot(x, y, color, "-");
  70. }
  71. /** Insert curve x vs.\ y into plot (specify color and line style).
  72. * Color specification:
  73. * - RGB (normalized decimal), {0.32, 0.1, 0.0}
  74. *
  75. * Line style specification:
  76. * - "-" solid line
  77. * - "." dotted line
  78. * - ".-" or "-." solid line with dots
  79. */
  80. template<typename Vec_t>
  81. void plot(const Vec_t &x, const Vec_t &y, const double col[3],
  82. const std::string linespec) {
  83. int i, N = x.size(), plot_points = 0, plot_lines = 0;
  84. vtkRectilinearGrid *curve;
  85. vtkFloatArray *yVal;
  86. vtkFloatArray *xVal;
  87. // determine line style
  88. if (linespec == "-")
  89. plot_lines = 1;
  90. else if (linespec == ".")
  91. plot_points = 1;
  92. else if (linespec == ".-" || linespec == "-.") {
  93. plot_points = 1;
  94. plot_lines = 1;
  95. }
  96. // put (x,y) into VTK FloatArrays
  97. xVal = vtkFloatArray::New();
  98. yVal = vtkFloatArray::New();
  99. for (i=0; i<N; i++) {
  100. xVal->InsertNextTuple1(x[i]);
  101. yVal->InsertNextTuple1(y[i]);
  102. }
  103. // Make a VTK Rectlinear grid from arrays
  104. curve = vtkRectilinearGrid::New();
  105. curve->SetDimensions(N, 1, 1);
  106. curve->SetXCoordinates(xVal);
  107. curve->GetPointData()->SetScalars(yVal);
  108. // attach gridfunction to plot
  109. xyplot->AddDataSetInput(curve);
  110. if (semilogx) {
  111. xyplot->LogxOn();
  112. }
  113. // VTK doesn't have this? :(
  114. //xyplot->LogyOn();
  115. // how to read data
  116. xyplot->SetXValuesToValue();
  117. // set attributes
  118. xyplot->SetPlotColor(plot_no, col[0], col[1], col[2]);
  119. xyplot->SetPlotLines(plot_no, plot_lines);
  120. xyplot->SetPlotPoints(plot_no, plot_points);
  121. plot_no++;
  122. xVal->Delete();
  123. yVal->Delete();
  124. curve->Delete();
  125. }
  126. void show();
  127. void draw_to_png(std::string filename);
  128. private:
  129. int xPix;
  130. int yPix;
  131. bool semilogx;
  132. int plot_no;
  133. vtkRenderer* rend;
  134. vtkXYPlotActor* xyplot;
  135. };
  136. /** Plot z = f(x,y) surface.
  137. * Inspired by "surf" in Matlab
  138. *
  139. * \author Dag Lindbo
  140. */
  141. class Surf_VTK
  142. {
  143. public:
  144. Surf_VTK(int px = 800, int py = 600);
  145. ~Surf_VTK();
  146. /** Create surface plot.
  147. * Matrix z, with respect to vectors x and y.
  148. */
  149. template<typename Vec_t, typename Mat_t>
  150. void surf(const Vec_t &x, const Vec_t &y, const Mat_t &z)
  151. {
  152. geometry(x, y, z);
  153. renderer(false, true, true, false);
  154. render_interactive(rend, xPix, yPix);
  155. }
  156. /** Create surface plot.
  157. * Matrix z, with respect to vectors x and y.
  158. *
  159. * Warp z-axis to produce better fit:
  160. * - do_warp = true
  161. */
  162. template<typename Vec_t, typename Mat_t>
  163. void surf(const Vec_t &x, const Vec_t &y, const Mat_t &z, bool do_warp)
  164. {
  165. geometry(x, y, z);
  166. renderer(false, true, true, do_warp);
  167. render_interactive(rend, xPix, yPix);
  168. }
  169. /** Create surface plot.
  170. * Matrix z, with respect to vectors x and y.
  171. *
  172. * Warp z-axis to produce better fit:
  173. * - do_warp = true
  174. *
  175. * Camera control:
  176. * - specify camera position: cam = { -15.0, 10.0, 12.0 }
  177. * - specify focal point: focal = { 0, 0, 0 }
  178. */
  179. template<typename Vec_t, typename Mat_t>
  180. void surf(const Vec_t &x, const Vec_t &y, const Mat_t &z, bool do_warp,
  181. double observer[3], double focal[3])
  182. {
  183. geometry(x, y, z);
  184. renderer(false, true, true, do_warp);
  185. render_interactive_cam(rend, xPix, yPix, observer, focal);
  186. }
  187. /** Create surface plot and render to file
  188. * Matrix z, with respect to vectors x and y.
  189. *
  190. * Warp z-axis to produce better fit:
  191. * - do_warp = true
  192. *
  193. * Camera control:
  194. * - specify camera position: cam = { -15.0, 10.0, 12.0 }
  195. * - specify focal point: focal = { 0, 0, 0 }
  196. */
  197. template<typename Vec_t, typename Mat_t>
  198. void surf_to_file(const Vec_t &x, const Vec_t &y, const Mat_t &z,
  199. bool do_warp, std::string fname, double observer[3],
  200. double focal[3])
  201. {
  202. geometry(x, y, z);
  203. renderer(false, true, true, do_warp);
  204. render_to_png_cam(rend, xPix, yPix, fname, observer, focal);
  205. }
  206. void purge();
  207. private:
  208. vtkStructuredGrid *gridfunc;
  209. vtkRenderer *rend;
  210. double Lxy, Lz;
  211. bool has_data;
  212. int xPix, yPix;
  213. template<typename Vec_t, typename Mat_t>
  214. void geometry(const Vec_t &x, const Vec_t &y, const Mat_t &z)
  215. {
  216. const unsigned int Nx = x.size(); //vec_dim(x);
  217. const unsigned int Ny = y.size(); //vec_dim(y);
  218. unsigned int i, j, k;
  219. // make sure the input is ok and that this surfaceplot is free
  220. assert(Nx == z.rows() );
  221. assert(Ny == z.cols() );
  222. assert(!has_data);
  223. // determine x-y range of data
  224. if (x(Nx-1)-x(0) > y(Ny-1)-y(0))
  225. Lxy = x(Nx-1)-x(0);
  226. else
  227. Lxy = y(Ny-1)-y(0);
  228. double z_low = 10000, z_upp = -10000;
  229. // put data, z, into a 2D structured grid
  230. gridfunc->SetDimensions(Nx, Ny, 1);
  231. vtkPoints *points = vtkPoints::New();
  232. for (j = 0; j < Ny; j++)
  233. {
  234. for (i = 0; i < Nx; i++)
  235. {
  236. points->InsertNextPoint(x(i), y(j), z(i, j));
  237. if (z(i, j)< z_low)
  238. z_low = z(i, j);
  239. if (z(i, j)> z_upp)
  240. z_upp = z(i, j);
  241. }
  242. }
  243. gridfunc->SetPoints(points);
  244. // get scalar field from z-values
  245. vtkFloatArray *colors = vtkFloatArray::New();
  246. colors->SetNumberOfComponents(1);
  247. colors->SetNumberOfTuples(Nx*Ny);
  248. k = 0;
  249. for (j = 0; j < Ny; j++)
  250. for (i = 0; i < Nx; i++)
  251. {
  252. colors->InsertComponent(k, 0, z(i, j));
  253. k++;
  254. }
  255. gridfunc->GetPointData()->SetScalars(colors);
  256. points->Delete();
  257. colors->Delete();
  258. has_data = true;
  259. Lz = z_upp-z_low;
  260. }
  261. void renderer(bool, bool, bool, bool);
  262. };
  263. /** Plot contour lines for a function in the plane.
  264. * Inspired by "contour" in Matlab
  265. *
  266. * \author Dag Lindbo
  267. */
  268. class Contour_VTK {
  269. public:
  270. Contour_VTK(int px = 800, int py = 600);
  271. // specify linear or log on bars
  272. Contour_VTK(int px = 800, int py = 600, SCALE xscale=LINEAR, SCALE
  273. yscale=LINEAR);
  274. ~Contour_VTK();
  275. /** Create contour plot.
  276. * Matrix z vs. vectors x and y.
  277. * Produces a default number of contour lines (10) and
  278. * colors the lines instead of the underlying surface.
  279. */
  280. template<typename Vec_t, typename Mat_t>
  281. void contour(const Vec_t &x, const Vec_t &y, const Mat_t &z) {
  282. geometry(x, y, z);
  283. renderer(true, false, 10);
  284. render_interactive(rend, xPix, yPix);
  285. }
  286. /**Create contour plot.
  287. * Matrix z vs. vectors x and y.
  288. *
  289. * Number of contour lines:
  290. * - num_lines
  291. *
  292. * Coloring:
  293. * - draw_surf = false: color contour lines and omits underlying surface
  294. * - draw_surf = true: draw contour lines white and color the
  295. * underlying surface
  296. */
  297. template<typename Vec_t, typename Mat_t>
  298. void contour(const Vec_t &x, const Vec_t &y, const Mat_t &z,
  299. bool draw_surf, int num_lines) {
  300. geometry(x, y, z);
  301. renderer(true, draw_surf, num_lines);
  302. render_interactive(rend, xPix, yPix);
  303. }
  304. /**Create contour plot and render to file.
  305. * Matrix z vs. vectors x and y.
  306. *
  307. * Number of contour lines:
  308. * - num_lines
  309. *
  310. * Coloring:
  311. * - draw_surf = false: color contour lines and omits underlying surface
  312. * - draw_surf = true: draw contour lines white and color the
  313. * underlying surface
  314. */
  315. template<typename Vec_t, typename Mat_t>
  316. void contour_to_file(const Vec_t &x, const Vec_t &y, const Mat_t &z,
  317. bool draw_surf, int num_lines, std::string fname)
  318. {
  319. geometry(x, y, z);
  320. renderer(true, draw_surf, num_lines);
  321. render_to_png(rend, xPix, yPix, fname);
  322. }
  323. void purge();
  324. void SetXLabel(const std::string &xlab);
  325. void SetYLabel(const std::string &ylab);
  326. private:
  327. SCALE XScale;
  328. SCALE YScale;
  329. vtkRectilinearGrid *gridfunc;
  330. vtkRenderer *rend;
  331. bool has_data;
  332. int xPix, yPix;
  333. double axscale;
  334. double ymin;
  335. double ymax;
  336. std::string xlabel;
  337. std::string ylabel;
  338. template<typename Vec_t, typename Mat_t>
  339. void geometry(const Vec_t &x, const Vec_t &y, const Mat_t &z) {
  340. const unsigned int Nx = x.size();
  341. const unsigned int Ny = y.size();
  342. unsigned int i, j, k;
  343. // make sure the input is ok and that this contourplot is free
  344. assert(Nx == z.rows() );
  345. assert(Ny == z.cols() );
  346. assert(!has_data);
  347. // x and y vectors go into vtkFloatArray
  348. vtkFloatArray *xcoord = vtkFloatArray::New();
  349. xcoord->SetNumberOfComponents(1);
  350. xcoord->SetNumberOfTuples(Nx);
  351. vtkFloatArray *ycoord = vtkFloatArray::New();
  352. ycoord->SetNumberOfComponents(1);
  353. ycoord->SetNumberOfTuples(Ny);
  354. // We want the two axis to be equal, not squashed
  355. // normalize axis ratio
  356. axscale = 1.;
  357. ymin = y.minCoeff();
  358. ymax = y.maxCoeff();
  359. if (YScale == LINEAR && XScale == LINEAR)
  360. axscale = (x.maxCoeff() - x.minCoeff()) /
  361. (y.maxCoeff() - y.minCoeff()) ;
  362. if (YScale == LOG10 && XScale == LINEAR)
  363. axscale = ( x.maxCoeff() - x.minCoeff() ) /
  364. (std::log10(y.maxCoeff()) - std::log10(y.minCoeff())) ;
  365. if (YScale == LOG10 && XScale == LOG10)
  366. axscale = ( (std::log10(x.maxCoeff()) - std::log10(x.minCoeff())) /
  367. (std::log10(y.maxCoeff()) - std::log10(y.minCoeff())) );
  368. if (YScale == LINEAR && XScale == LOG10)
  369. axscale = (std::log10(x.maxCoeff()) - std::log10(x.minCoeff())) /
  370. ( y.maxCoeff() - y.minCoeff() ) ;
  371. if (XScale == LINEAR) {
  372. for (i=0; i<Nx; i++)
  373. xcoord->InsertComponent(i, 0, x(i));
  374. } else {
  375. for (i=0; i<Nx; i++)
  376. xcoord->InsertComponent(i, 0, std::log10(x(i)));
  377. }
  378. if (YScale == LINEAR) {
  379. for (i=0; i<Ny; i++)
  380. ycoord->InsertComponent(i, 0, axscale*y(i));
  381. } else {
  382. for (i=0; i<Ny; i++)
  383. ycoord->InsertComponent(i, 0, axscale*std::log10(y(i)));
  384. }
  385. // Create rectilinear grid
  386. gridfunc->SetDimensions(Nx, Ny, 1);
  387. gridfunc->SetXCoordinates(xcoord);
  388. gridfunc->SetYCoordinates(ycoord);
  389. // add z-values as scalars to grid
  390. vtkFloatArray *colors = vtkFloatArray::New();
  391. colors->SetNumberOfComponents(1);
  392. colors->SetNumberOfTuples(Nx*Ny);
  393. k = 0;
  394. for (j = 0; j < Ny; j++)
  395. for (i = 0; i < Nx; i++) {
  396. colors->InsertComponent(k, 0, z(i, j));
  397. k++;
  398. }
  399. gridfunc->GetPointData()->SetScalars(colors);
  400. colors->Delete();
  401. xcoord->Delete();
  402. ycoord->Delete();
  403. has_data = true;
  404. }
  405. void renderer(bool, bool, int);
  406. };
  407. /** Plot vector-valued function in the plane.
  408. * Inspired by "quiver" in Matlab
  409. *
  410. * \author Dag Lindbo
  411. */
  412. class Quiver_VTK
  413. {
  414. public:
  415. Quiver_VTK(int px = 800, int py = 600);
  416. ~Quiver_VTK();
  417. /** Create vector arrow plot (quiver).
  418. * Pointwise vecotrs in matrices u and v, at grid
  419. * points given by vectors x and y. Color by magnitude.
  420. * Defaults to no scaling of arrow lengths.
  421. */
  422. template<typename Vec_t, typename Mat_t>
  423. void quiver(const Vec_t &x, const Vec_t &y, const Mat_t &u,
  424. const Mat_t &v)
  425. {
  426. geometry(x, y, u, v);
  427. renderer(1.0);
  428. render_interactive(rend, xPix, yPix);
  429. }
  430. /** Create vector arrow plot (quiver).
  431. * Pointwise vectors in matrices u and v, at grid
  432. * points given by vectors x and y. Color by magnitude.
  433. *
  434. * Scales arrows by a factor s.
  435. */
  436. template<typename Vec_t, typename Mat_t>
  437. void quiver(const Vec_t &x, const Vec_t &y, const Mat_t &u,
  438. const Mat_t &v, double s)
  439. {
  440. geometry(x, y, u, v);
  441. renderer(s);
  442. render_interactive(rend, xPix, yPix);
  443. }
  444. /** Create vector arrow plot (quiver) and render to file.
  445. * Pointwise vectors in matrices u and v, at grid
  446. * points given by vectors x and y. Color by magnitude
  447. *
  448. * Scales arrows by a factor s.
  449. */
  450. template<typename Vec_t, typename Mat_t>
  451. void quiver_to_file(const Vec_t &x, const Vec_t &y, const Mat_t &u,
  452. const Mat_t &v, double s, std::string filename)
  453. {
  454. geometry(x, y, u, v);
  455. renderer(s);
  456. render_to_png(rend, xPix, yPix, filename);
  457. }
  458. void purge();
  459. private:
  460. vtkRectilinearGrid *gridfunc;
  461. vtkRenderer *rend;
  462. bool has_data;
  463. int xPix, yPix;
  464. template<typename Vec_t, typename Mat_t>
  465. void geometry(const Vec_t& x, const Vec_t& y, const Mat_t& u,
  466. const Mat_t& v)
  467. {
  468. const unsigned int Nx = x.size(); //vec_dim(x);
  469. const unsigned int Ny = y.size(); //vec_dim(y);
  470. unsigned int i, j, k;
  471. // make sure the input is ok and that this contourplot is free
  472. assert(Nx == u.rows());
  473. assert(Ny == u.cols());
  474. assert(Nx == v.rows());
  475. assert(Ny == v.cols());
  476. assert(!has_data);
  477. // x and y vectors go into vtkFloatArray
  478. vtkFloatArray *xcoord = vtkFloatArray::New();
  479. xcoord->SetNumberOfComponents(1);
  480. xcoord->SetNumberOfTuples(Nx);
  481. vtkFloatArray *ycoord = vtkFloatArray::New();
  482. ycoord->SetNumberOfComponents(1);
  483. ycoord->SetNumberOfTuples(Ny);
  484. for (i=0; i<Nx; i++)
  485. xcoord->InsertComponent(i, 0, x(i));
  486. for (i=0; i<Ny; i++)
  487. ycoord->InsertComponent(i, 0, y(i));
  488. // Create rectilinear grid
  489. gridfunc->SetDimensions(Nx, Ny, 1);
  490. gridfunc->SetXCoordinates(xcoord);
  491. gridfunc->SetYCoordinates(ycoord);
  492. // add magnitude of (u,v) as scalars to grid
  493. vtkFloatArray *colors = vtkFloatArray::New();
  494. colors->SetNumberOfComponents(1);
  495. colors->SetNumberOfTuples(Nx*Ny);
  496. // add vector (u,v) to grid
  497. vtkFloatArray *vectors = vtkFloatArray::New();
  498. vectors->SetNumberOfComponents(3);
  499. vectors->SetNumberOfTuples(Nx*Ny);
  500. k = 0;
  501. for (j = 0; j < Ny; j++)
  502. for (i = 0; i < Nx; i++)
  503. {
  504. colors->InsertTuple1(k,sqrt(u(i,j)*u(i,j)+v(i,j)*v(i,j)));
  505. vectors->InsertTuple3(k, u(i, j), v(i, j), 0.0);
  506. k++;
  507. }
  508. gridfunc->GetPointData()->SetScalars(colors);
  509. gridfunc->GetPointData()->SetVectors(vectors);
  510. vectors->Delete();
  511. colors->Delete();
  512. xcoord->Delete();
  513. ycoord->Delete();
  514. has_data = true;
  515. }
  516. void renderer(double);
  517. };
  518. }
  519. #endif