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
53namespace 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)[i*dimensions->GetValue(1) + j];
543                     value[1] = (*yvalue)[i*dimensions->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.