Context Navigation

lscereader.cpp @ 241

Last change on this file since 241 was 241, checked in by hozdoba, 11 years ago

File size: 38.0 KB

Line
1
2	#ifndef LSCE_EXPORTS
3	# include "lscereader.hpp"
4	#else
5	# include "vtkLSCEReader.h"
6	#endif //LSCE_EXPORTS
7
8
9	// BibliothÃšque VTK
10	#include "vtkObjectFactory.h"
11	#include "vtkInformation.h"
12	#include "vtkInformationVector.h"
13	#include "vtkStreamingDemandDrivenPipeline.h"
14
15	#include "vtkPointLocator.h"
16
17	#include <vtkProperty.h>
18	#include <vtkDataSetMapper.h>
19	#include <vtkActor.h>
20	#include <vtkRenderWindow.h>
21	#include <vtkRenderer.h>
22	#include <vtkRenderWindowInteractor.h>
23	#include <vtkLegendScaleActor.h>
24
25	#include "vtkMath.h"
26
27	#include "inetcdf4_impl.hpp"
28	#include "inetcdf4_adv_impl.hpp"
29
30	// BibliothÃšque C++
31	#include <algorithm>
32
33	// Type Fortran
34	#define ARRAY(valuetype, numdims) boost::shared_ptr<boost::multi_array<valuetype, numdims> >
35
36	#define ARRAY_ASSIGN(value, valuetype, numdims, extent)\
37	value.reset(new boost::multi_array<valuetype, numdims>(boost::extents extent))
38
39	#define ARRAY_CREATE(value, valuetype, numdims, extent)\
40	ARRAY(valuetype, numdims) value = \
41	ARRAY(valuetype, numdims)(new boost::multi_array<valuetype, numdims>(boost::extents extent))
42
43	// Type C
44	#define ARRAY_C_ASSIGN(value, valuetype, numdims, extent)\
45	value = ARRAY(valuetype, numdims) \
46	(new boost::multi_array<valuetype, numdims>(boost::extents extent, c_storage_order()))
47
48	#define ARRAY_C_CREATE(value, valuetype, numdims, extent)\
49	ARRAY_C_ASSIGN(ARRAY(valuetype, numdims) value, valuetype, numdims, extent)
50
51
52	#ifndef LSCE_EXPORTS
53	namespace xmlioserver
54	{
55	namespace vtk
56	{
57	#endif //LSCE_EXPORTS
58
59	/// ////////////////////// DÃ©finitions ////////////////////// ///
60	vtkStandardNewMacro(vtkLSCEReader);
61
62	vtkLSCEReader::vtkLSCEReader(void)
63	: FileName(), CurGridType(RECTILINEAR), VarNames()
64	, A3D(true), ATemporal(true), ACell(true), AGridDef(false)
65	, Reader()
66	{
67	this->SetNumberOfInputPorts (0);
68	this->SetNumberOfOutputPorts(1);
69	}
70
71	vtkLSCEReader::~vtkLSCEReader(void)
72	{
73
74	}
75
76	///---------------------------------------------------------------
77
78	void vtkLSCEReader::SetFileName(const vtkStdString & fileName)
79	{
80	this->FileName = fileName;
81	this->Reader = boost::shared_ptr<xmlioserver::io::CINetCDF4Adv>
82	(new xmlioserver::io::CINetCDF4Adv(fileName));
83	this->Modified();
84	}
85
86	void vtkLSCEReader::SetVariable(const vtkStdString & variable)
87	{
88	this->RemoveAllSelectedVariables();
89	if (Reader->hasVariable(variable))
90	this->AddVariableToSelection(variable);
91	this->Modified();
92	}
93
94	const char * vtkLSCEReader::GetFileName(void) const
95	{
96	return (this->FileName.c_str());
97	}
98
99	void vtkLSCEReader::PrintSelf(ostream& os, vtkIndent indent)
100	{
101	this->Superclass::PrintSelf(os,indent);
102	}
103
104	//----------------------------------------------------------------
105	void vtkLSCEReader::AddVariableToSelection(const vtkStdString & varName)
106	{
107	if (this->Reader->isRectilinear(varName))
108	this->SetGridType (RECTILINEAR);
109	else if (this->Reader->isCurvilinear(varName))
110	this->SetGridType (CURVILINEAR);
111	else
112	this->SetGridType (UNSTRUCTURED);
113
114	this->VarNames.insert(varName);
115	}
116
117	void vtkLSCEReader::RemoveSelectedVariable(const vtkStdString & varName)
118	{
119	std::set<vtkStdString>::iterator it = this->VarNames.find(varName);
120	if (it != this->VarNames.end()) this->VarNames.erase(it);
121	}
122
123	void vtkLSCEReader::RemoveAllSelectedVariables(void)
124	{
125	this->VarNames.clear();
126	}
127
128	const std::set<vtkStdString> &
129	vtkLSCEReader::GetSelectedVariables(void) const
130	{
131	return (this->VarNames);
132	}
133
134	bool vtkLSCEReader::HasSelectedVariable(void) const
135	{
136	return (this->VarNames.size() != 0);
137	}
138
139	//----------------------------------------------------------------
140
141	void vtkLSCEReader::SetGridType (GridType type)
142	{
143	if ((this->IsUnstructured() && (type == CURVILINEAR)) \|\|
144	(this->IsUnstructured() && (type == RECTILINEAR)) \|\|
145	(this->IsCurvilinear() && (type == RECTILINEAR)))
146	this->RemoveAllSelectedVariables();
147	this->CurGridType = type;
148	}
149
150	bool vtkLSCEReader::IsUnstructured(void) const
151	{
152	return (this->CurGridType == UNSTRUCTURED);
153	}
154
155	bool vtkLSCEReader::IsCurvilinear(void) const
156	{
157	return (this->CurGridType == CURVILINEAR);
158	}
159
160	bool vtkLSCEReader::IsRectilinear(void) const
161	{
162	return (this->CurGridType == RECTILINEAR);
163	}
164
165	//----------------------------------------------------------------
166
167	void vtkLSCEReader::AcceptTemporalOnly(bool value)
168	{
169	this->ATemporal = value;
170	}
171
172	void vtkLSCEReader::Accept3DOnly(bool value)
173	{
174	this->A3D = value;
175	}
176
177	void vtkLSCEReader::AcceptCellOnly(bool value)
178	{
179	this->ACell = value;
180	}
181
182	//----------------------------------------------------------------
183
184	int vtkLSCEReader::ProcessRequest
185	(vtkInformation * request,
186	vtkInformationVector ** vtkNotUsed(inputVector),
187	vtkInformationVector * outputVector)
188	{
189	try
190	{
191
192	if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
193	{
194	return (this->RequestData(NULL, NULL, outputVector));
195	}
196
197	if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_OBJECT()))
198	{
199	// return (this->RequestDataObject(NULL, NULL, outputVector));
200	}
201
202	if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
203	{
204	return (this->RequestInformation(NULL, NULL, outputVector));
205	}
206
207	if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
208	{
209	return (this->RequestUpdateExtent(NULL, NULL, outputVector));
210	}
211	}
212	catch (xmlioserver::CException & _exception)
213	{
214	std::cerr << _exception.getMessage() << std::endl;
215	return (-1);
216	}
217	return (1);
218	}
219
220	//----------------------------------------------------------------
221
222	void vtkLSCEReader::ShowVariable(const vtkStdString & filename,
223	const vtkStdString & varname)
224	{
225	vtkSmartPointer<vtkLSCEReader> lscereader =
226	vtkSmartPointer<vtkLSCEReader>::New();
227
228	lscereader->SetFileName(filename);
229	lscereader->AddVariableToSelection(varname);
230	lscereader->Update();
231
232	// CrÃ©ation d'un mappeur et d'un acteur.
233	vtkSmartPointer<vtkDataSetMapper> mapper = vtkSmartPointer<vtkDataSetMapper>::New();
234	mapper->SetInputConnection(lscereader->GetOutput()->GetProducerPort());
235	mapper->SetScalarRange (1000, 10000);
236
237	vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
238	actor->SetMapper(mapper);
239	actor->GetProperty()->SetRepresentationToWireframe();
240
241	// CrÃ©ation d'une fenÃªtre de rendu.
242	vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
243	vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
244	renderWindow->AddRenderer(renderer);
245	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
246	vtkSmartPointer<vtkRenderWindowInteractor>::New();
247	renderWindowInteractor->SetRenderWindow(renderWindow);
248
249
250	vtkSmartPointer<vtkLegendScaleActor> legendScaleActor =
251	vtkSmartPointer<vtkLegendScaleActor>::New();
252
253	// Ajout de l'acteur Ã la scÃšne.
254	renderer->AddActor(actor);
255	renderer->AddActor(legendScaleActor);
256
257	renderer->GradientBackgroundOn();
258	renderer->SetBackground(1,1,1);
259	renderer->SetBackground2(0,0,0);
260
261	// Rendu et intÃ©raction
262	renderWindow->Render();
263	renderWindowInteractor->Start();
264	}
265
266	//----------------------------------------------------------------
267
268	void vtkLSCEReader::CreateRectilinearGrid(vtkRectilinearGrid * grid,
269	vtkInformation * outInfo,
270	vtkFloatArray * xspacing,
271	vtkFloatArray * yspacing,
272	vtkFloatArray * zspacing,
273	vtkIntArray * dimensions)
274	{
275	int extent[6] = { 0, dimensions->GetValue(0) - 1
276	, 0, dimensions->GetValue(1) - 1
277	, 0, dimensions->GetValue(2) - 1};
278	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
279	grid->SetExtent (extent);
280	grid->SetDimensions (dimensions->GetValue(0),
281	dimensions->GetValue(1),
282	dimensions->GetValue(2));
283	grid->SetXCoordinates (xspacing);
284	grid->SetYCoordinates (yspacing);
285	grid->SetZCoordinates (zspacing);
286	}
287
288	void vtkLSCEReader::CreateStructuredGrid(vtkStructuredGrid * grid,
289	vtkInformation * outInfo,
290	vtkPoints * points,
291	vtkIntArray * dimensions)
292	{
293	int extent[6] = { 0, dimensions->GetValue(0) - 1
294	, 0, dimensions->GetValue(1) - 1
295	, 0, dimensions->GetValue(2) - 1};
296	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
297	grid->SetExtent (extent);
298	grid->SetDimensions (dimensions->GetValue(0),
299	dimensions->GetValue(1),
300	dimensions->GetValue(2));
301	grid->SetPoints(points);
302	}
303
304	void vtkLSCEReader::CreateUnstructuredGrid(vtkUnstructuredGrid * grid,
305	vtkInformation * outInfo,
306	vtkPoints * points,
307	vtkCellArray * cells,
308	int celltype)
309	{
310	grid->SetPoints(points);
311	grid->SetCells(celltype, cells);
312	}
313
314	//----------------------------------------------------------------
315
316	void vtkLSCEReader::GetRectilinearConnectivity(int dimZ, int dimY, int dimX, vtkCellArray * cells)
317	{
318	if (dimZ == 1 && dimY > 1 && dimX > 1) // 2D
319	{
320	vtkIdType pts[4] ;
321
322	for (int i = 0, k = 0; i < (dimX - 1); i++)
323	for (int j = 0; j < (dimY - 1); j++, k++)
324	{
325	pts[0] = i * dimY + j;
326	pts[1] = i * dimY + j + 1;
327	pts[2] = (i+1) * dimY + j + 1;
328	pts[3] = (i+1) * dimY + j;
329	cells->InsertNextCell(4, pts);
330	}
331	}
332	else if (dimZ > 1 && dimY > 1 && dimX > 1) // 3D
333	{
334	vtkIdType pts[8];
335	for (int i = 0, l = 0; i < (dimX - 1); i++)
336	for (int j = 0; j < (dimY - 1); j++)
337	for (int k = 0; k < (dimZ - 1); k++, l++)
338	{
339	pts[0] = i * dimY * dimZ + j * dimZ + k;
340	pts[1] = i * dimY * dimZ + j * dimZ + k + 1;
341	pts[2] = (i+1) * dimY * dimZ + j * dimZ + k + 1;
342	pts[3] = (i+1) * dimY * dimZ + j * dimZ + k;
343
344	pts[4] = i * dimY * dimZ + (j + 1) * dimZ + k ;
345	pts[5] = i * dimY * dimZ + (j + 1) * dimZ + k +1;
346	pts[6] = (i+1) * dimY * dimZ + (j + 1) * dimZ + k + 1;
347	pts[7] = (i+1) * dimY * dimZ + (j + 1) * dimZ + k ;
348	cells->InsertNextCell(8, pts);
349	}
350	}
351	}
352
353	//----------------------------------------------------------------
354
355	void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
356	int yi, int yf,
357	int zi, int zf,
358	vtkFloatArray * xspacing,
359	vtkFloatArray * yspacing,
360	vtkFloatArray * zspacing,
361	vtkIntArray * dimensions)
362	{
363	dimensions->InsertNextValue (xf - xi + 1); // X
364	dimensions->InsertNextValue (yf - yi + 1); // Y
365	dimensions->InsertNextValue (zf - zi + 1); // Z
366	for (int i = xi; i <= xf; i++)
367	xspacing->InsertNextValue (i);
368	for (int j = yi; j <= yf; j++)
369	yspacing->InsertNextValue (j);
370	for (int k = zi; k <= zf; k++)
371	zspacing->InsertNextValue (k);
372	}
373
374	void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
375	int yi, int yf,
376	int zi, int zf,
377	vtkPoints * points,
378	vtkIntArray * dimensions)
379	{
380	dimensions->InsertNextValue (xf - xi + 1); // X
381	dimensions->InsertNextValue (yf - yi + 1); // Y
382	dimensions->InsertNextValue (zf - zi + 1); // Z
383	for (int i = xi; i <= xf; i++)
384	for (int j = yi; j <= yf; j++)
385	for (int k = zi; k <= zf; k++)
386	points->InsertNextPoint(i, j, k);
387	}
388
389	void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
390	int yi, int yf,
391	int zi, int zf,
392	vtkPoints * points,
393	vtkCellArray * cells,
394	vtkIntArray * dimensions)
395	{
396	this->CreateSimpleGrid (xi, xf, yi, yf, zi, zf, points, dimensions);
397	int dimX = dimensions->GetValue(0),
398	dimY = dimensions->GetValue(1),
399	dimZ = dimensions->GetValue(2);
400	this->GetRectilinearConnectivity(dimZ, dimY, dimX, cells);
401	}
402
403	//----------------------------------------------------------------
404	void vtkLSCEReader::GetSpacings(const vtkStdString & coordinate,
405	bool bounds, vtkFloatArray * spacing)
406	{
407	std::size_t size = (Reader->getDimensions(&coordinate).begin())->second;
408	ARRAY_CREATE(data_arr, float, 1, [1]);
409	if (Reader->hasVariable(coordinate))
410	{
411	if (bounds)
412	{
413	if (Reader->hasBounds(coordinate))
414	{ // CellData, variable de coordonnÃ©e avec attributes bounds.
415	vtkStdString boundid = Reader->getBoundsId(coordinate);
416	float f;
417	Reader->readData(*data_arr, boundid);
418	spacing->InsertNextValue (f = (*data_arr)[0]);
419	for (std::size_t j = 1; j < (2 * size); j++) // << Super laid
420	if (f != (*data_arr)[j])
421	spacing->InsertNextValue (f = (*data_arr)[j]);
422	}
423	else
424	{ // CellData, variable de coordonnÃ©e sans attributes bounds.
425	Reader->readData(*data_arr, coordinate);
426	spacing->InsertNextValue ((data_arr)[0] - ((data_arr)[1] - (*data_arr)[0])/2);
427	for (std::size_t j = 0; j < (size-1); j++)
428	spacing->InsertNextValue ((data_arr)[j] + ((data_arr)[j+1]-(*data_arr)[j])/2);
429	spacing->InsertNextValue ((data_arr)[size-2] + ((data_arr)[size-1] - (*data_arr)[size-2])/2);
430	}
431	}
432	else
433	{ // PointData, variable de coordonnÃ©e.
434	Reader->readData(*data_arr, coordinate);
435	for (std::size_t j = 0; j < (size-1); j++)
436	spacing->InsertNextValue ((*data_arr)[j]);
437	}
438	}
439	else
440	{
441	if (bounds) // CellData, dimension de coordonnÃ©e.
442	for (float i = -0.5; i < ((float)size + 1.); i = i + 1.)
443	spacing->InsertNextValue (i);
444	else // PointData, dimension de coordonnÃ©e.
445	for (float j = 0.; j < (float)size; j = j + 1.)
446	spacing->InsertNextValue (j);
447	}
448	}
449	//----------------------------------------------------------------
450
451	void vtkLSCEReader::AddPoint(vtkPoints * points, float * value, bool proj)
452	{
453	if (proj)
454	{
455	double h = 1.0;
456	double cartesianCoord[3];
457	value[0] = 0.017453292519943295 * value[0];
458	value[1] = 0.017453292519943295 * value[1];
459
460	cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
461	cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
462	cartesianCoord[2] = h * sin(value[1]);
463	points->InsertNextPoint (cartesianCoord);
464	}
465	else
466	{
467	points->InsertNextPoint (value);
468	}
469	}
470
471	//----------------------------------------------------------------
472
473	void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
474	const vtkStdString & ycoordinate,
475	const vtkStdString & zcoordinate,
476	vtkIntArray * dimensions, bool bounds)
477	{
478	int a = (bounds) ? 1 : 0;
479	dimensions->InsertNextValue
480	((Reader->getDimensions(&xcoordinate).begin())->second + a);
481	dimensions->InsertNextValue
482	((Reader->getDimensions(&ycoordinate).begin())->second + a);
483	if (zcoordinate.size() > 0)
484	{
485	dimensions->InsertNextValue
486	((Reader->getDimensions(&zcoordinate).begin())->second + a);
487	}
488	else
489	{
490	dimensions->InsertNextValue(1);
491	}
492	}
493
494	//----------------------------------------------------------------
495
496	void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
497	const vtkStdString & zcoordinate,
498	vtkIntArray * dimensions, bool bounds)
499	{
500	int a = (bounds) ? 1 : 0;
501	dimensions->InsertNextValue
502	(Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)] + a);
503	dimensions->InsertNextValue
504	(Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)] + a);
505	if (zcoordinate.size() > 0)
506	{
507	dimensions->InsertNextValue
508	((Reader->getDimensions(&zcoordinate).begin())->second + a);
509	}
510	else
511	{
512	dimensions->InsertNextValue(1);
513	}
514	}
515
516	//----------------------------------------------------------------
517	void vtkLSCEReader::GetPoints(const vtkStdString & xcoordinate,
518	const vtkStdString & ycoordinate,
519	const vtkStdString & zcoordinate,
520	bool bounds, bool proj,
521	vtkPoints * points, vtkIntArray * dimensions)
522	{
523	float value[3];
524	bool is2D = true;
525	ARRAY_CREATE(xvalue, float, 1, [1]);
526	ARRAY_CREATE(yvalue, float, 1, [1]);
527	vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New();
528	if (zcoordinate.size() > 0)
529	{
530	GetSpacings(zcoordinate, bounds, zspacing);
531	is2D = false;
532	}
533	else
534	zspacing->InsertNextValue(0);
535
536	dimensions->InsertNextValue
537	(Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)]);
538	dimensions->InsertNextValue
539	(Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)]);
540	dimensions->InsertNextValue(zspacing->GetNumberOfTuples());
541
542	if (!bounds)
543	{
544	Reader->readData(*xvalue, xcoordinate);
545	Reader->readData(*yvalue, ycoordinate);
546
547	for (int i = 0; i < dimensions->GetValue(0); i++)
548	for (int j = 0; j < dimensions->GetValue(1); j++)
549	for (int k = 0; k < dimensions->GetValue(2); k++)
550	{
551	value[0] = (xvalue)[idimensions->GetValue(1) + j];
552	value[1] = (yvalue)[idimensions->GetValue(1) + j];
553	value[2] = zspacing->GetValue(k);
554	AddPoint(points, value, proj);
555	}
556	return;
557	}
558
559	dimensions->SetValue(0, dimensions->GetValue(0) + 1);
560	dimensions->SetValue(1, dimensions->GetValue(1) + 1);
561
562	if (Reader->hasBounds(xcoordinate))
563	{
564	vtkStdString boundid = Reader->getBoundsId(xcoordinate);
565	Reader->readData(*xvalue, boundid);
566	}
567	if (Reader->hasBounds(ycoordinate))
568	{
569	vtkStdString boundid = Reader->getBoundsId(ycoordinate);
570	Reader->readData(*yvalue, boundid);
571	}
572	int ydim = dimensions->GetValue(0) - 1 ,
573	xdim = dimensions->GetValue(1) - 1;
574
575	for (int i = 0; i <= xdim; i++)
576	for (int j = 0; j <= ydim; j++)
577	for (int k = 0; k < dimensions->GetValue(2); k++)
578	{
579	if ((j != ydim) && (i != xdim))
580	{
581	value[0] = (xvalue)[4 (i * ydim + j) + 0];
582	value[1] = (yvalue)[4 (i * ydim + j) + 0];
583	value[2] = zspacing->GetValue(k);
584	AddPoint(points, value, proj);
585
586	continue;
587	}
588	else if ((j == ydim) && (i != xdim))
589	{
590	value[0] = (xvalue)[4 (i * ydim + j - 1) + 1];
591	value[1] = (yvalue)[4 (i * ydim + j - 1) + 1];
592	value[2] = zspacing->GetValue(k);
593	AddPoint(points, value, proj);
594
595	continue;
596	}
597	else if ((j != ydim) && (i == xdim))
598	{
599	value[0] = (xvalue)[4 ((i-1) * ydim + j) + 3];
600	value[1] = (yvalue)[4 ((i-1) * ydim + j) + 3];
601	value[2] = zspacing->GetValue(k);
602	AddPoint(points, value, proj);
603
604	continue;
605	}
606	else if ((j == ydim) && (i == xdim))
607	{
608	value[0] = (xvalue)[4 ((i-1) * ydim + j - 1) + 2];
609	value[1] = (yvalue)[4 ((i-1) * ydim + j - 1) + 2];
610	value[2] = zspacing->GetValue(k);
611	AddPoint(points, value, proj);
612
613	continue;
614	}
615	}
616	}
617
618	//----------------------------------------------------------------
619
620	void vtkLSCEReader::GetCellsAndPoints(const vtkStdString & xcoordinate,
621	const vtkStdString & ycoordinate,
622	const vtkStdString & zcoordinate,
623	bool bounds, bool proj, bool clean, std::size_t nbvertex,
624	vtkCellArray * cells, vtkPoints * points,
625	vtkIntArray * dimensions)
626	{
627	float value[3];
628	double h = 1.0;
629	double cartesianCoord[3], vbounds[6], range[2];
630	vtkIdType pts[nbvertex] ;
631	bool is2D = true;
632	ARRAY_CREATE(xvalue, float, 1, [1]);
633	ARRAY_CREATE(yvalue, float, 1, [1]);
634	vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New();
635	vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New();
636
637	if (zcoordinate.size() > 0)
638	{
639	GetSpacings(zcoordinate, bounds, zspacing);
640	is2D = false;
641	}
642	else
643	zspacing->InsertNextValue(0);
644
645	zspacing->GetRange (range);
646
647	std::size_t size = (*Reader->getDimensions(&xcoordinate).begin()).second;
648	if (Reader->hasBounds(xcoordinate))
649	{
650	vtkStdString boundid = Reader->getBoundsId(xcoordinate);
651	Reader->readData(*xvalue, boundid);
652	if (!proj)
653	{
654	vbounds[0] = *std::min_element(xvalue->begin(), xvalue->end());
655	vbounds[1] = *std::max_element(xvalue->begin(), xvalue->end());
656	}
657	else
658	{
659	vbounds[0] = -h;
660	vbounds[1] = h;
661	}
662	}
663	if (Reader->hasBounds(ycoordinate))
664	{
665	vtkStdString boundid = Reader->getBoundsId(ycoordinate);
666	Reader->readData(*yvalue, boundid);
667	if (!proj)
668	{
669	vbounds[2] = *std::min_element(yvalue->begin(), yvalue->end());
670	vbounds[3] = *std::max_element(yvalue->begin(), yvalue->end());
671	}
672	else
673	{
674	vbounds[2] = -h;
675	vbounds[3] = h;
676	}
677	}
678	if (!proj)
679	{
680	vbounds[4] = range[0];
681	vbounds[5] = range[1];
682	}
683	else
684	{
685	vbounds[4] = -h;
686	vbounds[5] = h;
687	}
688	locator->InitPointInsertion(points, vbounds, size * zspacing->GetNumberOfTuples() * nbvertex);
689
690	for (std::size_t u = 0; u < size; u++)
691	{
692	for (int v = 0; v < zspacing->GetNumberOfTuples(); v++)
693	{
694	for (std::size_t w = 0; w < nbvertex; w++)
695	{
696	value[0] = (xvalue)[nbvertex u + w];
697	value[1] = (yvalue)[nbvertex u + w];
698	value[2] = zspacing->GetValue(v);
699	pts[w] = nbvertex * u + w;
700
701	if (proj)
702	{
703	value[0] = 0.017453292519943295 * value[0];
704	value[1] = 0.017453292519943295 * value[1];
705
706	cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
707	cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
708	cartesianCoord[2] = h * sin(value[1]);
709	if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
710	else points->InsertPoint (pts[w], cartesianCoord);
711	}
712	else
713	{
714	cartesianCoord[0] = value[0];
715	cartesianCoord[1] = value[1];
716	cartesianCoord[2] = value[2];
717	if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
718	else points->InsertPoint (pts[w], cartesianCoord);
719	}
720	}
721	cells->InsertNextCell(nbvertex, pts);
722	}
723	}
724	}
725
726	//----------------------------------------------------------------
727
728	void vtkLSCEReader::AddScalarData(vtkDataSet * output,
729	const vtkStdString & varname,
730	std::size_t record, bool bounds)
731	{
732	vtkSmartPointer<vtkDoubleArray> data = vtkDoubleArray::New();
733	data->SetName(varname.c_str());
734	ARRAY_CREATE(data_arr, double, 1, [1]);
735	Reader->readData(*data_arr, varname, record);
736	Reader->replaceMissingValue(varname, *data_arr, vtkMath::Nan());
737
738
739	data->SetNumberOfValues (data_arr->size());
740	for (std::size_t i = 0; i < data_arr->size(); i++)
741	data->SetValue (i, (*data_arr)[i]);
742
743	if (bounds)
744	{
745	output->GetCellData()->AddArray(data);
746	output->GetCellData()->SetActiveScalars (varname.c_str());
747	}
748	else
749	{
750	output->GetPointData()->AddArray(data);
751	output->GetPointData()->SetActiveScalars (varname.c_str());
752	}
753	}
754
755	//----------------------------------------------------------------
756
757	int vtkLSCEReader::RequestData
758	(vtkInformation * vtkNotUsed(request),
759	vtkInformationVector ** vtkNotUsed(inputVector),
760	vtkInformationVector * outputVector)
761	{
762	vtkInformation * outInfo = outputVector->GetInformationObject(0);
763	vtkDataSet * output = vtkDataSet::GetData(outInfo);
764
765	if (!output)
766	{
767	switch (this->CurGridType)
768	{
769	case(RECTILINEAR) :
770	output = vtkRectilinearGrid::New();
771	break;
772	case(CURVILINEAR) :
773	output = vtkStructuredGrid::New();
774	break;
775	default :
776	output = vtkUnstructuredGrid::New();
777	break;
778	}
779	output->SetPipelineInformation(outInfo);
780	output->Delete();
781	}
782
783	vtkUnstructuredGrid * uoutput = vtkUnstructuredGrid::SafeDownCast(output);
784	vtkRectilinearGrid * routput = vtkRectilinearGrid::SafeDownCast(output);
785	vtkStructuredGrid * soutput = vtkStructuredGrid::SafeDownCast(output);
786
787	double time = 0.0;
788	if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()))
789	time = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())[0];
790
791	if (routput && !this->VarNames.empty())
792	{ // Grille rectiliniÃ©aire
793	vtkSmartPointer<vtkFloatArray>
794	xspacing = vtkSmartPointer<vtkFloatArray>::New(),
795	yspacing = vtkSmartPointer<vtkFloatArray>::New(),
796	zspacing = vtkSmartPointer<vtkFloatArray>::New();
797	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
798
799	//CreateSimpleGrid (-1, 10, -1, 10, -1, 10, xspacing, yspacing, zspacing, dims);
800	GetSpacings(Reader->getLonCoordName(*VarNames.begin()), true, xspacing);
801	GetSpacings(Reader->getLatCoordName(*VarNames.begin()), true, yspacing);
802	GetSpacings(Reader->getVertCoordName(*VarNames.begin()), true, zspacing);
803	/std::cout << zspacing->GetNumberOfTuples() << std::endl;/
804
805	dims->InsertNextValue(xspacing->GetNumberOfTuples());
806	dims->InsertNextValue(yspacing->GetNumberOfTuples());
807	dims->InsertNextValue(zspacing->GetNumberOfTuples());
808
809	CreateRectilinearGrid(routput, outInfo, xspacing, yspacing, zspacing, dims);
810	AddScalarData(routput, *VarNames.begin(), 0, true);
811
812	this->AGridDef = true;
813	}
814
815	if (soutput && !this->VarNames.empty())
816	{ // Grille structurÃ©e
817
818	vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
819	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
820
821	//CreateSimpleGrid (-1, 10, -1, 10, -1, 10, pts, dims);
822	GetPoints(Reader->getLonCoordName(*VarNames.begin()),
823	Reader->getLatCoordName(*VarNames.begin()),
824	Reader->getVertCoordName(*VarNames.begin()), true, true, pts, dims);
825
826	CreateStructuredGrid(soutput, outInfo, pts, dims);
827	AddScalarData(soutput, *VarNames.begin(), 0, true);
828
829	this->AGridDef = true;
830	}
831
832	if (uoutput && !this->VarNames.empty())
833	{ // Grille non-structurÃ©e
834	vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
835	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
836	vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
837	// VTK_QUAD, VTK_HEXAHEDRON, VTK_POLYGON
838	//CreateSimpleGrid (0, 12, 0, 12, 0, 12, pts, cells, dims);
839	//std::cout << Reader->getVertCoordName(*VarNames.begin()) << std::endl;
840
841	GetCellsAndPoints(Reader->getLonCoordName(*VarNames.begin()),
842	Reader->getLatCoordName(*VarNames.begin()),
843	Reader->getVertCoordName(*VarNames.begin()), true, true, true,
844	Reader->getNbVertex(*VarNames.begin()), cells, pts, dims);
845
846	CreateUnstructuredGrid(uoutput, outInfo, pts, cells, VTK_POLYGON);
847
848	AddScalarData(uoutput, VarNames.begin(), this->GetRecord(VarNames.begin(), time), true);
849
850	this->AGridDef = true;
851	}
852	return (1);
853	}
854
855	//----------------------------------------------------------------
856
857	void vtkLSCEReader::GetTimeInformation
858	(const vtkStdString & variable, vtkDoubleArray * values, double * timeRange)
859	{
860	std::string timename = this->Reader->getTimeCoordName(variable);
861	if (this->Reader->hasVariable(timename))
862	{
863	ARRAY_CREATE(data_arr, double, 1, [1]);
864	Reader->readData(*data_arr, timename);
865
866	values->SetNumberOfValues (data_arr->size());
867	for (std::size_t i = 0; i < data_arr->size(); i++)
868	values->SetValue (i, (*data_arr)[i]);
869	}
870	else
871	{
872	std::size_t nbtimestep = this->Reader->getNbOfTimestep();
873	values->SetNumberOfValues (nbtimestep);
874	for (std::size_t i = 0; i < nbtimestep; i++)
875	values->SetValue (i, i);
876	}
877	timeRange[0] = values->GetValue(0);
878	timeRange[1] = values->GetValue(values->GetNumberOfTuples()-1);
879	}
880
881	std::size_t vtkLSCEReader::GetRecord(const vtkStdString & variable, const double & value)
882	{
883	std::string timename = this->Reader->getTimeCoordName(variable);
884	if (this->Reader->hasVariable(timename))
885	{
886	ARRAY_CREATE(data_arr, double, 1, [1]);
887	Reader->readData(*data_arr, timename);
888	return (std::find(data_arr->begin(), data_arr->end(), value) - data_arr->begin());
889	}
890	else
891	{
892	return (value);
893	}
894	}
895
896	//----------------------------------------------------------------
897
898	int vtkLSCEReader::RequestInformation
899	(vtkInformation * vtkNotUsed(request),
900	vtkInformationVector ** vtkNotUsed(inputVector),
901	vtkInformationVector * outputVector)
902	{
903	vtkInformation * outInfo = outputVector->GetInformationObject(0);
904	vtkDataSet * output = vtkDataSet::GetData(outInfo);
905
906	this->AGridDef = false;
907
908	if (!this->VarNames.empty())
909	{
910	switch (this->CurGridType)
911	{
912	case(RECTILINEAR) :
913	output = vtkRectilinearGrid::New();
914	break;
915	case(CURVILINEAR) :
916	output = vtkStructuredGrid::New();
917	break;
918	default :
919	output = vtkUnstructuredGrid::New();
920	break;
921	}
922	output->SetPipelineInformation(outInfo);
923	output->Delete();
924	}
925
926	vtkRectilinearGrid * routput = vtkRectilinearGrid::SafeDownCast(output);
927	vtkStructuredGrid * soutput = vtkStructuredGrid::SafeDownCast(output);
928
929	if (!this->VarNames.empty())
930	{ // Informations temporelles
931	if (this->Reader->hasTemporalDim() && this->Reader->isTemporal(*VarNames.begin()))
932	{
933	vtkSmartPointer<vtkDoubleArray> timeStepValues = vtkSmartPointer<vtkDoubleArray>::New();
934	double timeRange[2];
935	this->GetTimeInformation(*VarNames.begin(), timeStepValues, timeRange);
936
937	outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
938	timeStepValues->GetPointer(0),
939	timeStepValues->GetNumberOfTuples());
940	outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
941	}
942	else
943	{
944	outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
945	outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
946	}
947
948	}
949
950	if (routput && !this->VarNames.empty())
951	{ // Grille rectiliniÃ©aire
952	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
953	this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
954	Reader->getLatCoordName(*VarNames.begin()),
955	Reader->getVertCoordName(*VarNames.begin()),
956	dims, true);
957
958	int extent[6] = { 0, dims->GetValue(0) - 1
959	, 0, dims->GetValue(1) - 1
960	, 0, dims->GetValue(2) - 1};
961	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
962	}
963
964	if (soutput && !this->VarNames.empty())
965	{ // Grille structurÃ©e
966
967	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
968	this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
969	Reader->getVertCoordName(*VarNames.begin()),
970	dims, true);
971
972	int extent[6] = { 0, dims->GetValue(0) - 1
973	, 0, dims->GetValue(1) - 1
974	, 0, dims->GetValue(2) - 1};
975	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
976	}
977
978	return (1);
979	}
980
981	///---------------------------------------------------------------
982
983	#ifndef LSCE_EXPORTS
984	} // namespace vtk
985	} // namespace xmlioserver
986	#endif //LSCE_EXPORTS

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: XMLIO_V2/dev/dev_rv/src4/xmlio/vtk/lscereader.cpp @ 241

Download in other formats: