XMLIOSERVER 0.4
Serveur d'Entrées/Sorties parallèles
lscereader.cpp
Aller à la documentation de ce fichier.
00001 
00002 #ifndef LSCE_EXPORTS
00003 #       include "lscereader.hpp"
00004 #else
00005 #       include "vtkLSCEReader.h"
00006 #endif //LSCE_EXPORTS
00007 
00008 
00009 // Bibliothèque VTK
00010 #include "vtkObjectFactory.h"
00011 #include "vtkInformation.h"
00012 #include "vtkInformationVector.h"
00013 #include "vtkStreamingDemandDrivenPipeline.h"
00014 
00015 #include "vtkPointLocator.h"
00016 
00017 #include <vtkProperty.h>
00018 #include <vtkDataSetMapper.h>
00019 #include <vtkActor.h>
00020 #include <vtkRenderWindow.h>
00021 #include <vtkRenderer.h>
00022 #include <vtkRenderWindowInteractor.h>
00023 #include <vtkLegendScaleActor.h>
00024 
00025 #include "vtkMath.h"
00026 
00027 #include "inetcdf4_impl.hpp"
00028 #include "inetcdf4_adv_impl.hpp"
00029 
00030 // Bibliothèque C++
00031 #include <algorithm>
00032 
00033 // Type Fortran
00034 #define ARRAY(valuetype, numdims) boost::shared_ptr<boost::multi_array<valuetype, numdims> >
00035 
00036 #define ARRAY_ASSIGN(value, valuetype, numdims, extent)\
00037    value.reset(new boost::multi_array<valuetype, numdims>(boost::extents extent)) 
00038 
00039 #define ARRAY_CREATE(value, valuetype, numdims, extent)\
00040    ARRAY(valuetype, numdims) value =                   \
00041    ARRAY(valuetype, numdims)(new boost::multi_array<valuetype, numdims>(boost::extents extent))
00042 
00043 // Type C
00044 #define ARRAY_C_ASSIGN(value, valuetype, numdims, extent)\
00045    value = ARRAY(valuetype, numdims)                     \
00046    (new boost::multi_array<valuetype, numdims>(boost::extents extent, c_storage_order()))
00047 
00048 #define ARRAY_C_CREATE(value, valuetype, numdims, extent)\
00049    ARRAY_C_ASSIGN(ARRAY(valuetype, numdims) value, valuetype, numdims, extent)
00050 
00051 
00052 #ifndef LSCE_EXPORTS
00053 namespace xmlioserver
00054 {
00055    namespace vtk
00056    {
00057 #endif //LSCE_EXPORTS
00058        
00060       vtkStandardNewMacro(vtkLSCEReader);
00061 
00062       vtkLSCEReader::vtkLSCEReader(void)
00063          : FileName(), CurGridType(RECTILINEAR), VarNames()
00064          , A3D(true), ATemporal(true), ACell(true), AGridDef(false)
00065          , Reader()
00066       {
00067          this->SetNumberOfInputPorts (0);
00068          this->SetNumberOfOutputPorts(1);
00069       }
00070 
00071       vtkLSCEReader::~vtkLSCEReader(void)
00072       {
00073 
00074       }
00075 
00077 
00078       void vtkLSCEReader::SetFileName(const vtkStdString & fileName)
00079       {
00080          this->FileName = fileName;
00081          this->Reader   = boost::shared_ptr<xmlioserver::io::CINetCDF4Adv>
00082                                        (new xmlioserver::io::CINetCDF4Adv(fileName));
00083          this->Modified();
00084       }
00085       
00086       void vtkLSCEReader::SetVariable(const vtkStdString & variable)
00087       {
00088          this->RemoveAllSelectedVariables();
00089          if (Reader->hasVariable(variable))
00090             this->AddVariableToSelection(variable);
00091          this->Modified();
00092       }
00093 
00094       const char * vtkLSCEReader::GetFileName(void) const
00095       {
00096          return (this->FileName.c_str());
00097       }
00098 
00099       void vtkLSCEReader::PrintSelf(ostream& os, vtkIndent indent)
00100       {
00101         this->Superclass::PrintSelf(os,indent);
00102       }
00103 
00104       //----------------------------------------------------------------
00105       void vtkLSCEReader::AddVariableToSelection(const vtkStdString & varName)
00106       {
00107          if (this->Reader->isRectilinear(varName))
00108             this->SetGridType (RECTILINEAR);
00109          else if (this->Reader->isCurvilinear(varName))
00110             this->SetGridType (CURVILINEAR);
00111          else
00112             this->SetGridType (UNSTRUCTURED);
00113 
00114          this->VarNames.insert(varName);
00115       }
00116 
00117       void vtkLSCEReader::RemoveSelectedVariable(const vtkStdString & varName)
00118       {
00119          std::set<vtkStdString>::iterator it = this->VarNames.find(varName);
00120          if (it != this->VarNames.end()) this->VarNames.erase(it);
00121       }
00122 
00123       void vtkLSCEReader::RemoveAllSelectedVariables(void)
00124       {
00125          this->VarNames.clear();
00126       }
00127 
00128       const std::set<vtkStdString> &
00129             vtkLSCEReader::GetSelectedVariables(void) const
00130       {
00131          return (this->VarNames);
00132       }
00133 
00134       bool vtkLSCEReader::HasSelectedVariable(void) const
00135       {
00136          return (this->VarNames.size() != 0);
00137       }
00138 
00139       //----------------------------------------------------------------
00140 
00141       void vtkLSCEReader::SetGridType (GridType type)
00142       {
00143          if ((this->IsUnstructured() && (type == CURVILINEAR)) ||
00144              (this->IsUnstructured() && (type == RECTILINEAR)) ||
00145              (this->IsCurvilinear()  && (type == RECTILINEAR)))
00146             this->RemoveAllSelectedVariables();
00147          this->CurGridType = type;
00148       }
00149 
00150       bool vtkLSCEReader::IsUnstructured(void) const
00151       {
00152          return (this->CurGridType == UNSTRUCTURED);
00153       }
00154 
00155       bool vtkLSCEReader::IsCurvilinear(void) const
00156       {
00157          return (this->CurGridType == CURVILINEAR);
00158       }
00159 
00160       bool vtkLSCEReader::IsRectilinear(void) const
00161       {
00162          return (this->CurGridType == RECTILINEAR);
00163       }
00164 
00165       //----------------------------------------------------------------
00166 
00167       void vtkLSCEReader::AcceptTemporalOnly(bool value)
00168       {
00169          this->ATemporal = value;
00170       }
00171 
00172       void vtkLSCEReader::Accept3DOnly(bool value)
00173       {
00174          this->A3D = value;
00175       }
00176 
00177       void vtkLSCEReader::AcceptCellOnly(bool value)
00178       {
00179          this->ACell = value;
00180       }
00181 
00182       //----------------------------------------------------------------
00183 
00184       int vtkLSCEReader::ProcessRequest
00185                            (vtkInformation       *  request,
00186                             vtkInformationVector ** vtkNotUsed(inputVector),
00187                             vtkInformationVector *  outputVector)
00188       {
00189          try
00190          {
00191           
00192             if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
00193             {
00194                return (this->RequestData(NULL, NULL, outputVector));
00195             }
00196 
00197             if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_OBJECT()))
00198             {
00199             //   return (this->RequestDataObject(NULL, NULL, outputVector));
00200             }
00201 
00202             if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
00203             {
00204                  return (this->RequestInformation(NULL, NULL, outputVector));
00205             }
00206 
00207             if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
00208             {
00209                return (this->RequestUpdateExtent(NULL, NULL, outputVector));
00210             }
00211          }
00212          catch (xmlioserver::CException & _exception)
00213          {
00214             std::cerr << _exception.getMessage() << std::endl;
00215             return (-1);
00216          }
00217          return (1);
00218       }
00219 
00220       //----------------------------------------------------------------
00221 
00222       void vtkLSCEReader::ShowVariable(const vtkStdString & filename,
00223                                        const vtkStdString & varname)
00224       {
00225          vtkSmartPointer<vtkLSCEReader> lscereader =
00226                vtkSmartPointer<vtkLSCEReader>::New();
00227 
00228          lscereader->SetFileName(filename);
00229          lscereader->AddVariableToSelection(varname);
00230          lscereader->Update();
00231 
00232          // Création d'un mappeur et d'un acteur.
00233          vtkSmartPointer<vtkDataSetMapper> mapper = vtkSmartPointer<vtkDataSetMapper>::New();
00234          mapper->SetInputConnection(lscereader->GetOutput()->GetProducerPort());
00235          mapper->SetScalarRange  (1000, 10000);
00236 
00237          vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
00238          actor->SetMapper(mapper);
00239          actor->GetProperty()->SetRepresentationToWireframe();
00240 
00241          // Création d'une fenêtre de rendu.
00242          vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
00243          vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
00244          renderWindow->AddRenderer(renderer);
00245          vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
00246             vtkSmartPointer<vtkRenderWindowInteractor>::New();
00247          renderWindowInteractor->SetRenderWindow(renderWindow);
00248 
00249 
00250          vtkSmartPointer<vtkLegendScaleActor> legendScaleActor =
00251             vtkSmartPointer<vtkLegendScaleActor>::New();
00252 
00253          // Ajout de l'acteur à la scène.
00254          renderer->AddActor(actor);
00255          renderer->AddActor(legendScaleActor);
00256 
00257          renderer->GradientBackgroundOn();
00258          renderer->SetBackground(1,1,1);
00259          renderer->SetBackground2(0,0,0);
00260 
00261          // Rendu et intéraction
00262          renderWindow->Render();
00263          renderWindowInteractor->Start();
00264       }
00265 
00266       //----------------------------------------------------------------
00267 
00268       void vtkLSCEReader::CreateRectilinearGrid(vtkRectilinearGrid * grid,
00269                                                 vtkInformation     * outInfo,
00270                                                 vtkFloatArray      * xspacing,
00271                                                 vtkFloatArray      * yspacing,
00272                                                 vtkFloatArray      * zspacing,
00273                                                 vtkIntArray        * dimensions)
00274       {
00275             int extent[6] = { 0, dimensions->GetValue(0) - 1
00276                             , 0, dimensions->GetValue(1) - 1
00277                             , 0, dimensions->GetValue(2) - 1};
00278             outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
00279             grid->SetExtent (extent);
00280             grid->SetDimensions (dimensions->GetValue(0),
00281                                  dimensions->GetValue(1),
00282                                  dimensions->GetValue(2));
00283             grid->SetXCoordinates (xspacing);
00284             grid->SetYCoordinates (yspacing);
00285             grid->SetZCoordinates (zspacing);
00286       }
00287 
00288       void vtkLSCEReader::CreateStructuredGrid(vtkStructuredGrid * grid,
00289                                                vtkInformation    * outInfo,
00290                                                vtkPoints         * points,
00291                                                vtkIntArray       * dimensions)
00292       {
00293             int extent[6] = { 0, dimensions->GetValue(0) - 1
00294                             , 0, dimensions->GetValue(1) - 1
00295                             , 0, dimensions->GetValue(2) - 1};
00296             outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
00297             grid->SetExtent (extent);
00298             grid->SetDimensions (dimensions->GetValue(0),
00299                                  dimensions->GetValue(1),
00300                                  dimensions->GetValue(2));
00301             grid->SetPoints(points);
00302       }
00303 
00304       void vtkLSCEReader::CreateUnstructuredGrid(vtkUnstructuredGrid * grid,
00305                                                  vtkInformation      * outInfo,
00306                                                  vtkPoints           * points,
00307                                                  vtkCellArray        * cells,
00308                                                  int                   celltype)
00309       {
00310          grid->SetPoints(points);
00311          grid->SetCells(celltype, cells);
00312       }
00313 
00314       //----------------------------------------------------------------
00315 
00316       void vtkLSCEReader::GetRectilinearConnectivity(int dimZ, int dimY, int dimX, vtkCellArray * cells)
00317       {
00318          if (dimZ == 1 && dimY >  1 && dimX >  1) // 2D
00319          {
00320             vtkIdType pts[4] ;
00321 
00322             for (int i = 0, k = 0; i < (dimX - 1); i++)
00323                for (int j = 0; j < (dimY - 1); j++, k++)
00324             {
00325                pts[0] =  i * dimY + j;
00326                pts[1] =  i * dimY + j + 1;
00327                pts[2] =  (i+1) * dimY + j + 1;
00328                pts[3] =  (i+1) * dimY + j;
00329                cells->InsertNextCell(4,  pts);
00330             }
00331          }
00332          else if (dimZ > 1 && dimY >  1 && dimX >  1) // 3D
00333          {
00334             vtkIdType pts[8];
00335             for (int i = 0, l = 0; i < (dimX - 1); i++)
00336                for (int j = 0; j < (dimY - 1); j++)
00337                   for (int k = 0; k < (dimZ - 1); k++, l++)
00338             {
00339                pts[0] =  i * dimY * dimZ + j * dimZ + k;
00340                pts[1] =  i * dimY * dimZ + j * dimZ + k + 1;
00341                pts[2] =  (i+1) * dimY * dimZ + j * dimZ + k + 1;
00342                pts[3] =  (i+1) * dimY * dimZ + j * dimZ + k;
00343 
00344                pts[4] =  i * dimY * dimZ  + (j + 1) * dimZ + k ;
00345                pts[5] =  i * dimY * dimZ  + (j + 1) * dimZ + k +1;
00346                pts[6] =  (i+1) * dimY * dimZ + (j + 1) * dimZ + k + 1;
00347                pts[7] =  (i+1) * dimY * dimZ + (j + 1) * dimZ + k ;
00348                cells->InsertNextCell(8,  pts);
00349             }
00350          }
00351       }
00352 
00353       //----------------------------------------------------------------
00354 
00355       void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
00356                                             int yi, int yf,
00357                                             int zi, int zf,
00358                                             vtkFloatArray * xspacing,
00359                                             vtkFloatArray * yspacing,
00360                                             vtkFloatArray * zspacing,
00361                                             vtkIntArray   * dimensions)
00362       {
00363          dimensions->InsertNextValue (xf - xi + 1); // X
00364          dimensions->InsertNextValue (yf - yi + 1); // Y
00365          dimensions->InsertNextValue (zf - zi + 1); // Z
00366          for (int i = xi; i <= xf; i++)
00367             xspacing->InsertNextValue (i);
00368          for (int j = yi; j <= yf; j++)
00369             yspacing->InsertNextValue (j);
00370          for (int k = zi; k <= zf; k++)
00371             zspacing->InsertNextValue (k);
00372       }
00373 
00374       void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
00375                                             int yi, int yf,
00376                                             int zi, int zf,
00377                                             vtkPoints   * points,
00378                                             vtkIntArray * dimensions)
00379       {
00380          dimensions->InsertNextValue (xf - xi + 1); // X
00381          dimensions->InsertNextValue (yf - yi + 1); // Y
00382          dimensions->InsertNextValue (zf - zi + 1); // Z
00383          for (int i = xi; i <= xf; i++)
00384             for (int j = yi; j <= yf; j++)
00385                for (int k = zi; k <= zf; k++)
00386                   points->InsertNextPoint(i, j, k);
00387       }
00388 
00389       void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
00390                                             int yi, int yf,
00391                                             int zi, int zf,
00392                                             vtkPoints    * points,
00393                                             vtkCellArray * cells,
00394                                             vtkIntArray  * dimensions)
00395       {
00396          this->CreateSimpleGrid (xi, xf, yi, yf, zi, zf, points, dimensions);
00397          int dimX = dimensions->GetValue(0),
00398              dimY = dimensions->GetValue(1),
00399              dimZ = dimensions->GetValue(2);
00400          this->GetRectilinearConnectivity(dimZ, dimY, dimX, cells);
00401       }
00402 
00403       //----------------------------------------------------------------
00404       void vtkLSCEReader::GetSpacings(const vtkStdString & coordinate,
00405                                       bool bounds, vtkFloatArray * spacing)
00406       {
00407          std::size_t size = (Reader->getDimensions(&coordinate).begin())->second;
00408          ARRAY_CREATE(data_arr, float, 1, [1]);
00409          if (Reader->hasVariable(coordinate))
00410          {
00411             if (bounds)
00412             {
00413                if (Reader->hasBounds(coordinate))
00414                { // CellData, variable de coordonnée avec attributes bounds.
00415                   vtkStdString boundid = Reader->getBoundsId(coordinate);
00416                   float f;
00417                   Reader->readData(*data_arr, boundid);
00418                   spacing->InsertNextValue (f = (*data_arr)[0]);
00419                   for (std::size_t j = 1; j < (2 * size); j++) // << Super laid
00420                      if (f != (*data_arr)[j])
00421                         spacing->InsertNextValue (f = (*data_arr)[j]);
00422                }
00423                else
00424                { // CellData, variable de coordonnée sans attributes bounds.
00425                   Reader->readData(*data_arr, coordinate);
00426                   spacing->InsertNextValue ((*data_arr)[0] - ((*data_arr)[1] - (*data_arr)[0])/2);
00427                   for (std::size_t j = 0; j < (size-1); j++)
00428                      spacing->InsertNextValue ((*data_arr)[j] + ((*data_arr)[j+1]-(*data_arr)[j])/2);
00429                   spacing->InsertNextValue ((*data_arr)[size-2] + ((*data_arr)[size-1] - (*data_arr)[size-2])/2);
00430                }
00431             }
00432             else
00433             { // PointData, variable de coordonnée.
00434                Reader->readData(*data_arr, coordinate);
00435                for (std::size_t j = 0; j < (size-1); j++)
00436                   spacing->InsertNextValue ((*data_arr)[j]);
00437             }
00438          }
00439          else
00440          {
00441             if (bounds) // CellData, dimension de coordonnée.
00442                for (float i = -0.5; i < ((float)size + 1.); i = i + 1.)
00443                   spacing->InsertNextValue (i);
00444             else // PointData, dimension de coordonnée.
00445                for (float j = 0.; j < (float)size; j = j + 1.)
00446                   spacing->InsertNextValue (j);
00447          }
00448       }
00449       //----------------------------------------------------------------
00450 
00451       void vtkLSCEReader::AddPoint(vtkPoints * points, float *  value, bool proj)
00452       {
00453          if (proj)
00454          {
00455             double h = 1.0;
00456             double cartesianCoord[3];
00457             value[0] =  0.017453292519943295 * value[0];
00458             value[1] =  0.017453292519943295 * value[1];
00459 
00460             cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
00461             cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
00462             cartesianCoord[2] = h * sin(value[1]);
00463             points->InsertNextPoint (cartesianCoord);
00464          }
00465          else
00466          {
00467             points->InsertNextPoint (value);
00468          }
00469       }
00470 
00471       //----------------------------------------------------------------
00472       
00473       void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
00474                                         const vtkStdString & ycoordinate,
00475                                         const vtkStdString & zcoordinate,
00476                                         vtkIntArray * dimensions, bool bounds)
00477       {
00478          int a = (bounds) ? 1 : 0;
00479          dimensions->InsertNextValue
00480             ((Reader->getDimensions(&xcoordinate).begin())->second + a);
00481          dimensions->InsertNextValue
00482             ((Reader->getDimensions(&ycoordinate).begin())->second + a);
00483          if (zcoordinate.size() > 0)
00484          {
00485             dimensions->InsertNextValue
00486                ((Reader->getDimensions(&zcoordinate).begin())->second + a);
00487          }
00488          else
00489          {
00490             dimensions->InsertNextValue(1);
00491          }
00492       }
00493 
00494       //----------------------------------------------------------------
00495       
00496       void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
00497                                         const vtkStdString & zcoordinate,
00498                                         vtkIntArray * dimensions, bool bounds)
00499       {
00500          int a = (bounds) ? 1 : 0;
00501          dimensions->InsertNextValue
00502             (Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)] + a);
00503          dimensions->InsertNextValue
00504             (Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)] + a);
00505          if (zcoordinate.size() > 0)
00506          {
00507             dimensions->InsertNextValue
00508                ((Reader->getDimensions(&zcoordinate).begin())->second + a);
00509          }
00510          else
00511          {
00512             dimensions->InsertNextValue(1);
00513          }
00514       }
00515 
00516       //----------------------------------------------------------------
00517       void vtkLSCEReader::GetPoints(const vtkStdString & xcoordinate,
00518                                     const vtkStdString & ycoordinate,
00519                                     const vtkStdString & zcoordinate,
00520                                     bool bounds, bool proj,
00521                                     vtkPoints * points, vtkIntArray * dimensions)
00522       {
00523          float value[3];
00524          bool is2D = true;
00525          ARRAY_CREATE(xvalue, float, 1, [1]);
00526          ARRAY_CREATE(yvalue, float, 1, [1]);
00527          vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New();
00528          if (zcoordinate.size() > 0)
00529          {
00530             GetSpacings(zcoordinate, bounds, zspacing);
00531             is2D = false;
00532          }
00533          else
00534             zspacing->InsertNextValue(0);
00535 
00536          dimensions->InsertNextValue
00537             (Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)]);
00538          dimensions->InsertNextValue
00539             (Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)]);
00540          dimensions->InsertNextValue(zspacing->GetNumberOfTuples());
00541 
00542          if (!bounds)
00543          {
00544             Reader->readData(*xvalue, xcoordinate);
00545             Reader->readData(*yvalue, ycoordinate);
00546 
00547             for (int i = 0; i < dimensions->GetValue(0); i++)
00548                for (int j = 0; j < dimensions->GetValue(1); j++)
00549                   for (int k = 0; k < dimensions->GetValue(2); k++)
00550                   {
00551                      value[0] = (*xvalue)[i*dimensions->GetValue(1) + j];
00552                      value[1] = (*yvalue)[i*dimensions->GetValue(1) + j];
00553                      value[2] = zspacing->GetValue(k);
00554                      AddPoint(points, value, proj);
00555                   }
00556             return;
00557          }
00558 
00559          dimensions->SetValue(0, dimensions->GetValue(0) + 1);
00560          dimensions->SetValue(1, dimensions->GetValue(1) + 1);
00561 
00562          if (Reader->hasBounds(xcoordinate))
00563          {
00564             vtkStdString boundid = Reader->getBoundsId(xcoordinate);
00565             Reader->readData(*xvalue, boundid);
00566          }
00567          if (Reader->hasBounds(ycoordinate))
00568          {
00569             vtkStdString boundid = Reader->getBoundsId(ycoordinate);
00570             Reader->readData(*yvalue, boundid);
00571          }
00572          int ydim = dimensions->GetValue(0) - 1 ,
00573              xdim = dimensions->GetValue(1) - 1;
00574 
00575          for (int i = 0; i <= xdim; i++)
00576             for (int j = 0; j <= ydim; j++)
00577                for (int k = 0; k < dimensions->GetValue(2); k++)
00578                {
00579                   if ((j != ydim) && (i != xdim))
00580                   {
00581                      value[0] = (*xvalue)[4 * (i * ydim + j) + 0];
00582                      value[1] = (*yvalue)[4 * (i * ydim + j) + 0];
00583                      value[2] = zspacing->GetValue(k);
00584                      AddPoint(points, value, proj);
00585 
00586                      continue;
00587                   }
00588                   else if ((j == ydim) && (i != xdim))
00589                   {
00590                      value[0] = (*xvalue)[4 * (i * ydim + j - 1) + 1];
00591                      value[1] = (*yvalue)[4 * (i * ydim + j - 1) + 1];
00592                      value[2] = zspacing->GetValue(k);
00593                      AddPoint(points, value, proj);
00594 
00595                      continue;
00596                   }
00597                   else if ((j != ydim) && (i == xdim))
00598                   {
00599                      value[0] = (*xvalue)[4 * ((i-1) * ydim + j) + 3];
00600                      value[1] = (*yvalue)[4 * ((i-1) * ydim + j) + 3];
00601                      value[2] = zspacing->GetValue(k);
00602                      AddPoint(points, value, proj);
00603 
00604                      continue;
00605                   }
00606                   else if ((j == ydim) && (i == xdim))
00607                   {
00608                      value[0] = (*xvalue)[4 * ((i-1) * ydim + j - 1) + 2];
00609                      value[1] = (*yvalue)[4 * ((i-1) * ydim + j - 1) + 2];
00610                      value[2] = zspacing->GetValue(k);
00611                      AddPoint(points, value, proj);
00612 
00613                      continue;
00614                   }
00615                }
00616       }
00617 
00618       //----------------------------------------------------------------
00619 
00620       void vtkLSCEReader::GetCellsAndPoints(const vtkStdString & xcoordinate,
00621                                             const vtkStdString & ycoordinate,
00622                                             const vtkStdString & zcoordinate,
00623                                             bool bounds, bool proj, bool clean, std::size_t nbvertex,
00624                                             vtkCellArray * cells, vtkPoints * points,
00625                                             vtkIntArray * dimensions)
00626       {
00627          float  value[3];
00628          double h = 1.0;
00629          double cartesianCoord[3], vbounds[6], range[2];
00630          vtkIdType pts[nbvertex] ;
00631          bool is2D = true;
00632          ARRAY_CREATE(xvalue, float, 1, [1]);
00633          ARRAY_CREATE(yvalue, float, 1, [1]);
00634          vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New();
00635          vtkSmartPointer<vtkFloatArray> zspacing  = vtkSmartPointer<vtkFloatArray>::New();
00636          
00637          if (zcoordinate.size() > 0)
00638          {
00639             GetSpacings(zcoordinate, bounds, zspacing);
00640             is2D = false;
00641          }
00642          else
00643             zspacing->InsertNextValue(0);
00644          
00645          zspacing->GetRange (range);
00646 
00647          std::size_t size = (*Reader->getDimensions(&xcoordinate).begin()).second;
00648          if (Reader->hasBounds(xcoordinate))
00649          {
00650             vtkStdString boundid = Reader->getBoundsId(xcoordinate);
00651             Reader->readData(*xvalue, boundid);
00652             if (!proj)
00653             {
00654                vbounds[0] = *std::min_element(xvalue->begin(), xvalue->end());
00655                vbounds[1] = *std::max_element(xvalue->begin(), xvalue->end());
00656             }
00657             else
00658             {
00659                vbounds[0] = -h;
00660                vbounds[1] = h;
00661             }
00662          }
00663          if (Reader->hasBounds(ycoordinate))
00664          {
00665             vtkStdString boundid = Reader->getBoundsId(ycoordinate);
00666             Reader->readData(*yvalue, boundid);
00667             if (!proj)
00668             {
00669                vbounds[2] = *std::min_element(yvalue->begin(), yvalue->end());
00670                vbounds[3] = *std::max_element(yvalue->begin(), yvalue->end());
00671             }
00672             else
00673             {
00674                vbounds[2] = -h;
00675                vbounds[3] = h;
00676             }
00677          }
00678          if (!proj)
00679          {
00680             vbounds[4] = range[0];
00681             vbounds[5] = range[1];
00682          }
00683          else
00684          {
00685             vbounds[4] = -h;
00686             vbounds[5] = h;
00687          }
00688          locator->InitPointInsertion(points, vbounds, size * zspacing->GetNumberOfTuples() * nbvertex); 
00689          
00690          for (std::size_t u = 0; u < size; u++)
00691          {
00692             for (int v = 0; v < zspacing->GetNumberOfTuples(); v++)
00693             {
00694                for (std::size_t w = 0; w < nbvertex; w++)
00695                {
00696                   value[0] = (*xvalue)[nbvertex * u + w];
00697                   value[1] = (*yvalue)[nbvertex * u + w];
00698                   value[2] = zspacing->GetValue(v);
00699                   pts[w] = nbvertex * u + w;
00700 
00701                   if (proj)
00702                   {
00703                      value[0] = 0.017453292519943295 * value[0];
00704                      value[1] = 0.017453292519943295 * value[1];
00705 
00706                      cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
00707                      cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
00708                      cartesianCoord[2] = h * sin(value[1]);
00709                      if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
00710                      else points->InsertPoint (pts[w], cartesianCoord);
00711                   }
00712                   else
00713                   {
00714                      cartesianCoord[0] = value[0];
00715                      cartesianCoord[1] = value[1];
00716                      cartesianCoord[2] = value[2];
00717                      if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
00718                      else points->InsertPoint (pts[w], cartesianCoord);
00719                   }
00720                }
00721                cells->InsertNextCell(nbvertex,  pts);
00722             }
00723          }
00724       }
00725 
00726       //----------------------------------------------------------------
00727 
00728       void vtkLSCEReader::AddScalarData(vtkDataSet * output,
00729                                         const vtkStdString & varname,
00730                                         std::size_t record, bool bounds)
00731       {
00732          vtkSmartPointer<vtkDoubleArray> data = vtkDoubleArray::New();
00733          data->SetName(varname.c_str());
00734          ARRAY_CREATE(data_arr, double, 1, [1]);
00735          Reader->readData(*data_arr, varname, record);
00736          Reader->replaceMissingValue(varname, *data_arr, vtkMath::Nan());
00737          
00738          
00739          data->SetNumberOfValues (data_arr->size());
00740          for (std::size_t i = 0; i < data_arr->size(); i++)
00741             data->SetValue (i, (*data_arr)[i]);
00742 
00743          if (bounds)
00744          {
00745             output->GetCellData()->AddArray(data);
00746             output->GetCellData()->SetActiveScalars (varname.c_str());
00747          }
00748          else
00749          {
00750             output->GetPointData()->AddArray(data);
00751             output->GetPointData()->SetActiveScalars (varname.c_str());
00752          }
00753       }
00754 
00755       //----------------------------------------------------------------
00756       
00757       int vtkLSCEReader::RequestData
00758                            (vtkInformation       *  vtkNotUsed(request),
00759                             vtkInformationVector ** vtkNotUsed(inputVector),
00760                             vtkInformationVector *  outputVector)
00761       {
00762          vtkInformation *  outInfo = outputVector->GetInformationObject(0);
00763          vtkDataSet     *  output  = vtkDataSet::GetData(outInfo);
00764          
00765          if (!output)
00766          {
00767             switch (this->CurGridType)
00768             {
00769                case(RECTILINEAR)  :
00770                   output = vtkRectilinearGrid::New();
00771                   break;
00772                case(CURVILINEAR)  :
00773                   output = vtkStructuredGrid::New();
00774                   break;
00775                default :
00776                   output = vtkUnstructuredGrid::New();
00777                   break;
00778             }
00779             output->SetPipelineInformation(outInfo);
00780             output->Delete();
00781          }
00782 
00783          vtkUnstructuredGrid * uoutput  = vtkUnstructuredGrid::SafeDownCast(output);
00784          vtkRectilinearGrid  * routput  = vtkRectilinearGrid::SafeDownCast(output);
00785          vtkStructuredGrid   * soutput  = vtkStructuredGrid::SafeDownCast(output);
00786          
00787         double time = 0.0;
00788         if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()))
00789            time = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())[0];
00790         
00791          if (routput && !this->VarNames.empty())
00792          { // Grille rectiliniéaire
00793             vtkSmartPointer<vtkFloatArray>
00794                            xspacing = vtkSmartPointer<vtkFloatArray>::New(),
00795                            yspacing = vtkSmartPointer<vtkFloatArray>::New(),
00796                            zspacing = vtkSmartPointer<vtkFloatArray>::New();
00797             vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
00798 
00799             //CreateSimpleGrid (-1, 10, -1, 10, -1, 10, xspacing, yspacing, zspacing, dims);
00800             GetSpacings(Reader->getLonCoordName(*VarNames.begin()), true, xspacing);
00801             GetSpacings(Reader->getLatCoordName(*VarNames.begin()), true, yspacing);
00802             GetSpacings(Reader->getVertCoordName(*VarNames.begin()), true, zspacing);
00803             /*std::cout << zspacing->GetNumberOfTuples() << std::endl;*/
00804 
00805             dims->InsertNextValue(xspacing->GetNumberOfTuples());
00806             dims->InsertNextValue(yspacing->GetNumberOfTuples());
00807             dims->InsertNextValue(zspacing->GetNumberOfTuples());
00808             
00809             CreateRectilinearGrid(routput, outInfo, xspacing, yspacing, zspacing, dims);
00810             AddScalarData(routput, *VarNames.begin(), 0, true);
00811 
00812             this->AGridDef = true;
00813          }
00814 
00815          if (soutput && !this->VarNames.empty())
00816          { // Grille structurée
00817              
00818             vtkSmartPointer<vtkPoints>   pts  = vtkSmartPointer<vtkPoints>::New();
00819             vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
00820             
00821             //CreateSimpleGrid (-1, 10, -1, 10, -1, 10, pts, dims);
00822             GetPoints(Reader->getLonCoordName(*VarNames.begin()),
00823                       Reader->getLatCoordName(*VarNames.begin()),
00824                       Reader->getVertCoordName(*VarNames.begin()), true, true, pts, dims);
00825 
00826             CreateStructuredGrid(soutput, outInfo, pts, dims);
00827             AddScalarData(soutput, *VarNames.begin(), 0, true);
00828 
00829             this->AGridDef = true;
00830          }
00831 
00832          if (uoutput  && !this->VarNames.empty())
00833          { // Grille non-structurée
00834             vtkSmartPointer<vtkPoints>    pts   = vtkSmartPointer<vtkPoints>::New();
00835             vtkSmartPointer<vtkIntArray>  dims  = vtkSmartPointer<vtkIntArray>::New();
00836             vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
00837             // VTK_QUAD, VTK_HEXAHEDRON, VTK_POLYGON
00838             //CreateSimpleGrid (0, 12, 0, 12, 0, 12, pts, cells, dims);
00839             //std::cout << Reader->getVertCoordName(*VarNames.begin()) << std::endl;
00840 
00841             GetCellsAndPoints(Reader->getLonCoordName(*VarNames.begin()),
00842                               Reader->getLatCoordName(*VarNames.begin()),
00843                               Reader->getVertCoordName(*VarNames.begin()), true, true, true,
00844                               Reader->getNbVertex(*VarNames.begin()), cells, pts, dims);
00845 
00846             CreateUnstructuredGrid(uoutput, outInfo, pts, cells, VTK_POLYGON);
00847          
00848             AddScalarData(uoutput, *VarNames.begin(), this->GetRecord(*VarNames.begin(), time), true);
00849             
00850             this->AGridDef = true;
00851          }
00852          return (1);
00853       }
00854 
00855       //----------------------------------------------------------------
00856       
00857       void vtkLSCEReader::GetTimeInformation
00858                (const vtkStdString & variable, vtkDoubleArray * values, double * timeRange)
00859       {
00860          std::string timename = this->Reader->getTimeCoordName(variable);
00861          if (this->Reader->hasVariable(timename))
00862          {
00863             ARRAY_CREATE(data_arr, double, 1, [1]);
00864             Reader->readData(*data_arr, timename);
00865             
00866             values->SetNumberOfValues (data_arr->size());
00867             for (std::size_t i = 0; i < data_arr->size(); i++)
00868                values->SetValue (i, (*data_arr)[i]);
00869          }
00870          else
00871          {
00872             std::size_t nbtimestep = this->Reader->getNbOfTimestep();
00873             values->SetNumberOfValues (nbtimestep);
00874             for (std::size_t i = 0; i < nbtimestep; i++)
00875                values->SetValue (i, i);
00876          }
00877          timeRange[0] = values->GetValue(0);
00878          timeRange[1] = values->GetValue(values->GetNumberOfTuples()-1);
00879       }
00880       
00881       std::size_t vtkLSCEReader::GetRecord(const vtkStdString & variable, const double & value)
00882       {
00883          std::string timename = this->Reader->getTimeCoordName(variable);
00884          if (this->Reader->hasVariable(timename))
00885          {
00886             ARRAY_CREATE(data_arr, double, 1, [1]);
00887             Reader->readData(*data_arr, timename);
00888             return (std::find(data_arr->begin(), data_arr->end(), value) - data_arr->begin());
00889          }
00890          else
00891          {
00892             return (value);
00893          }
00894       }
00895 
00896       //----------------------------------------------------------------
00897       
00898       int vtkLSCEReader::RequestInformation
00899                            (vtkInformation       *  vtkNotUsed(request),
00900                             vtkInformationVector ** vtkNotUsed(inputVector),
00901                             vtkInformationVector *  outputVector)
00902       {
00903          vtkInformation *  outInfo = outputVector->GetInformationObject(0);
00904          vtkDataSet     *  output  = vtkDataSet::GetData(outInfo);
00905 
00906          this->AGridDef = false;
00907          
00908          if (!this->VarNames.empty())
00909          {
00910             switch (this->CurGridType)
00911             {
00912                case(RECTILINEAR)  :
00913                   output = vtkRectilinearGrid::New();
00914                   break;
00915                case(CURVILINEAR)  :
00916                   output = vtkStructuredGrid::New();
00917                   break;
00918                default :
00919                   output = vtkUnstructuredGrid::New();
00920                   break;
00921             }
00922             output->SetPipelineInformation(outInfo);
00923             output->Delete();
00924          }
00925 
00926          vtkRectilinearGrid  * routput  = vtkRectilinearGrid::SafeDownCast(output);
00927          vtkStructuredGrid   * soutput  = vtkStructuredGrid::SafeDownCast(output);
00928          
00929          if (!this->VarNames.empty())
00930          { // Informations temporelles
00931             if (this->Reader->hasTemporalDim() && this->Reader->isTemporal(*VarNames.begin()))
00932             {
00933                vtkSmartPointer<vtkDoubleArray> timeStepValues = vtkSmartPointer<vtkDoubleArray>::New();
00934                double timeRange[2];
00935                this->GetTimeInformation(*VarNames.begin(), timeStepValues, timeRange);
00936 
00937                outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
00938                             timeStepValues->GetPointer(0),
00939                             timeStepValues->GetNumberOfTuples());
00940                outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
00941             }
00942             else
00943             {
00944                outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
00945                outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
00946             }
00947             
00948          }
00949          
00950          if (routput && !this->VarNames.empty())
00951          { // Grille rectiliniéaire
00952             vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
00953             this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
00954                                 Reader->getLatCoordName(*VarNames.begin()),
00955                                 Reader->getVertCoordName(*VarNames.begin()),
00956                                 dims, true);
00957             
00958             int extent[6] = { 0, dims->GetValue(0) - 1
00959                             , 0, dims->GetValue(1) - 1
00960                             , 0, dims->GetValue(2) - 1};
00961             outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
00962          }
00963          
00964          if (soutput && !this->VarNames.empty())
00965          { // Grille structurée
00966             
00967             vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
00968             this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
00969                                 Reader->getVertCoordName(*VarNames.begin()),
00970                                 dims, true);
00971             
00972             int extent[6] = { 0, dims->GetValue(0) - 1
00973                             , 0, dims->GetValue(1) - 1
00974                             , 0, dims->GetValue(2) - 1};
00975             outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
00976          }
00977             
00978          return (1);
00979       }
00980 
00982 
00983 #ifndef LSCE_EXPORTS
00984    } // namespace vtk
00985 } // namespace xmlioserver
00986 #endif //LSCE_EXPORTS
 Tout Classes Espaces de nommage Fichiers Fonctions Variables Définition de type Énumérations Valeurs énumérées Amis Macros