Context Navigation

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

Last change on this file since 239 was 239, checked in by hozdoba, 10 years ago

File size: 34.7 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)
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
190	if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
191	{
192	return (this->RequestData(NULL, NULL, outputVector));
193	}
194
195	if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_OBJECT()))
196	{
197	// return (this->RequestDataObject(NULL, NULL, outputVector));
198	}
199
200	if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
201	{
202	return (this->RequestInformation(NULL, NULL, outputVector));
203	}
204
205	if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
206	{
207	return (this->RequestUpdateExtent(NULL, NULL, outputVector));
208	}
209
210	return (1);
211	}
212
213	//----------------------------------------------------------------
214
215	void vtkLSCEReader::ShowVariable(const vtkStdString & filename,
216	const vtkStdString & varname)
217	{
218	vtkSmartPointer<vtkLSCEReader> lscereader =
219	vtkSmartPointer<vtkLSCEReader>::New();
220
221	lscereader->SetFileName(filename);
222	lscereader->AddVariableToSelection(varname);
223	lscereader->Update();
224
225	// CrÃ©ation d'un mappeur et d'un acteur.
226	vtkSmartPointer<vtkDataSetMapper> mapper = vtkSmartPointer<vtkDataSetMapper>::New();
227	mapper->SetInputConnection(lscereader->GetOutput()->GetProducerPort());
228	mapper->SetScalarRange (1000, 10000);
229
230	vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
231	actor->SetMapper(mapper);
232	actor->GetProperty()->SetRepresentationToWireframe();
233
234	// CrÃ©ation d'une fenÃªtre de rendu.
235	vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
236	vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
237	renderWindow->AddRenderer(renderer);
238	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
239	vtkSmartPointer<vtkRenderWindowInteractor>::New();
240	renderWindowInteractor->SetRenderWindow(renderWindow);
241
242
243	vtkSmartPointer<vtkLegendScaleActor> legendScaleActor =
244	vtkSmartPointer<vtkLegendScaleActor>::New();
245
246	// Ajout de l'acteur Ã la scÃšne.
247	renderer->AddActor(actor);
248	renderer->AddActor(legendScaleActor);
249
250	renderer->GradientBackgroundOn();
251	renderer->SetBackground(1,1,1);
252	renderer->SetBackground2(0,0,0);
253
254	// Rendu et intÃ©raction
255	renderWindow->Render();
256	renderWindowInteractor->Start();
257	}
258
259	//----------------------------------------------------------------
260
261	void vtkLSCEReader::CreateRectilinearGrid(vtkRectilinearGrid * grid,
262	vtkInformation * outInfo,
263	vtkFloatArray * xspacing,
264	vtkFloatArray * yspacing,
265	vtkFloatArray * zspacing,
266	vtkIntArray * dimensions)
267	{
268	int extent[6] = { 0, dimensions->GetValue(0) - 1
269	, 0, dimensions->GetValue(1) - 1
270	, 0, dimensions->GetValue(2) - 1};
271	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
272	grid->SetExtent (extent);
273	grid->SetDimensions (dimensions->GetValue(0),
274	dimensions->GetValue(1),
275	dimensions->GetValue(2));
276	grid->SetXCoordinates (xspacing);
277	grid->SetYCoordinates (yspacing);
278	grid->SetZCoordinates (zspacing);
279	}
280
281	void vtkLSCEReader::CreateStructuredGrid(vtkStructuredGrid * grid,
282	vtkInformation * outInfo,
283	vtkPoints * points,
284	vtkIntArray * dimensions)
285	{
286	int extent[6] = { 0, dimensions->GetValue(0) - 1
287	, 0, dimensions->GetValue(1) - 1
288	, 0, dimensions->GetValue(2) - 1};
289	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
290	grid->SetExtent (extent);
291	grid->SetDimensions (dimensions->GetValue(0),
292	dimensions->GetValue(1),
293	dimensions->GetValue(2));
294	grid->SetPoints(points);
295	}
296
297	void vtkLSCEReader::CreateUnstructuredGrid(vtkUnstructuredGrid * grid,
298	vtkInformation * outInfo,
299	vtkPoints * points,
300	vtkCellArray * cells,
301	int celltype)
302	{
303	grid->SetPoints(points);
304	grid->SetCells(celltype, cells);
305	}
306
307	//----------------------------------------------------------------
308
309	void vtkLSCEReader::GetRectilinearConnectivity(int dimZ, int dimY, int dimX, vtkCellArray * cells)
310	{
311	if (dimZ == 1 && dimY > 1 && dimX > 1) // 2D
312	{
313	vtkIdType pts[4] ;
314
315	for (int i = 0, k = 0; i < (dimX - 1); i++)
316	for (int j = 0; j < (dimY - 1); j++, k++)
317	{
318	pts[0] = i * dimY + j;
319	pts[1] = i * dimY + j + 1;
320	pts[2] = (i+1) * dimY + j + 1;
321	pts[3] = (i+1) * dimY + j;
322	cells->InsertNextCell(4, pts);
323	}
324	}
325	else if (dimZ > 1 && dimY > 1 && dimX > 1) // 3D
326	{
327	vtkIdType pts[8];
328	for (int i = 0, l = 0; i < (dimX - 1); i++)
329	for (int j = 0; j < (dimY - 1); j++)
330	for (int k = 0; k < (dimZ - 1); k++, l++)
331	{
332	pts[0] = i * dimY * dimZ + j * dimZ + k;
333	pts[1] = i * dimY * dimZ + j * dimZ + k + 1;
334	pts[2] = (i+1) * dimY * dimZ + j * dimZ + k + 1;
335	pts[3] = (i+1) * dimY * dimZ + j * dimZ + k;
336
337	pts[4] = i * dimY * dimZ + (j + 1) * dimZ + k ;
338	pts[5] = i * dimY * dimZ + (j + 1) * dimZ + k +1;
339	pts[6] = (i+1) * dimY * dimZ + (j + 1) * dimZ + k + 1;
340	pts[7] = (i+1) * dimY * dimZ + (j + 1) * dimZ + k ;
341	cells->InsertNextCell(8, pts);
342	}
343	}
344	}
345
346	//----------------------------------------------------------------
347
348	void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
349	int yi, int yf,
350	int zi, int zf,
351	vtkFloatArray * xspacing,
352	vtkFloatArray * yspacing,
353	vtkFloatArray * zspacing,
354	vtkIntArray * dimensions)
355	{
356	dimensions->InsertNextValue (xf - xi + 1); // X
357	dimensions->InsertNextValue (yf - yi + 1); // Y
358	dimensions->InsertNextValue (zf - zi + 1); // Z
359	for (int i = xi; i <= xf; i++)
360	xspacing->InsertNextValue (i);
361	for (int j = yi; j <= yf; j++)
362	yspacing->InsertNextValue (j);
363	for (int k = zi; k <= zf; k++)
364	zspacing->InsertNextValue (k);
365	}
366
367	void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
368	int yi, int yf,
369	int zi, int zf,
370	vtkPoints * points,
371	vtkIntArray * dimensions)
372	{
373	dimensions->InsertNextValue (xf - xi + 1); // X
374	dimensions->InsertNextValue (yf - yi + 1); // Y
375	dimensions->InsertNextValue (zf - zi + 1); // Z
376	for (int i = xi; i <= xf; i++)
377	for (int j = yi; j <= yf; j++)
378	for (int k = zi; k <= zf; k++)
379	points->InsertNextPoint(i, j, k);
380	}
381
382	void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
383	int yi, int yf,
384	int zi, int zf,
385	vtkPoints * points,
386	vtkCellArray * cells,
387	vtkIntArray * dimensions)
388	{
389	this->CreateSimpleGrid (xi, xf, yi, yf, zi, zf, points, dimensions);
390	int dimX = dimensions->GetValue(0),
391	dimY = dimensions->GetValue(1),
392	dimZ = dimensions->GetValue(2);
393	this->GetRectilinearConnectivity(dimZ, dimY, dimX, cells);
394	}
395
396	//----------------------------------------------------------------
397	void vtkLSCEReader::GetSpacings(const vtkStdString & coordinate,
398	bool bounds, vtkFloatArray * spacing)
399	{
400	std::size_t size = (Reader->getDimensions(&coordinate).begin())->second;
401	ARRAY_CREATE(data_arr, float, 1, [1]);
402	if (Reader->hasVariable(coordinate))
403	{
404	if (bounds)
405	{
406	if (Reader->hasBounds(coordinate))
407	{ // CellData, variable de coordonnÃ©e avec attributes bounds.
408	vtkStdString boundid = Reader->getBoundsId(coordinate);
409	float f;
410	Reader->readData(*data_arr, boundid);
411	spacing->InsertNextValue (f = (*data_arr)[0]);
412	for (std::size_t j = 1; j < (2 * size); j++) // << Super laid
413	if (f != (*data_arr)[j])
414	spacing->InsertNextValue (f = (*data_arr)[j]);
415	}
416	else
417	{ // CellData, variable de coordonnÃ©e sans attributes bounds.
418	Reader->readData(*data_arr, coordinate);
419	spacing->InsertNextValue ((data_arr)[0] - ((data_arr)[1] - (*data_arr)[0])/2);
420	for (std::size_t j = 0; j < (size-1); j++)
421	spacing->InsertNextValue ((data_arr)[j] + ((data_arr)[j+1]-(*data_arr)[j])/2);
422	spacing->InsertNextValue ((data_arr)[size-2] + ((data_arr)[size-1] - (*data_arr)[size-2])/2);
423	}
424	}
425	else
426	{ // PointData, variable de coordonnÃ©e.
427	Reader->readData(*data_arr, coordinate);
428	for (std::size_t j = 0; j < (size-1); j++)
429	spacing->InsertNextValue ((*data_arr)[j]);
430	}
431	}
432	else
433	{
434	if (bounds) // CellData, dimension de coordonnÃ©e.
435	for (float i = -0.5; i < ((float)size + 1.); i = i + 1.)
436	spacing->InsertNextValue (i);
437	else // PointData, dimension de coordonnÃ©e.
438	for (float j = 0.; j < (float)size; j = j + 1.)
439	spacing->InsertNextValue (j);
440	}
441	}
442	//----------------------------------------------------------------
443
444	void vtkLSCEReader::AddPoint(vtkPoints * points, float * value, bool proj)
445	{
446	if (proj)
447	{
448	double h = 1.0;
449	double cartesianCoord[3];
450	value[0] = 0.017453292519943295 * value[0];
451	value[1] = 0.017453292519943295 * value[1];
452
453	cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
454	cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
455	cartesianCoord[2] = h * sin(value[1]);
456	points->InsertNextPoint (cartesianCoord);
457	}
458	else
459	{
460	points->InsertNextPoint (value);
461	}
462	}
463
464	//----------------------------------------------------------------
465
466	void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
467	const vtkStdString & ycoordinate,
468	const vtkStdString & zcoordinate,
469	vtkIntArray * dimensions, bool bounds)
470	{
471	int a = (bounds) ? 1 : 0;
472	dimensions->InsertNextValue
473	((Reader->getDimensions(&xcoordinate).begin())->second + a);
474	dimensions->InsertNextValue
475	((Reader->getDimensions(&ycoordinate).begin())->second + a);
476	if (zcoordinate.size() > 0)
477	{
478	dimensions->InsertNextValue
479	((Reader->getDimensions(&zcoordinate).begin())->second + a);
480	}
481	else
482	{
483	dimensions->InsertNextValue(1);
484	}
485	}
486
487	void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
488	const vtkStdString & zcoordinate,
489	vtkIntArray * dimensions, bool bounds)
490	{
491	int a = (bounds) ? 1 : 0;
492	dimensions->InsertNextValue
493	(Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)] + a);
494	dimensions->InsertNextValue
495	(Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)] + a);
496	if (zcoordinate.size() > 0)
497	{
498	dimensions->InsertNextValue
499	((Reader->getDimensions(&zcoordinate).begin())->second + a);
500	}
501	else
502	{
503	dimensions->InsertNextValue(1);
504	}
505	}
506
507	//----------------------------------------------------------------
508	void vtkLSCEReader::GetPoints(const vtkStdString & xcoordinate,
509	const vtkStdString & ycoordinate,
510	const vtkStdString & zcoordinate,
511	bool bounds, bool proj,
512	vtkPoints * points, vtkIntArray * dimensions)
513	{
514	float value[3];
515	bool is2D = true;
516	ARRAY_CREATE(xvalue, float, 1, [1]);
517	ARRAY_CREATE(yvalue, float, 1, [1]);
518	vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New();
519	if (zcoordinate.size() > 0)
520	{
521	GetSpacings(zcoordinate, bounds, zspacing);
522	is2D = false;
523	}
524	else
525	zspacing->InsertNextValue(0);
526
527	dimensions->InsertNextValue
528	(Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)]);
529	dimensions->InsertNextValue
530	(Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)]);
531	dimensions->InsertNextValue(zspacing->GetNumberOfTuples());
532
533	if (!bounds)
534	{
535	Reader->readData(*xvalue, xcoordinate);
536	Reader->readData(*yvalue, ycoordinate);
537
538	for (int i = 0; i < dimensions->GetValue(0); i++)
539	for (int j = 0; j < dimensions->GetValue(1); j++)
540	for (int k = 0; k < dimensions->GetValue(2); k++)
541	{
542	value[0] = (xvalue)[idimensions->GetValue(1) + j];
543	value[1] = (yvalue)[idimensions->GetValue(1) + j];
544	value[2] = zspacing->GetValue(k);
545	AddPoint(points, value, proj);
546	}
547	return;
548	}
549
550	dimensions->SetValue(0, dimensions->GetValue(0)+1);
551	dimensions->SetValue(1, dimensions->GetValue(1)+1);
552
553	if (Reader->hasBounds(xcoordinate))
554	{
555	vtkStdString boundid = Reader->getBoundsId(xcoordinate);
556	Reader->readData(*xvalue, boundid);
557	}
558	if (Reader->hasBounds(ycoordinate))
559	{
560	vtkStdString boundid = Reader->getBoundsId(ycoordinate);
561	Reader->readData(*yvalue, boundid);
562	}
563	int ydim = dimensions->GetValue(0) - 1 ,
564	xdim = dimensions->GetValue(1) - 1;
565
566	for (int i = 0; i <= xdim; i++)
567	for (int j = 0; j <= ydim; j++)
568	for (int k = 0; k < dimensions->GetValue(2); k++)
569	{
570	if ((j != ydim) && (i != xdim))
571	{
572	value[0] = (xvalue)[4 (i * ydim + j) + 0];
573	value[1] = (yvalue)[4 (i * ydim + j) + 0];
574	value[2] = zspacing->GetValue(k);
575	AddPoint(points, value, proj);
576
577	continue;
578	}
579	else if ((j == ydim) && (i != xdim))
580	{
581	value[0] = (xvalue)[4 (i * ydim + j - 1) + 1];
582	value[1] = (yvalue)[4 (i * ydim + j - 1) + 1];
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-1) * ydim + j) + 3];
591	value[1] = (yvalue)[4 ((i-1) * ydim + j) + 3];
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 - 1) + 2];
600	value[1] = (yvalue)[4 ((i-1) * ydim + j - 1) + 2];
601	value[2] = zspacing->GetValue(k);
602	AddPoint(points, value, proj);
603
604	continue;
605	}
606	}
607	}
608
609	//----------------------------------------------------------------
610
611	void vtkLSCEReader::GetCellsAndPoints(const vtkStdString & xcoordinate,
612	const vtkStdString & ycoordinate,
613	const vtkStdString & zcoordinate,
614	bool bounds, bool proj, bool clean, std::size_t nbvertex,
615	vtkCellArray * cells, vtkPoints * points,
616	vtkIntArray * dimensions)
617	{
618	float value[3], range[2];
619	double h = 1.0;
620	double cartesianCoord[3], vbounds[6];
621	vtkIdType pts[nbvertex] ;
622	bool is2D = true;
623	ARRAY_CREATE(xvalue, float, 1, [1]);
624	ARRAY_CREATE(yvalue, float, 1, [1]);
625	vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New();
626	vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New();
627
628	if (zcoordinate.size() > 0)
629	{
630	GetSpacings(zcoordinate, bounds, zspacing);
631	is2D = false;
632	}
633	else
634	zspacing->InsertNextValue(0);
635
636	zspacing->GetValueRange (range);
637
638	std::size_t size = (*Reader->getDimensions(&xcoordinate).begin()).second;
639	if (Reader->hasBounds(xcoordinate))
640	{
641	vtkStdString boundid = Reader->getBoundsId(xcoordinate);
642	Reader->readData(*xvalue, boundid);
643	if (!proj)
644	{
645	vbounds[0] = *std::min_element(xvalue->begin(), xvalue->end());
646	vbounds[1] = *std::max_element(xvalue->begin(), xvalue->end());
647	}
648	else
649	{
650	vbounds[0] = -h;
651	vbounds[1] = h;
652	}
653	}
654	if (Reader->hasBounds(ycoordinate))
655	{
656	vtkStdString boundid = Reader->getBoundsId(ycoordinate);
657	Reader->readData(*yvalue, boundid);
658	if (!proj)
659	{
660	vbounds[2] = *std::min_element(yvalue->begin(), yvalue->end());
661	vbounds[3] = *std::max_element(yvalue->begin(), yvalue->end());
662	}
663	else
664	{
665	vbounds[2] = -h;
666	vbounds[3] = h;
667	}
668	}
669	if (!proj)
670	{
671	vbounds[4] = range[0];
672	vbounds[5] = range[1];
673	}
674	else
675	{
676	vbounds[4] = -h;
677	vbounds[5] = h;
678	}
679	locator->InitPointInsertion(points, vbounds, size * zspacing->GetNumberOfTuples() * nbvertex);
680
681	for (std::size_t u = 0; u < size; u++)
682	{
683	for (int v = 0; v < zspacing->GetNumberOfTuples(); v++)
684	{
685	for (std::size_t w = 0; w < nbvertex; w++)
686	{
687	value[0] = (xvalue)[nbvertex u + w];
688	value[1] = (yvalue)[nbvertex u + w];
689	value[2] = zspacing->GetValue(v);
690	pts[w] = nbvertex * u + w;
691
692	if (proj)
693	{
694	value[0] = 0.017453292519943295 * value[0];
695	value[1] = 0.017453292519943295 * value[1];
696
697	cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
698	cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
699	cartesianCoord[2] = h * sin(value[1]);
700	if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
701	else points->InsertPoint (pts[w], cartesianCoord);
702	}
703	else
704	{
705	cartesianCoord[0] = value[0];
706	cartesianCoord[1] = value[1];
707	cartesianCoord[2] = value[2];
708	if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
709	else points->InsertPoint (pts[w], cartesianCoord);
710	}
711	}
712	cells->InsertNextCell(nbvertex, pts);
713	}
714	}
715	}
716
717	//----------------------------------------------------------------
718
719	void vtkLSCEReader::AddScalarData(vtkDataSet * output,
720	const vtkStdString & varname,
721	std::size_t record, bool bounds)
722	{
723	vtkSmartPointer<vtkDoubleArray> data = vtkDoubleArray::New();
724	data->SetName(varname.c_str());
725	ARRAY_CREATE(data_arr, double, 1, [1]);
726	Reader->readData(*data_arr, varname);
727	Reader->replaceMissingValue(varname, *data_arr, vtkMath::Nan());
728
729
730	data->SetNumberOfValues (data_arr->size());
731	for (std::size_t i = 0; i < data_arr->size(); i++)
732	data->SetValue (i, (*data_arr)[i]);
733
734	if (bounds)
735	{
736	output->GetCellData()->AddArray(data);
737	output->GetCellData()->SetActiveScalars (varname.c_str());
738	}
739	else
740	{
741	output->GetPointData()->AddArray(data);
742	output->GetPointData()->SetActiveScalars (varname.c_str());
743	}
744	}
745
746	//----------------------------------------------------------------
747
748	int vtkLSCEReader::RequestData
749	(vtkInformation * vtkNotUsed(request),
750	vtkInformationVector ** vtkNotUsed(inputVector),
751	vtkInformationVector * outputVector)
752	{
753	vtkInformation * outInfo = outputVector->GetInformationObject(0);
754	vtkDataSet * output = vtkDataSet::GetData(outInfo);
755
756	if (!output)
757	{
758	switch (this->CurGridType)
759	{
760	case(RECTILINEAR) :
761	output = vtkRectilinearGrid::New();
762	break;
763	case(CURVILINEAR) :
764	output = vtkStructuredGrid::New();
765	break;
766	default :
767	output = vtkUnstructuredGrid::New();
768	break;
769	}
770	output->SetPipelineInformation(outInfo);
771	output->Delete();
772	}
773
774	vtkUnstructuredGrid * uoutput = vtkUnstructuredGrid::SafeDownCast(output);
775	vtkRectilinearGrid * routput = vtkRectilinearGrid::SafeDownCast(output);
776	vtkStructuredGrid * soutput = vtkStructuredGrid::SafeDownCast(output);
777
778	if (routput && !this->VarNames.empty())
779	{ // Grille rectiliniÃ©aire
780	vtkSmartPointer<vtkFloatArray>
781	xspacing = vtkSmartPointer<vtkFloatArray>::New(),
782	yspacing = vtkSmartPointer<vtkFloatArray>::New(),
783	zspacing = vtkSmartPointer<vtkFloatArray>::New();
784	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
785
786	//CreateSimpleGrid (-1, 10, -1, 10, -1, 10, xspacing, yspacing, zspacing, dims);
787	GetSpacings(Reader->getLonCoordName(*VarNames.begin()), true, xspacing);
788	GetSpacings(Reader->getLatCoordName(*VarNames.begin()), true, yspacing);
789	GetSpacings(Reader->getVertCoordName(*VarNames.begin()), true, zspacing);
790	/std::cout << zspacing->GetNumberOfTuples() << std::endl;/
791
792	dims->InsertNextValue(xspacing->GetNumberOfTuples());
793	dims->InsertNextValue(yspacing->GetNumberOfTuples());
794	dims->InsertNextValue(zspacing->GetNumberOfTuples());
795
796	AddScalarData(routput, *VarNames.begin(), 0, true);
797	CreateRectilinearGrid(routput, outInfo, xspacing, yspacing, zspacing, dims);
798	}
799
800	if (soutput && !this->VarNames.empty())
801	{ // Grille structurÃ©e
802
803	vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
804	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
805
806	//CreateSimpleGrid (-1, 10, -1, 10, -1, 10, pts, dims);
807	GetPoints(Reader->getLonCoordName(*VarNames.begin()),
808	Reader->getLatCoordName(*VarNames.begin()),
809	Reader->getVertCoordName(*VarNames.begin()), true, true, pts, dims);
810
811	AddScalarData(soutput, *VarNames.begin(), 0, true);
812	CreateStructuredGrid(soutput, outInfo, pts, dims);
813	}
814
815	if (uoutput && !this->VarNames.empty())
816	{ // Grille non-structurÃ©e
817
818	vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
819	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
820	vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
821	// VTK_QUAD, VTK_HEXAHEDRON, VTK_POLYGON
822	//CreateSimpleGrid (0, 12, 0, 12, 0, 12, pts, cells, dims);
823	//std::cout << Reader->getVertCoordName(*VarNames.begin()) << std::endl;
824
825	GetCellsAndPoints(Reader->getLonCoordName(*VarNames.begin()),
826	Reader->getLatCoordName(*VarNames.begin()),
827	Reader->getVertCoordName(*VarNames.begin()), true, false, true,
828	Reader->getNbVertex(*VarNames.begin()), cells, pts, dims);
829
830	AddScalarData(uoutput, *VarNames.begin(), 0, true);
831	CreateUnstructuredGrid(uoutput, outInfo, pts, cells, VTK_POLYGON);
832
833
834	}
835	return (1);
836	}
837
838	int vtkLSCEReader::RequestInformation
839	(vtkInformation * vtkNotUsed(request),
840	vtkInformationVector ** vtkNotUsed(inputVector),
841	vtkInformationVector * outputVector)
842	{
843	vtkInformation * outInfo = outputVector->GetInformationObject(0);
844	vtkDataSet * output = vtkDataSet::GetData(outInfo);
845
846	if (true)
847	{
848	switch (this->CurGridType)
849	{
850	case(RECTILINEAR) :
851	output = vtkRectilinearGrid::New();
852	break;
853	case(CURVILINEAR) :
854	output = vtkStructuredGrid::New();
855	break;
856	default :
857	output = vtkUnstructuredGrid::New();
858	break;
859	}
860	output->SetPipelineInformation(outInfo);
861	output->Delete();
862	}
863
864	vtkRectilinearGrid * routput = vtkRectilinearGrid::SafeDownCast(output);
865	vtkStructuredGrid * soutput = vtkStructuredGrid::SafeDownCast(output);
866
867	if (routput && !this->VarNames.empty())
868	{ // Grille rectiliniÃ©aire
869	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
870	this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
871	Reader->getLatCoordName(*VarNames.begin()),
872	Reader->getVertCoordName(*VarNames.begin()),
873	dims, true);
874
875	int extent[6] = { 0, dims->GetValue(0) - 1
876	, 0, dims->GetValue(1) - 1
877	, 0, dims->GetValue(2) - 1};
878	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
879	}
880
881	if (soutput && !this->VarNames.empty())
882	{ // Grille structurÃ©e
883
884	vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
885	this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
886	Reader->getVertCoordName(*VarNames.begin()),
887	dims, true);
888
889	int extent[6] = { 0, dims->GetValue(0) - 1
890	, 0, dims->GetValue(1) - 1
891	, 0, dims->GetValue(2) - 1};
892	outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
893	}
894
895	return (1);
896	}
897
898	///---------------------------------------------------------------
899
900	#ifndef LSCE_EXPORTS
901	} // namespace vtk
902	} // namespace xmlioserver
903	#endif //LSCE_EXPORTS

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

Download in other formats: