_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
#include <vtkHorizontalAverage.h>
#include <vtkDataSet.h>
#include <vtkStructuredGrid.h>
#include <vtkStructuredPoints.h>
#include <vtkRectilinearGrid.h>
#include <vtkDataArray.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
#include <vtkFloatArray.h>
vtkStandardNewMacro(vtkHorizontalAverage);
vtkHorizontalAverage::vtkHorizontalAverage()
{
this -> Normalize = 0;
}
void vtkHorizontalAverage::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Normalize: "
<< (this->Normalize ? this->Normalize : 0) << endl;
}
void vtkHorizontalAverage::Execute()
{vtkDataSet *input = this -> GetInput();
vtkStructuredGrid *output = this ->GetOutput();
int DataSetType;
DataSetType = input -> GetDataObjectType();
int dims[3],i,j,k;
// pointer for casting
vtkRectilinearGrid *rect_input = vtkRectilinearGrid::New();
vtkStructuredGrid *structgrid_input = vtkStructuredGrid::New();
vtkPointData *iPtr,*oPtr;
double *inPtr, *outPtr;
vtkDataArray *scaldat;
vtkDataArray *new_scaldat;
vtkFloatArray *scalars = vtkFloatArray::New();
vtkFloatArray *new_scalars = vtkFloatArray::New();
int numPts = input -> GetNumberOfPoints();
output -> Allocate(numPts);
if((DataSetType == VTK_RECTILINEAR_GRID) && (rect_input = vtkRectilinearGrid::SafeDownCast(input)))
{
rect_input -> GetDimensions(dims);
//get the scalar attribute Data
iPtr = rect_input->GetPointData();
scaldat = iPtr -> GetScalars();
oPtr = output->GetPointData();
new_scaldat = oPtr -> GetScalars();
scalars = vtkFloatArray::SafeDownCast(scaldat);
new_scalars = vtkFloatArray::SafeDownCast(new_scaldat);
if(!scaldat && !new_scaldat) vtkErrorMacro("Downcast to vtkFloatArray failed");
// inPtr = rect_input->GetPoint(0);
// outPtr = output->GetPoint(0);
}
else if((DataSetType == VTK_STRUCTURED_GRID )&& (structgrid_input = vtkStructuredGrid::SafeDownCast(input)))
{vtkPointData *iPtr = structgrid_input->GetPointData();
vtkPointData *oPtr = output->GetPointData();
iPtr = structgrid_input->GetPointData();
scaldat = iPtr -> GetScalars();
oPtr = output->GetPointData();
new_scaldat = oPtr -> GetScalars();
scalars = vtkFloatArray::SafeDownCast(scaldat);
new_scalars = vtkFloatArray::SafeDownCast(new_scaldat);
if(!scaldat && !new_scaldat) vtkErrorMacro("Downcast to vtkFloatArray failed");
structgrid_input->GetDimensions(dims);
}
else vtkErrorMacro("unsupported Dataset Format");
double *horizontal = new double[dims[0]];
double *horizontal_sigma = new double[dims[0]];
int index;
// compute the horizontal average
for (i=0; i<dims[0]; i++) {
horizontal[i] = 0.0;
for (j=0; j<dims[1]; j++) {
for (k=0; k<dims[2]; k++) {
index = k*(dims[1]*dims[0]) + j*dims[0] + i;
horizontal[i] += scalars -> GetComponent(index, 1); //horizontal[i] += inPtr[index];
}
}
horizontal[i] = horizontal[i]/(dims[1]*dims[2]);
}
//subtract horizontal average
for (i=0; i<dims[0]; i++) {
for (j=0; j<dims[1]; j++) {
for (k=0; k<dims[2]; k++) {
index = k*(dims[1]*dims[0]) + j*dims[0] + i;
new_scalars -> InsertNextTuple ((scalars -> GetScalar(index)) - horizontal[i]);
}
}
}
if (this->GetNormalize()) {
// compute sigma
for (i=0; i<dims[0]; i++) {
horizontal_sigma[i] = 0.0;
for (j=0; j<dims[1]; j++) {
for (k=0; k<dims[2]; k++) {
index = k*(dims[1]*dims[0]) + j*dims[0] + i;
horizontal_sigma[i] += (scalars -> GetComponent(index, 1))*(GetComponent(index, 1));
}
}
horizontal_sigma[i] /= (dims[1]*dims[2]-1);
horizontal_sigma[i] = sqrt(horizontal_sigma[i]);
}
// now normalize
for (i=0; i<dims[0]; i++) {
for (j=0; j<dims[1]; j++) {
for (k=0; k<dims[2]; k++) {
index = k*(dims[1]*dims[0]) + j*dims[0] + i;
//outPtr[index] /= horizontal_sigma[i];
new_scalars -> InsertNextTuple((new_scalars -> GetScalar(index)) /= horizontal_sigma[i])
}
}
}
}
delete [] horizontal;
delete [] horizontal_sigma;
//now set output
output -> CopyStructure(input);
output -> SetScalars (new_scalars);
iPtr->Delete();
oPtr->Delete();
delete [] inPtr;
delete [] outPtr;
}
#ifndef __vtkHorizontalAverage_h
#define __vtkHorizontalAverage_h
#include <vtkDataSetToStructuredGridFilter.h>
class VTK_EXPORT vtkHorizontalAverage : public vtkDataSetToStructuredGridFilter
{ public:
static vtkHorizontalAverage *New();
// vtkStandardNewMacro(vtkHorizontalAverage);
vtkTypeMacro(vtkHorizontalAverage,vtkDataSetToStructuredGridFilter);
vtkSetMacro(Normalize,int);
vtkGetMacro(Normalize,int);
vtkBooleanMacro(Normalize, int);
void PrintSelf(ostream&, vtkIndent);
protected: vtkHorizontalAverage();
~vtkHorizontalAverage(){};
int Normalize;
void operator=(const vtkHorizontalAverage&) {};
void Execute();
};
#endif
_______________________________________________
ParaView mailing list
[email protected]
http://www.paraview.org/mailman/listinfo/paraview
