vtkDSPFilterGroup.cxx 17.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkDSPFilterGroup.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkDSPFilterGroup.h"

#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
#include "vtkFloatArray.h"
#include "vtkIdList.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkStructuredGrid.h"
#include "vtkObjectFactory.h"
#include "vtkDSPFilterDefinition.h"

Sean McBride's avatar
Sean McBride committed
35
#include <cctype>
36
37
#include <vector>
#include <string>
38
39
40
41
42
43
44

vtkStandardNewMacro(vtkDSPFilterGroup);


class vtkDSPFilterGroupVectorIntSTLCloak
{
public:
45
  std::vector<int> m_vector;
46
47
48
49
};
class vtkDSPFilterGroupVectorVectorIntSTLCloak
{
public:
50
  std::vector< std::vector<int> > m_vector;
51
52
53
54
55
};

class vtkDSPFilterGroupVectorArraySTLCloak
{
public:
56
  std::vector<vtkFloatArray *> m_vector;
57
58
59
60
};
class vtkDSPFilterGroupVectorVectorArraySTLCloak
{
public:
61
  std::vector< std::vector<vtkFloatArray *> > m_vector;
62
63
64
65
};
class vtkDSPFilterGroupVectorStringSTLCloak
{
public:
66
  std::vector<std::string> m_vector;
67
68
69
70
71
};

class vtkDSPFilterGroupVectorDefinitionSTLCloak
{
public:
72
  std::vector<vtkDSPFilterDefinition *> m_vector;
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
};


//----------------------------------------------------------------------------
vtkDSPFilterGroup::vtkDSPFilterGroup()
{
  this->FilterDefinitions = new vtkDSPFilterGroupVectorDefinitionSTLCloak;
  this->CachedInputs = new vtkDSPFilterGroupVectorArraySTLCloak;
  this->CachedInputNames = new vtkDSPFilterGroupVectorStringSTLCloak;
  this->CachedInputTimesteps = new vtkDSPFilterGroupVectorIntSTLCloak;
  this->CachedOutputs = new vtkDSPFilterGroupVectorVectorArraySTLCloak;
  this->CachedOutputTimesteps = new vtkDSPFilterGroupVectorVectorIntSTLCloak;

  this->FilterDefinitions->m_vector.resize(0);
  this->CachedInputs->m_vector.resize(0);
  this->CachedInputNames->m_vector.resize(0);
  this->CachedInputTimesteps->m_vector.resize(0);
  this->CachedOutputs->m_vector.resize(0);
  this->CachedOutputTimesteps->m_vector.resize(0);
}

//----------------------------------------------------------------------------
vtkDSPFilterGroup::~vtkDSPFilterGroup()
{
  this->FilterDefinitions->m_vector.resize(0);
  this->CachedInputs->m_vector.resize(0);
  this->CachedInputNames->m_vector.resize(0);
  this->CachedInputTimesteps->m_vector.resize(0);
  this->CachedOutputs->m_vector.resize(0);
  this->CachedOutputTimesteps->m_vector.resize(0);

  delete this->FilterDefinitions;
  delete this->CachedInputs;
  delete this->CachedInputNames;
  delete this->CachedInputTimesteps;
  delete this->CachedOutputs;
  delete this->CachedOutputTimesteps;
}


//----------------------------------------------------------------------------
void vtkDSPFilterGroup::AddFilter(vtkDSPFilterDefinition *filter)
{
Kyle Lutz's avatar
Kyle Lutz committed
116
  //XXX can't just add this filter, need to check for duplicates and removals?
117
118
119
120
121
122
123

  vtkDSPFilterDefinition *thefilter = vtkDSPFilterDefinition::New();
  thefilter->Copy(filter);


  this->FilterDefinitions->m_vector.push_back( thefilter );

124
  std::vector<vtkFloatArray *> l_cachedOutsForThisFilter;
125
126
127
  l_cachedOutsForThisFilter.resize(0);
  this->CachedOutputs->m_vector.push_back( l_cachedOutsForThisFilter );

128
  std::vector<int> l_cachedOutTimesForThisFilter;
129
130
131
  l_cachedOutTimesForThisFilter.resize(0);
  this->CachedOutputTimesteps->m_vector.push_back(l_cachedOutTimesForThisFilter);

132

133
134
135
#if 0
  printf("**********************FILTERS AFTER ADDING FILTER***********************\n");
  for(int i=0;i<this->GetNumFilters();i++)
136
  {
137
138
139
140
141
      vtkDSPFilterDefinition *filterfromlist = this->FilterDefinitions->m_vector[i];
      printf("vtkDSPFilterGroup::AddFilter %d of %d  input=%s output=%s nums=%d dens=%d forwardnums=%d    this=%p\n",
       i,this->GetNumFilters(),
       filterfromlist->GetInputVariableName(),
       filterfromlist->GetOutputVariableName(),
142
143
       filterfromlist->GetNumNumeratorWeights(),
       filterfromlist->GetNumDenominatorWeights(),
144
145
       filterfromlist->GetNumForwardNumeratorWeights(),
       this);
146
  }
147
148
149
150
151
  printf("************************************************************************\n");
#endif
}

//----------------------------------------------------------------------------
152
void vtkDSPFilterGroup::RemoveFilter(const char *a_outputVariableName)
153
{
154
155
156
  std::vector<vtkDSPFilterDefinition *>::iterator l_iter;
  std::vector< std::vector<vtkFloatArray *> >::iterator l_cachedOutputsIter = this->CachedOutputs->m_vector.begin();
  std::vector< std::vector<int> >::iterator l_cachedOutputTimesIter = this->CachedOutputTimesteps->m_vector.begin();
157

158
  for(l_iter=this->FilterDefinitions->m_vector.begin();l_iter!=this->FilterDefinitions->m_vector.end();++l_iter)
159
  {
160
161
      if(!strcmp(a_outputVariableName,(*l_iter)->GetOutputVariableName()))
      {
162
163
    //this is the filter to delete
    this->FilterDefinitions->m_vector.erase(l_iter);
164
    if(l_cachedOutputsIter!=this->CachedOutputs->m_vector.end())
165
      this->CachedOutputs->m_vector.erase(l_cachedOutputsIter);
166
    if(l_cachedOutputTimesIter!=this->CachedOutputTimesteps->m_vector.end())
167
168
      this->CachedOutputTimesteps->m_vector.erase(l_cachedOutputTimesIter);
    break;
169
      }
170
171
      ++l_cachedOutputsIter;
      ++l_cachedOutputTimesIter;
172
  }
173

174
175
176
#if 0
  printf("**********************FILTERS AFTER REMOVING FILTER*********************\n");
  for(int i=0;i<this->GetNumFilters();i++)
177
  {
178
179
180
181
182
      vtkDSPFilterDefinition *filterfromlist = this->FilterDefinitions[i];
      printf("vtkDSPFilterGroup::RemoveFilter %d of %d  input=%s output=%s nums=%d dens=%d    this=%p\n",
       i,this->GetNumFilters(),
       filterfromlist->GetInputVariableName(),
       filterfromlist->GetOutputVariableName(),
183
184
       filterfromlist->GetNumNumeratorWeights(),
       filterfromlist->GetNumDenominatorWeights(),
185
       this);
186
  }
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
  printf("************************************************************************\n");
#endif
}


//----------------------------------------------------------------------------
void vtkDSPFilterGroup::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
}


//----------------------------------------------------------------------------
const char *vtkDSPFilterGroup::GetInputVariableName( int a_whichFilter )
{
  return this->FilterDefinitions->m_vector[a_whichFilter]->GetInputVariableName();
}
//----------------------------------------------------------------------------
bool vtkDSPFilterGroup::IsThisInputVariableInstanceNeeded( const char *a_name, int a_timestep, int a_outputTimestep )
{
  for(int i=0; i<this->GetNumFilters(); i++)
  {
209
      if( !strcmp(this->FilterDefinitions->m_vector[i]->GetInputVariableName(),a_name) )
210
      {
211
212
    if( this->FilterDefinitions->m_vector[i]->IsThisInputVariableInstanceNeeded(a_timestep,a_outputTimestep) )
    {
213
        return(true);
214
    }
215
216
217
218
219
220
221
222
223
      }
  }
  return(false);
}
//----------------------------------------------------------------------------
bool vtkDSPFilterGroup::IsThisInputVariableInstanceCached( const char *a_name, int a_timestep )
{
  for(int i=0;i<(int)this->CachedInputTimesteps->m_vector.size();i++)
  {
224
      if(this->CachedInputTimesteps->m_vector[i]==a_timestep)
225
      {
226
227
    if( this->CachedInputNames->m_vector[i]==a_name )
    {
228
        return(true);
229
    }
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
      }
  }
  return(false);
}
//----------------------------------------------------------------------------
void vtkDSPFilterGroup::AddInputVariableInstance( const char *a_name, int a_timestep, vtkFloatArray *a_data )
{
  //This assumes that the instance is not already cached! perhaps should check anyway?

  this->CachedInputTimesteps->m_vector.push_back(a_timestep);
  this->CachedInputNames->m_vector.push_back(a_name);

  vtkFloatArray *l_array = vtkFloatArray::New();
  l_array->DeepCopy(a_data);
  this->CachedInputs->m_vector.push_back(l_array);
}

//----------------------------------------------------------------------------
vtkFloatArray *vtkDSPFilterGroup::GetCachedInput( int a_whichFilter, int a_whichTimestep )
{
250
  std::string l_inputName = this->FilterDefinitions->m_vector[a_whichFilter]->GetInputVariableName();
251
252
  for(int i=0;i<(int)this->CachedInputTimesteps->m_vector.size();i++)
  {
253
      if(this->CachedInputTimesteps->m_vector[i]==a_whichTimestep)
254
      {
255
256
    if( this->CachedInputNames->m_vector[i]==l_inputName )
    {
257
        return(this->CachedInputs->m_vector[i]);
258
    }
259
260
      }
  }
261
  return(nullptr);
262
263
264
265
266
267
268
269
}


//----------------------------------------------------------------------------
vtkFloatArray *vtkDSPFilterGroup::GetCachedOutput( int a_whichFilter, int a_whichTimestep )
{
  for(int i=0;i<(int)this->CachedOutputs->m_vector[a_whichFilter].size();i++)
  {
270
271
      if(a_whichTimestep==this->CachedOutputTimesteps->m_vector[a_whichFilter][i])
      {
272
273
274
    vtkFloatArray *l_tmp = (this->CachedOutputs->m_vector[a_whichFilter])[i];
    if( !strcmp( l_tmp->GetName(),
           this->FilterDefinitions->m_vector[a_whichFilter]->GetOutputVariableName() ))
275
    {
276
277
        //printf("vtkDSPFilterGroup::GetCachedOutput found time %d output in cache\n",a_whichTimestep);
        return(l_tmp);
278
    }
279
280
281

    //else printf("vtkDSPFilterGroup::GetCachedOutput DID NOT FIND time %d output in cache %s %s\n",a_whichTimestep,
    //        l_tmp->GetName(), this->FilterDefinitions[a_whichFilter]->OutputVariableName.c_str() );
282
      }
283
284
  }

285
  return(nullptr);
286
287
288
289
290
291
292
293
294
295
296
297
}


//----------------------------------------------------------------------------
void vtkDSPFilterGroup::Copy( vtkDSPFilterGroup *other )
{
  this->FilterDefinitions->m_vector = other->FilterDefinitions->m_vector;
}



//----------------------------------------------------------------------------
298
int vtkDSPFilterGroup::GetNumFilters( )
299
{
300
  return static_cast<int>(this->FilterDefinitions->m_vector.size());
301
302
303
304
305
306
307
308
309
310
311
}


//----------------------------------------------------------------------------
vtkDSPFilterDefinition * vtkDSPFilterGroup::GetFilter(int a_whichFilter)
{
  return this->FilterDefinitions->m_vector[a_whichFilter];
}


//----------------------------------------------------------------------------
312
vtkFloatArray *vtkDSPFilterGroup::GetOutput( int a_whichFilter, int a_whichTimestep, int &a_instancesCalculated )
313
314
315
316
317
318
319
{
  int i,j,k;
  int l_numFilters = this->GetNumFilters();



  if( (int)this->CachedOutputs->m_vector.size() < l_numFilters )
320
  {
Kyle Lutz's avatar
Kyle Lutz committed
321
      //this shouldn't happen with saf. Should happen 1 time with exodus.
322
323
324
325
      //printf("vtkDSPFilterGroup::GetOutput resizing cache vector\n");

      int l_numNow=(int)this->CachedOutputs->m_vector.size();
      for(i=l_numNow;i<l_numFilters;i++)
326
      {
327
    std::vector<vtkFloatArray *> l_cachedOutsForThisFilter;
328
329
330
    l_cachedOutsForThisFilter.resize(0);
    this->CachedOutputs->m_vector.push_back( l_cachedOutsForThisFilter );

331
    std::vector<int> l_cachedOutTimesForThisFilter;
332
333
    l_cachedOutTimesForThisFilter.resize(0);
    this->CachedOutputTimesteps->m_vector.push_back(l_cachedOutTimesForThisFilter);
334
      }
335
336
337
338
339
  }

  //is this output array already cached?
  vtkFloatArray *l_tmp = this->GetCachedOutput( a_whichFilter, a_whichTimestep );
  if(l_tmp)
340
  {
341
342
      //printf("vtkDSPFilterGroup::GetOutput found time %d output in cache\n",a_whichTimestep);
      return(l_tmp);
343
  }
344
345
346
347
348
349
350
351
352
353
  //else printf("vtkDSPFilterGroup::GetOutput DID NOT FIND time %d output in cache (%d cache slots)\n",
  //      a_whichTimestep,(int)this->CachedOutputs[a_whichFilter].size() );


  vtkFloatArray *l_output = vtkFloatArray::New();
  l_output->SetName( FilterDefinitions->m_vector[a_whichFilter]->GetOutputVariableName() );

  int l_numNumerators = (int)FilterDefinitions->m_vector[a_whichFilter]->GetNumNumeratorWeights();
  int l_numForwardNumerators = (int)FilterDefinitions->m_vector[a_whichFilter]->GetNumForwardNumeratorWeights();
  if(!l_numNumerators && !l_numForwardNumerators)
354
  {
355
      printf("vtkDSPFilterGroup::GetOutput there are no numerator filter weights?\n");
356
      return(nullptr);
357
  }
358
359
360
361
  int l_numDenominators = (int)FilterDefinitions->m_vector[a_whichFilter]->GetNumDenominatorWeights();

  double l_a1 = 1.0;
  if(l_numDenominators)
362
  {
363
      l_a1 = FilterDefinitions->m_vector[a_whichFilter]->GetDenominatorWeight(0);
364
  }
365
366
367
368
369
370
371
372
373
374


  //printf("vtkDSPFilterGroup::GetOutput numerators=%d forwardnums=%d dens=%d\n",
  //   l_numNumerators,l_numForwardNumerators,l_numDenominators);


  //There should always be a valid input at the same time as an output
  vtkFloatArray *l_firstInput = this->GetCachedInput(a_whichFilter,a_whichTimestep);

  if(!l_firstInput)
375
  {
376
      printf("\n  vtkDSPFilterGroup::GetOutput error time %d has no input\n\n",a_whichTimestep);
377
      return(nullptr);
378
  }
379
380
381
382
383

  const int l_numEntries = l_firstInput->GetNumberOfTuples();
  const int l_numComponents = l_firstInput->GetNumberOfComponents();

  if(!l_numEntries || !l_numComponents)
384
  {
385
386
      printf("\n  vtkDSPFilterGroup::GetOutput error time %d, l_numEntries=%d, l_numComponents=%d\n\n",
       a_whichTimestep,l_numEntries,l_numComponents);
387
      return(nullptr);
388
  }
389
390
391
392
393
394
395
396

  //printf("vtkDSPFilterGroup::GetOutput first input entries=%d comps=%d\n",l_numEntries,l_numComponents);


  l_output->SetNumberOfComponents(l_numComponents);
  l_output->SetNumberOfTuples(l_numEntries);

  for( i=0; i<l_numNumerators; i++ )
397
  {
398
399
400
401
402
403
404
405
406
407
408
      int l_useThisTimestep = a_whichTimestep-i;
      double l_weight = this->FilterDefinitions->m_vector[a_whichFilter]->GetNumeratorWeight(i)/l_a1;

      if(l_useThisTimestep < 0) l_useThisTimestep=0; //pre-time is considered infinite procession of input value at time 0

      //printf("vtkDSPFilterGroup::GetOutput numerator weight %d is %e (incl a1=%e) time=%d\n",i,l_weight,l_a1,l_useThisTimestep);

      vtkFloatArray *l_input = this->GetCachedInput(a_whichFilter,l_useThisTimestep);
      float *l_outPtr = (float *)l_output->GetVoidPointer(0);

      if(!i)
409
      {
410
    for(j=0;j<l_numEntries*l_numComponents;j++) l_outPtr[i]=0;
411
      }
412
413

      if(l_input)
414
      {
415
416
417
    float *l_inPtr = (float *)l_input->GetVoidPointer(0);
    for(j=0;j<l_numEntries;j++)
    {
418
419
        for(k=0;k<l_numComponents;k++)
        {
420
421
422
      l_outPtr[0] += l_weight * l_inPtr[0];
      l_inPtr++;
      l_outPtr++;
423
        }
424
425
426
    }
      }
      else
427
      {
Kyle Lutz's avatar
Kyle Lutz committed
428
    printf("error vtkDSPFilterGroup::GetOutput can't get input %d\n",l_useThisTimestep);
429
      }
430
431
432
433
434
  }



  for( i=1; i<l_numDenominators; i++ )
435
  {
436
437
438
439
      double l_weight = this->FilterDefinitions->m_vector[a_whichFilter]->GetDenominatorWeight(i)/l_a1;


      if(a_whichTimestep-i < 0) break;//pre-time outputs are considered to be zero
440

441
      //printf("vtkDSPFilterGroup::GetOutput denominator weight %d is %e (incl a1=%e) time=%d\n",i,l_weight,l_a1,a_whichTimestep-i);
442

443
444
445
446
447
      vtkFloatArray *l_input = this->GetOutput( a_whichFilter, a_whichTimestep-i, a_instancesCalculated );

      float *l_outPtr = (float *)l_output->GetVoidPointer(0);

      if(l_input)
448
      {
449
450
451
    float *l_inPtr = (float *)l_input->GetVoidPointer(0);
    for(j=0;j<l_numEntries;j++)
    {
452
453
        for(k=0;k<l_numComponents;k++)
        {
454
455
456
      l_outPtr[0] -= l_weight * l_inPtr[0];
      l_inPtr++;
      l_outPtr++;
457
        }
458
459
460
461
462
463
464
    }
      }
  }


  //Handle forward inputs
  for( i=0; i<l_numForwardNumerators; i++ )
465
  {
466
467
468
469
470
      int l_useThisTimestep = a_whichTimestep+i+1;
      double l_weight = this->FilterDefinitions->m_vector[a_whichFilter]->GetForwardNumeratorWeight(i)/l_a1;


      float *l_outPtr = (float *)l_output->GetVoidPointer(0);
471

472
473
474
      vtkFloatArray *l_input = this->GetCachedInput(a_whichFilter,l_useThisTimestep);

      while(!l_input && l_useThisTimestep>=0)
475
      {
476
477
478
479
480
    //printf("         time %d failed......trying prev time.....\n",l_useThisTimestep);

    //Try the timestep before: all post-time inputs are considered to be the same as the last input
    l_useThisTimestep--;
    l_input = this->GetCachedInput(a_whichFilter,l_useThisTimestep);
481
      }
482
483

      if(l_input)
484
      {
485
486
487
488
489
490

    //printf("vtkDSPFilterGroup::GetOutput forward numerator weight %d is %e (incl a1=%e) time=%d\n",i,l_weight,l_a1,l_useThisTimestep);

    float *l_inPtr = (float *)l_input->GetVoidPointer(0);
    for(j=0;j<l_numEntries;j++)
    {
491
492
        for(k=0;k<l_numComponents;k++)
        {
493
494
495
      l_outPtr[0] += l_weight * l_inPtr[0];
      l_inPtr++;
      l_outPtr++;
496
        }
497
498
499
    }
      }
      else
500
      {
Kyle Lutz's avatar
Kyle Lutz committed
501
    printf("\nerror vtkDSPFilterGroup::GetOutput can't get forward input %d\n\n",l_useThisTimestep);
502
      }
503
504
505
506
  }


#if 0 //debug print
507
  {
508
509
510
511
512
     float *l_outPtr = (float *)l_output->GetVoidPointer(0);
     float *l_inPtr = (float *)l_firstInput->GetVoidPointer(0);
     float l_maxDiff=0;
     for(j=0;j<l_numEntries;j++)
     {
513
514
   for(k=0;k<l_numComponents;k++)
   {
515
       if( l_inPtr[0] - l_outPtr[0] )
516
       {
517
518
     if( fabs(l_inPtr[0] - l_outPtr[0]) > l_maxDiff ) l_maxDiff = fabs(l_inPtr[0] - l_outPtr[0]);

519

520
521
522
523
     printf("j=%d k=%d \t in=%f \t out=%f \t diff=%e   maxdiff=%e   diffperc=%f\n",j,k,
               l_inPtr[0],l_outPtr[0],l_inPtr[0] - l_outPtr[0],l_maxDiff,
      fabs(l_inPtr[0] - l_outPtr[0]) / fabs(l_inPtr[0]) );

524
       }
525
526
       l_inPtr++;
       l_outPtr++;
527
   }
528
     }
529
  }
530
531
532
533
534

#endif


  a_instancesCalculated++;
535

536
537
538
  //printf("****vtkDSPFilterGroup::GetOutput calculated  filter=%d time=%d entries=%d comps=%d***    out cache was %d slots\n",a_whichFilter,
  // a_whichTimestep,l_numEntries,l_numComponents,
  // this->CachedOutputs[a_whichFilter].size()  );
539

540
541
542
543
544
545
546
547

  this->CachedOutputs->m_vector[a_whichFilter].push_back(l_output);
  this->CachedOutputTimesteps->m_vector[a_whichFilter].push_back(a_whichTimestep);


  return(l_output);
}