Commit 8abfcf02 authored by Mathieu Westphal's avatar Mathieu Westphal Committed by Kitware Robot

Merge topic 'SymmetricTensor6ComponentSupport2'

555019bc Adding support for 6 component symmetric tensor array
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Reviewed-by: Berk Geveci's avatarBerk Geveci <berk.geveci@kitware.com>
Merge-request: !2082
parents 0f3fc6b8 555019bc
......@@ -80,6 +80,7 @@ static int TestJacobiN();
static int TestClampValue();
static int TestClampValues();
static int TestClampAndNormalizeValue();
static int TestTensorFromSymmetricTensor();
static int TestGetScalarTypeFittingRange();
static int TestGetAdjustedScalarRange();
static int TestExtentIsWithinOtherExtent();
......@@ -155,6 +156,7 @@ int UnitTestMath(int,char *[])
status += TestClampValue();
status += TestClampValues();
status += TestClampAndNormalizeValue();
status += TestTensorFromSymmetricTensor();
status += TestGetScalarTypeFittingRange();
status += TestGetAdjustedScalarRange();
status += TestExtentIsWithinOtherExtent();
......@@ -3279,6 +3281,54 @@ int TestClampAndNormalizeValue()
return status;
}
// Validate by checking symmetric tensor values
// are in correct places
int TestTensorFromSymmetricTensor()
{
int status = 0;
std::cout << "TensorFromSymmetricTensor..";
double symmTensor[9];
for (int i = 0; i < 6; i++)
{
symmTensor[i] = vtkMath::Random();
}
double tensor[9];
vtkMath::TensorFromSymmetricTensor(symmTensor, tensor);
if (tensor[0] != symmTensor[0] ||
tensor[1] != symmTensor[3] ||
tensor[2] != symmTensor[5] ||
tensor[3] != symmTensor[3] ||
tensor[4] != symmTensor[1] ||
tensor[5] != symmTensor[4] ||
tensor[6] != symmTensor[5] ||
tensor[7] != symmTensor[4] ||
tensor[8] != symmTensor[2])
{
std::cout << " Unexpected results from TensorFromSymmetricTensor " << std::endl;
++status;
}
vtkMath::TensorFromSymmetricTensor(symmTensor);
for (int i = 0; i < 9; i++)
{
if (symmTensor[i] != tensor[i])
{
std::cout << " Unexpected results from in place TensorFromSymmetricTensor " << std::endl;
++status;
break;
}
}
if (status)
{
std::cout << "..FAILED" << std::endl;
}
else
{
std::cout << ".PASSED" << std::endl;
}
return status;
}
// Validate by checking ranges with numeric_limits
int TestGetScalarTypeFittingRange()
{
......
......@@ -49,7 +49,8 @@ template <class Scalar> void vtkAngularPeriodicDataArray<Scalar>
return;
}
if (data->GetNumberOfComponents() != 3 && data->GetNumberOfComponents() != 9)
if (data->GetNumberOfComponents() != 3 && data->GetNumberOfComponents() != 6 &&
data->GetNumberOfComponents() != 9)
{
vtkWarningMacro(<< "Original data has " << data->GetNumberOfComponents() <<
" components, Expecting 3 or 9.");
......@@ -130,13 +131,18 @@ Transform(Scalar* pos) const
vtkMath::Normalize(pos);
}
}
else if (this->NumberOfComponents == 9)
else if (this->NumberOfComponents == 9 || this->NumberOfComponents == 6)
{
// Template type force a copy to a double array for tensor
double localPos [9];
double tmpMat [9];
double tmpMat2 [9];
std::copy(pos, pos + 9, localPos);
double localPos[9];
double tmpMat[9];
double tmpMat2[9];
std::copy(pos, pos + this->NumberOfComponents, localPos);
if (this->NumberOfComponents == 6)
{
vtkMath::TensorFromSymmetricTensor(localPos);
}
vtkMatrix3x3::Transpose(this->RotationMatrix->GetData(), tmpMat);
vtkMatrix3x3::Multiply3x3(this->RotationMatrix->GetData(), localPos, tmpMat2);
vtkMatrix3x3::Multiply3x3(tmpMat2, tmpMat, localPos);
......
......@@ -910,7 +910,6 @@ public:
static int SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int xOrder,
double **mt);
/**
* Solves for the least squares best fit matrix for the equation X'M' = Y'.
* Uses pseudoinverse to get the ordinary least squares.
......@@ -1095,6 +1094,21 @@ public:
static double ClampAndNormalizeValue(double value,
const double range[2]);
/**
* Convert a 6-Component symmetric tensor into a 9-Component tensor, no allocation performed.
* Symmetric tensor is expected to have the following order : XX, YY, ZZ, XY, YZ, XZ
*/
template<class T1, class T2>
static void TensorFromSymmetricTensor(T1 symmTensor[6], T2 tensor[9]);
/**
* Convert a 6-Component symmetric tensor into a 9-Component tensor, overwriting
* the tensor input.
* Symmetric tensor is expected to have the following order : XX, YY, ZZ, XY, YZ, XZ
*/
template<class T>
static void TensorFromSymmetricTensor(T tensor[9]);
/**
* Return the scalar type that is most likely to have enough precision
* to store a given range of data once it has been scaled and shifted
......@@ -1178,7 +1192,6 @@ public:
* Test if a number has finite value i.e. it is normal, subnormal or zero, but not infinite or Nan.
*/
static bool IsFinite(double x);
protected:
vtkMath() {}
~vtkMath() VTK_OVERRIDE {}
......@@ -1506,6 +1519,32 @@ inline double vtkMath::ClampAndNormalizeValue(double value,
return result;
}
//-----------------------------------------------------------------------------
template<class T1, class T2>
inline void vtkMath::TensorFromSymmetricTensor(T1 symmTensor[9], T2 tensor[9])
{
for (int i = 0; i < 3; i++)
{
tensor[4*i] = symmTensor[i];
}
tensor[1] = tensor[3] = symmTensor[3];
tensor[2] = tensor[6] = symmTensor[5];
tensor[5] = tensor[7] = symmTensor[4];
}
//-----------------------------------------------------------------------------
template<class T>
inline void vtkMath::TensorFromSymmetricTensor(T tensor[9])
{
tensor[6] = tensor[5]; // XZ
tensor[7] = tensor[4]; // YZ
tensor[8] = tensor[2]; // ZZ
tensor[4] = tensor[1]; // YY
tensor[5] = tensor[7]; // YZ
tensor[2] = tensor[6]; // XZ
tensor[1] = tensor[3]; // XY
}
namespace vtk_detail
{
// Can't specialize templates inside a template class, so we move the impl here.
......
......@@ -1263,14 +1263,15 @@ int vtkDataSetAttributes::CheckNumberOfComponents(vtkAbstractArray* aa,
}
else if ( vtkDataSetAttributes::AttributeLimits[attributeType] == EXACT )
{
if ( numComp !=
vtkDataSetAttributes::NumberOfAttributeComponents[attributeType] )
if (numComp ==
vtkDataSetAttributes::NumberOfAttributeComponents[attributeType] ||
(numComp == 6 && attributeType == TENSORS)) // TENSORS6 support
{
return 0;
return 1;
}
else
{
return 1;
return 0;
}
}
else if ( vtkDataSetAttributes::AttributeLimits[attributeType] == NOLIMIT )
......@@ -1633,7 +1634,6 @@ void vtkDataSetAttributes::InternalCopyAllocate(
newDA->SetLookupTable(list.LUT[i]);
}
// If attribute data, do something extra
if ( i < NUM_ATTRIBUTES )
{
......
1173cc950f18506abe0c78029830e0c7
4232484d090d519c4bcc90f2c7f7d27d
......@@ -85,8 +85,14 @@ class TestTensorGlyph(Testing.vtkTest):
g4.Update()
g4.SetPosition((2.0, 2.0, 0.0))
# 6Components symetric tensor
g5 = SimpleGlyph(reader)
g5.glyph.SetInputArrayToProcess(0, 0, 0, 0, "symTensors1")
g5.SetPosition((4.0, 2.0, 0.0))
g5.Update()
ren = vtk.vtkRenderer()
for i in (g1, g2, g3, g4):
for i in (g1, g2, g3, g4, g5):
for j in i.GetActors():
ren.AddActor(j)
......@@ -95,7 +101,7 @@ class TestTensorGlyph(Testing.vtkTest):
cam = ren.GetActiveCamera()
cam.Azimuth(-20)
cam.Elevation(20)
cam.Zoom(1.5)
cam.Zoom(1.1)
ren.SetBackground(0.5, 0.5, 0.5)
......
......@@ -893,9 +893,16 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
int normalize[9])
{
int i, normalizeAny, updated=0;
int numComp = 9;
vtkDataArray *fieldArray[9];
for (i=0; i<9; i++)
// Check for symmetric tensor input
if (arrayComp[6] == -1 || arrays[6] == NULL)
{
numComp = 6;
}
for (i = 0; i < numComp; i++)
{
if ( arrays[i] == NULL )
{
......@@ -903,7 +910,7 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
}
}
for ( i=0; i < 9; i++ )
for ( i=0; i < numComp; i++ )
{
fieldArray[i] = this->GetFieldArray(fd, arrays[i], arrayComp[i]);
......@@ -914,7 +921,7 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
}
}
for (normalizeAny=i=0; i < 9; i++)
for (normalizeAny = i = 0; i < numComp; i++)
{
updated |= this->UpdateComponentRange(fieldArray[i], componentRange[i]);
if ( num != (componentRange[i][1] - componentRange[i][0] + 1) )
......@@ -926,7 +933,7 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
}
vtkDataArray *newTensors;
for (i=1; i < 9; i++) //see whether all the data is from the same array
for (i = 1; i < numComp; i++) //see whether all the data is from the same array
{
if ( fieldArray[i] != fieldArray[i-1] )
{
......@@ -935,7 +942,7 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
}
// see whether we can reuse the data array from the field
if ( i >= 9 && fieldArray[0]->GetNumberOfComponents() == 9 &&
if ( i >= numComp && fieldArray[0]->GetNumberOfComponents() == numComp &&
fieldArray[0]->GetNumberOfTuples() == num && !normalizeAny )
{
newTensors = fieldArray[0];
......@@ -943,11 +950,11 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
}
else //have to copy data into created array
{
newTensors = vtkDataArray::CreateDataArray(this->GetComponentsType(9, fieldArray));
newTensors->SetNumberOfComponents(9);
newTensors = vtkDataArray::CreateDataArray(this->GetComponentsType(numComp, fieldArray));
newTensors->SetNumberOfComponents(numComp);
newTensors->SetNumberOfTuples(num);
for ( i=0; i < 9; i++ )
for ( i=0; i < numComp; i++ )
{
if ( this->ConstructArray(newTensors, i, fieldArray[i], arrayComp[i],
componentRange[i][0], componentRange[i][1],
......@@ -963,7 +970,7 @@ void vtkFieldDataToAttributeDataFilter::ConstructTensors(int num, vtkFieldData *
newTensors->Delete();
if ( updated ) //reset for next execution pass
{
for (i=0; i < 9; i++)
for (i=0; i < numComp; i++)
{
componentRange[i][0] = componentRange[i][1] = -1;
}
......
......@@ -1405,15 +1405,17 @@ void vtkQuadricDecimation::ComputeNumberOfComponents(void)
// Tensors attributes
if (pd->GetTensors() != NULL && this->TensorsAttribute)
{
for (j = 0; j < 9; j++)
vtkDataArray* inTensors = pd->GetTensors();
int nComp = inTensors->GetNumberOfComponents();
for (j = 0; j < nComp; j++)
{
pd->GetTensors()->GetRange(range, j);
inTensors->GetRange(range, j);
maxRange = (maxRange < (range[1] - range[0]) ?
(range[1] - range[0]) : maxRange);
}
if (maxRange != 0.0)
{
this->NumberOfComponents += 9;
this->NumberOfComponents += nComp;
pd->CopyTensorsOn();
this->AttributeScale[4] = this->TensorsWeight/maxRange;
}
......
......@@ -269,8 +269,12 @@ int vtkTensorGlyph::RequestData(
ptIncr = numDirs * inPtId * numSourcePts;
// Translation is postponed
// Symmetric tensor support
inTensors->GetTuple(inPtId, tensor);
if (inTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
// compute orientation vectors and scale factors from tensor
if ( this->ExtractEigenvalues ) // extract appropriate eigenfunctions
......
......@@ -146,6 +146,10 @@ int vtkExtractTensorComponents::RequestData(
for (ptId=0; ptId < numPts; ptId++)
{
inTensors->GetTuple(ptId, tensor);
if (inTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
if ( this->ExtractScalars )
{
......
......@@ -437,6 +437,10 @@ int vtkHyperStreamline::RequestData(
for (k=0; k < cell->GetNumberOfPoints(); k++)
{
cellTensors->GetTuple(k, tensor);
if (cellTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
for (j=0; j<3; j++)
{
for (i=0; i<3; i++)
......@@ -522,6 +526,10 @@ int vtkHyperStreamline::RequestData(
for (k=0; k < cell->GetNumberOfPoints(); k++)
{
cellTensors->GetTuple(k, tensor);
if (cellTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
for (j=0; j<3; j++)
{
for (i=0; i<3; i++)
......@@ -602,6 +610,10 @@ int vtkHyperStreamline::RequestData(
for (k=0; k < cell->GetNumberOfPoints(); k++)
{
cellTensors->GetTuple(k, tensor);
if (cellTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
for (j=0; j<3; j++)
{
for (i=0; i<3; i++)
......
......@@ -334,8 +334,10 @@ void vtkAngularPeriodicFilter::ComputeAngularPeriodicData(
int attribute = data->IsArrayAnAttribute(i);
vtkDataArray* array = data->GetArray(i);
vtkDataArray* transformedArray;
// Perdiodic copy of vector (3 components) or tensor (9 components) data
if (array->GetNumberOfComponents() == 3 || array->GetNumberOfComponents() == 9)
// Perdiodic copy of vector (3 components) or symmectric tensor (6 component, converted to 9 )
// or tensor (9 components) data
int numComp = array->GetNumberOfComponents();
if (numComp == 3 || numComp == 6 || numComp == 9)
{
transformedArray = this->TransformDataArray(array, angle, false,
attribute == vtkDataSetAttributes::NORMALS);
......
......@@ -143,15 +143,27 @@ int vtkMatrixMathFilter::RequestData
{
case DETERMINANT :
{
double tensor[9];
inTensors->GetTuple(i, tensor);
if (inTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
double const q = vtkMath::Determinant3x3(
reinterpret_cast<double(*)[3]>(inTensors->GetTuple(i)));
reinterpret_cast<double(*)[3]>(tensor));
quality->SetTuple(i, &q);
break;
}
case EIGENVALUE :
case EIGENVECTOR :
{
double* d = inTensors->GetTuple(i);
double d[9];
inTensors->GetTuple(i, d);
if (inTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(d);
}
double w[3]={0}, v[9]={0}, t[]={d[1]-d[3], d[2]-d[6], d[5]-d[7]};
// Use Jacobi iterative method only if the matrix is real symmetric.
......@@ -178,8 +190,15 @@ int vtkMatrixMathFilter::RequestData
}
case INVERSE :
{
double tensor[9];
inTensors->GetTuple(i, tensor);
if (inTensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(tensor);
}
double AI [3][3] = {{0}}, (*A) [3]
= reinterpret_cast<double(*)[3]>(inTensors->GetTuple(i));
= reinterpret_cast<double(*)[3]>(tensor);
// vtkMath::Invert3x3 should quite fit here, unfortunately, it does not
// check for matrix singularity which in the worest case leads to divide
......
......@@ -647,6 +647,16 @@ int vtkDataReader::ReadCellData(vtkDataSet *ds, int numCells)
}
}
//
// read 3x2 symmetric tensor data
//
else if ( ! strncmp(line, "tensors6", 8) )
{
if ( ! this->ReadTensorData(a, numCells, 6) )
{
return 0;
}
}
//
// read 3x3 tensor data
//
else if ( ! strncmp(line, "tensors", 7) )
......@@ -795,6 +805,16 @@ int vtkDataReader::ReadPointData(vtkDataSet *ds, int numPts)
}
}
//
// read 3x2 symmetric tensor data
//
else if ( ! strncmp(line, "tensors6", 8) )
{
if ( ! this->ReadTensorData(a, numPts, 6) )
{
return 0;
}
}
//
// read 3x3 tensor data
//
else if ( ! strncmp(line, "tensors", 7) )
......@@ -953,6 +973,16 @@ int vtkDataReader::ReadVertexData(vtkGraph *g, int numVertices)
}
}
//
// read 3x2 symmetric tensor data
//
else if ( ! strncmp(line, "tensors6", 8) )
{
if ( ! this->ReadTensorData(a, numVertices, 6) )
{
return 0;
}
}
//
// read 3x3 tensor data
//
else if ( ! strncmp(line, "tensors", 7) )
......@@ -1100,6 +1130,16 @@ int vtkDataReader::ReadEdgeData(vtkGraph *g, int numEdges)
}
}
//
// read 3x2 symmetric tensor data
//
else if ( ! strncmp(line, "tensors6", 8) )
{
if ( ! this->ReadTensorData(a, numEdges, 6) )
{
return 0;
}
}
//
// read 3x3 tensor data
//
else if ( ! strncmp(line, "tensors", 7) )
......@@ -1245,6 +1285,16 @@ int vtkDataReader::ReadRowData(vtkTable *t, int numEdges)
}
}
//
// read 3x2 symmetric tensor data
//
else if ( ! strncmp(line, "tensors6", 8) )
{
if ( ! this->ReadTensorData(a, numEdges, 6) )
{
return 0;
}
}
//
// read 3x3 tensor data
//
else if ( ! strncmp(line, "tensors", 7) )
......@@ -2300,7 +2350,7 @@ int vtkDataReader::ReadNormalData(vtkDataSetAttributes *a, int numPts)
}
// Read tensor point attributes. Return 0 if error.
int vtkDataReader::ReadTensorData(vtkDataSetAttributes *a, int numPts)
int vtkDataReader::ReadTensorData(vtkDataSetAttributes *a, int numPts, int numComp)
{
int skipTensor=0;
char line[256], name[256];
......@@ -2323,7 +2373,7 @@ int vtkDataReader::ReadTensorData(vtkDataSetAttributes *a, int numPts)
}
data = vtkArrayDownCast<vtkDataArray>(
this->ReadArray(line, numPts, 9));
this->ReadArray(line, numPts, numComp));
if ( data != NULL )
{
data->SetName(name);
......
......@@ -475,7 +475,7 @@ protected:
int ReadScalarData(vtkDataSetAttributes *a, int num);
int ReadVectorData(vtkDataSetAttributes *a, int num);
int ReadNormalData(vtkDataSetAttributes *a, int num);
int ReadTensorData(vtkDataSetAttributes *a, int num);
int ReadTensorData(vtkDataSetAttributes *a, int num, int numComp = 9);
int ReadCoScalarData(vtkDataSetAttributes *a, int num);
int ReadLutData(vtkDataSetAttributes *a);
int ReadTCoordsData(vtkDataSetAttributes *a, int num);
......
......@@ -1675,11 +1675,18 @@ int vtkDataWriter::WriteTensorData(ostream *fp, vtkDataArray *tensors, int num)
this->EncodeString(tensorsName, this->TensorsName, true);
}
*fp << "TENSORS ";
*fp << "TENSORS";
int numComp = 9;
if (tensors->GetNumberOfComponents() == 6)
{
*fp << "6";
numComp = 6;
}
*fp << " ";
sprintf(format, "%s %s\n", tensorsName, "%s");
delete[] tensorsName;
return this->WriteArray(fp, tensors->GetDataType(), tensors, format, num, 9);
return this->WriteArray(fp, tensors->GetDataType(), tensors, format, num, numComp);
}
int vtkDataWriter::WriteGlobalIdData(ostream *fp, vtkDataArray *globalIds, int num)
......
......@@ -122,9 +122,13 @@ void vtkEllipsoidTensorProbeRepresentation
{
tensors->GetTuple( this->ProbeCellId, t1 );
tensors->GetTuple( this->ProbeCellId+1, t2 );
if (tensors->GetNumberOfComponents() == 6)
{
vtkMath::TensorFromSymmetricTensor(t1);
vtkMath::TensorFromSymmetricTensor(t2);
}
}
// NN interpolation ?
// if ( r < 0.5 )
// {
......
be704b99bfd92b741b579176ba879fd0
06a49bf39aa3f3974ccb22473c150688
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment