diff --git a/Anima/diffusion/odf/odf_average/.vscode/settings.json b/Anima/diffusion/odf/odf_average/.vscode/settings.json new file mode 100644 index 000000000..6457ba81f --- /dev/null +++ b/Anima/diffusion/odf/odf_average/.vscode/settings.json @@ -0,0 +1,55 @@ +{ + "files.associations": { + "array": "cpp", + "atomic": "cpp", + "bit": "cpp", + "*.tcc": "cpp", + "cctype": "cpp", + "charconv": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "compare": "cpp", + "concepts": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdint": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "cwctype": "cpp", + "deque": "cpp", + "string": "cpp", + "unordered_map": "cpp", + "vector": "cpp", + "exception": "cpp", + "algorithm": "cpp", + "functional": "cpp", + "iterator": "cpp", + "memory": "cpp", + "memory_resource": "cpp", + "numeric": "cpp", + "optional": "cpp", + "random": "cpp", + "string_view": "cpp", + "system_error": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "format": "cpp", + "fstream": "cpp", + "initializer_list": "cpp", + "iosfwd": "cpp", + "istream": "cpp", + "limits": "cpp", + "new": "cpp", + "numbers": "cpp", + "ostream": "cpp", + "span": "cpp", + "stdexcept": "cpp", + "streambuf": "cpp", + "text_encoding": "cpp", + "typeinfo": "cpp", + "variant": "cpp" + } +} \ No newline at end of file diff --git a/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.cxx b/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.cxx index 43829ed67..68afbfd63 100644 --- a/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.cxx +++ b/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.cxx @@ -4,11 +4,37 @@ #include + namespace anima { + /** + * animaODFAverage computes the average of input ODFs voxel by voxel by using a Riemannian framework. + * Here we denote by ODF either a diffusion ODF, a fiber ODF or a TOD (see Dhollander et al., 2014 for the latter). + * This code has been deviced only for situations where at most 1 input ODF is not normalized (normalized = integrates to 1 over the sphere). + * It's the case when we compute an average between several TOD images (all ODF are normalized), or an average between a fiber ODF image and a TOD image (only the ODFs from the fiber ODF image are not normalized), which are the only two situations necessary for now. + * But using for example this code to compute the average between 2 fiber ODF images doesn't enter in the scope of the following algorithm (even if an output image can still be computed), as then more than 1 ODF is not normalized (in all voxels). + */ + + ODFAverageImageFilter::ODFAverageImageFilter() + { + m_EpsValueTestNormalized = std::pow(10, -6); + m_EpsValueTestNull = std::pow(10, -12); + m_VectorLength = 0; + m_ODFSHBasis = nullptr; + m_SamplePoints = g_SphericalDesignPoints; + //Put sample points (ie spherical designs) in spherical coordinates + std::for_each(m_SamplePoints.begin(), m_SamplePoints.end(), [](VectorType &v){anima::TransformCartesianToSphericalCoordinates(v,v);}); + } + + + ODFAverageImageFilter::~ODFAverageImageFilter() = default; + void ODFAverageImageFilter::AddWeightImage(const unsigned int i, const WeightImagePointer &weightImage) { + /** + * Adds a newly read weight image to the vector m_WeightImages + */ if (i == m_WeightImages.size()) { m_WeightImages.push_back(weightImage); @@ -16,157 +42,344 @@ namespace anima } if (i > m_WeightImages.size()) + { itkExceptionMacro("Weight images must be added contiguously."); + } m_WeightImages[i] = weightImage; } + void ODFAverageImageFilter::BeforeThreadedGenerateData() { + /** + * Preliminary computations + */ Superclass::BeforeThreadedGenerateData(); - unsigned int odfSHOrder = std::round(-1.5 + 0.5 * std::sqrt(8.0 * static_cast(this->GetInput(0)->GetVectorLength()) + 1.0)); - m_VectorLength = (odfSHOrder + 1) * (odfSHOrder + 2) / 2; + m_VectorLength = this->GetInput(0)->GetVectorLength(); + //Get the maximum spherical harmonic (SH) order from the vector length + unsigned int odfSHOrder = std::round(-1.5 + 0.5 * std::sqrt(8.0 * static_cast(m_VectorLength) + 1.0)); - m_SpherHarm.set_size(m_NbSamplesPhi * m_NbSamplesTheta, m_VectorLength); + //Build m_SpherHarm (it will contain the values of each spherical harmonic (columns) for each sample point (lines)) + m_SpherHarm.set_size(NB_SAMPLE_POINTS, m_VectorLength); m_ODFSHBasis = new anima::ODFSphericalHarmonicBasis(odfSHOrder); - - // Discretize SH double sqrt2 = std::sqrt(2.0); - double deltaPhi = 2.0 * M_PI / (static_cast(m_NbSamplesPhi) - 1.0); - double deltaTheta = M_PI / (static_cast(m_NbSamplesTheta) - 1.0); - - unsigned int k = 0; - for (unsigned int i = 0; i < m_NbSamplesTheta; ++i) + VectorType currentPoint; + for (unsigned int i = 0; i < NB_SAMPLE_POINTS; i++) { - double theta = static_cast(i) * deltaTheta; - - for (unsigned int j = 0; j < m_NbSamplesPhi; ++j) - { - double phi = static_cast(j) * deltaPhi; - unsigned int c = 0; - for (double l = 0; l <= odfSHOrder; l += 2) - for (double m = -l; m <= l; m++) - m_SpherHarm.put(k, c++, m_ODFSHBasis->getNthSHValueAtPosition(l, m, theta, phi)); - k++; - } + currentPoint = m_SamplePoints[i]; + unsigned int c = 0; + for (double l = 0; l <= odfSHOrder; l += 2) + for (double m = -l; m <= l; m++) + m_SpherHarm.put(i, c++, m_ODFSHBasis->getNthSHValueAtPosition(l, m, currentPoint[0], currentPoint[1])); } + //Pseudo-inverse computation (will be used to move from sampled representation to representation by coefficients) m_SolveSHMatrix = vnl_matrix_inverse(m_SpherHarm.transpose() * m_SpherHarm).as_matrix() * m_SpherHarm.transpose(); } + void ODFAverageImageFilter::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) { + /** + *Main function. For each voxel of the output image, it computes the average of the ODFs from the corresponding input voxels in input images. + */ + + //Initializations using InputImageIteratorType = itk::ImageRegionConstIterator; using InputWeightIteratorType = itk::ImageRegionConstIterator; using OutputImageIteratorType = itk::ImageRegionIterator; - unsigned int numInputs = this->GetNumberOfIndexedInputs(); + unsigned int numInputs = this->GetNumberOfIndexedInputs(); //number of input ODFs - std::vector inItrs(numInputs); - std::vector weightItrs(numInputs); + std::vector inItrs(numInputs); //vector of input images (which are here iterators on the input images) + std::vector weightItrs(numInputs); //vector of iterators on weight images for (unsigned int i = 0; i < numInputs; ++i) { inItrs[i] = InputImageIteratorType(this->GetInput(i), outputRegionForThread); weightItrs[i] = InputWeightIteratorType(m_WeightImages[i], outputRegionForThread); } - - OutputImageIteratorType outItr(this->GetOutput(), outputRegionForThread); - - VectorType weightValues(numInputs); - VectorType workValue(m_NbSamplesPhi * m_NbSamplesTheta); - HistoArrayType workValues(numInputs); - - InputPixelType inputValue(m_VectorLength); + OutputImageIteratorType outItr(this->GetOutput(), outputRegionForThread); //iterator on output image + + VectorType coefODFValue(m_VectorLength); + VectorType sampledODFValue(NB_SAMPLE_POINTS); + Vector2DType currentODFs(numInputs); + Vector2DType notNullCurrentODFs; + VectorType currentWeights(numInputs); + VectorType currentWeightsNotNullODFs; + VectorType inputValue(m_VectorLength); OutputPixelType outputValue(m_VectorLength); - double epsValue = std::sqrt(std::numeric_limits::epsilon()); - + //Average computation while (!outItr.IsAtEnd()) { double weightSum = 0.0; + double normValue = 1.0; + int numNotNullODFs; + int indexODFtoNormalize; + + //Compute sum of weights of all input images for the current voxel for (unsigned int i = 0; i < numInputs; ++i) { double weightValue = weightItrs[i].Get(); weightSum += weightValue; - weightValues[i] = weightValue; + currentWeights[i] = weightValue; } - if (weightSum < epsValue) + if (weightSum < m_EpsValueTestNull) + //If in a voxel, the sum of weights is 0, then we output a null distribution { outputValue.Fill(0.0); outItr.Set(outputValue); + //Move to next voxel for (unsigned int i = 0; i < numInputs; ++i) { ++inItrs[i]; ++weightItrs[i]; } - ++outItr; + this->IncrementNumberOfProcessedPoints(); //to see the progression + } + + else + { + for (unsigned int i = 0; i < numInputs; ++i) + { + //Read input ODFs in the current voxel position and change their type + this->ConvertODF(inItrs[i].Get(), inputValue); + currentODFs[i] = inputValue; + } + + //Remove null ODFs and their associated weights to avoid useless computations + numNotNullODFs = this->RemoveNullODFs(currentODFs, currentWeights, numInputs, notNullCurrentODFs, currentWeightsNotNullODFs); + switch (numNotNullODFs) + { + case 0: + //Here, all input ODFs are null, so we simply output a null distribution + outputValue.Fill(0.0); + break; + + case 1: + //Here, only 1 input ODF is not null, so we simply output this ODF + for (unsigned int i = 0; i < m_VectorLength; i++) + { + outputValue[i] = notNullCurrentODFs[0][i]; + } + break; + + default: + //Here, more than 1 ODF is not null, so we need to compute an average with the Riemannian framework + indexODFtoNormalize = this->GetIndexODFtoNormalize(notNullCurrentODFs, numNotNullODFs); + for (unsigned int i = 0; i < numNotNullODFs; i++) + { + this->DiscretizeODF(notNullCurrentODFs[i], sampledODFValue); //it enables to compute easily the square root + coefODFValue = this->GetSquareRootODFCoef(sampledODFValue); //the averaging needs to be performed on square roots of ODFs (see Goh et al., 2011) + if (i == indexODFtoNormalize) + { + //Here, the current ODF needs to be normalized (it's the only one) + normValue = this->NormalizeODF(coefODFValue, notNullCurrentODFs[i]); + } + else + { + notNullCurrentODFs[i] = coefODFValue; + } + } + //notNullCurrentODFs now contains the ODFs in the correct form for averaging + coefODFValue = this->ComputeAverage(notNullCurrentODFs, currentWeightsNotNullODFs, numNotNullODFs); + if (indexODFtoNormalize != -1) + { + //Here, one input ODF has been normalized, so we reverse the normalization to retrieve the initial ODF amplitude + this->InverseNormalization(normValue, coefODFValue); + } + this->DiscretizeODF(coefODFValue, sampledODFValue); //it enables to compute easily the square + outputValue = this->GetSquareODFCoef(sampledODFValue); //reverse the square root operation + } + + //Reset vectors + notNullCurrentODFs.clear(); + currentWeightsNotNullODFs.clear(); + + //Set output value and move to next voxel + outItr.Set(outputValue); + for (unsigned int i = 0; i < numInputs; ++i) + { + ++inItrs[i]; + ++weightItrs[i]; + } + ++outItr; this->IncrementNumberOfProcessedPoints(); - continue; } + } + } + - for (unsigned int i = 0; i < numInputs; ++i) + void ODFAverageImageFilter::AfterThreadedGenerateData() + { + /** + * Deletes pointers which are not SmartPointers + */ + delete m_ODFSHBasis; + } + + + void ODFAverageImageFilter::ConvertODF(const InputPixelType &inputODF, VectorType &convertedODF) + { + /** + * Converts an ODF from InputPixelType to VectorType + */ + for (unsigned int i = 0; i < m_VectorLength; i++) + { + convertedODF[i] = inputODF[i]; + } + } + + + int ODFAverageImageFilter::RemoveNullODFs(const Vector2DType &allODFs, const VectorType &allWeights, int nbODFs, Vector2DType &selectedODFs, VectorType &selectedWeights) + { + /** + * allODFs contains all input ODFs in a given voxel position and allWeights contains the associated weights. nbODFs is the number of input ODFs. + * selectedODFs will contain only the input ODFs which are not null, and we will keep only their associated weights in selectedWeights. + * The function also returns the number of not null ODFs. + */ + int numNotNullODFs = 0; + for (unsigned int i = 0; i < nbODFs; i++) + { + VectorType odf = allODFs[i]; + + //Check is the ODF is null (namely all its coefficients are 0) + bool currentODFisNull = true; + int j = 0; + while (currentODFisNull && jDiscretizeODF(inputValue, workValue); - workValues[i] = this->GetSquareRootODFCoef(workValue); + if (!std::abs(odf[j])GetAverageHisto(workValues, weightValues, outputValue); - this->DiscretizeODF(outputValue, workValue); - outputValue = this->GetSquareODFCoef(workValue); + //Keep the ODF if it's not null + if (!currentODFisNull) + { + selectedODFs.push_back(odf); + selectedWeights.push_back(allWeights[i]); + numNotNullODFs++; + } + } + return numNotNullODFs; + } - outItr.Set(outputValue); - for (unsigned int i = 0; i < numInputs; ++i) + int ODFAverageImageFilter::GetIndexODFtoNormalize(Vector2DType &odfs, int &numNotNullODFs) + { + /** + * Gets the position of the potential ODF that needs to be normalized inside the vector of ODFs odfs. + * Returns -1 if all input ODFs are already normalized. + * + */ + int index = -1; + int i=0; + while (index == -1 && i < numNotNullODFs) + { + if (std::abs(odfs[i][0]-1/std::sqrt(4*M_PI)) >= m_EpsValueTestNormalized) + //An ODF is normalized if and only if its first coefficient is equal to 1/sqrt(4*PI) { - ++inItrs[i]; - ++weightItrs[i]; + index = i; } + i++; + } + return index; + + } - ++outItr; - this->IncrementNumberOfProcessedPoints(); + double ODFAverageImageFilter::NormalizeODF(const VectorType &inputODF, VectorType &normalizedODF) + { + /** + * Normalizes an ODF so that its square integrates to 1 over the sphere + * (This function is applied to square root of input ODFs, so the aim is now to have a square integral equal to 1) + */ + double sumSquareCoefs = 0; + double normValue; + + //Compute the square integral of the ODF (which is simply the square sum of its coefficients, since SH form an orthornormal basis) + for (unsigned int i = 0; i < m_VectorLength; i++) + { + sumSquareCoefs+=std::pow(inputODF[i],2); } + //If the square sum is null, we output null vector + if (std::abs(sumSquareCoefs)getValueAtPosition(coefODF, currentPoint[0], currentPoint[1]); //currentPoint[0] : theta coordinate. currentPoint[1] : phi coordinate + } } - void ODFAverageImageFilter::DiscretizeODF(const InputPixelType &modelValue, VectorType &odf) - { - unsigned int k = 0; - double deltaPhi = 2.0 * M_PI / (static_cast(m_NbSamplesPhi) - 1.0); - double deltaTheta = M_PI / (static_cast(m_NbSamplesTheta) - 1.0); - for (unsigned int i = 0; i < m_NbSamplesTheta; ++i) + void ODFAverageImageFilter::InverseNormalization(const double &normValue, VectorType &odf) + { + /**If 1 (and so exactly 1) input ODF has been normalized by multiplication by a scalar normValue, we here divide the output average ODF by normValue to retrieve the initial ODF amplitude. + * In our case, if an input ODF has been normalized, it's the fiber ODF before an averaging with a TOD. And indeed, in that case, the amplitudes of the output average ODF and its peaks must be the ones of the input fiber ODF. + */ + if (!std::abs(normValue) < m_EpsValueTestNull) { - double theta = static_cast(i) * deltaTheta; - for (unsigned int j = 0; j < m_NbSamplesPhi; ++j) + for (unsigned int j = 0; j < m_VectorLength; j++) { - double phi = static_cast(j) * deltaPhi; - odf[k] = m_ODFSHBasis->getValueAtPosition(modelValue, theta, phi); - k++; - } + odf[j] = odf[j]/normValue; + } } } + ODFAverageImageFilter::VectorType ODFAverageImageFilter::GetSquareRootODFCoef(const VectorType &odf) { - MatrixType squareRootOdf(m_NbSamplesTheta * m_NbSamplesPhi, 1); + /** + * Returns the square root of an input ODF + * Output ODF is represented by coefficients + * Input ODF is in sample form + */ + + MatrixType squareRootOdf(NB_SAMPLE_POINTS, 1); MatrixType squareRootCoef(m_VectorLength, 1); - for (unsigned int i = 0; i < m_NbSamplesTheta * m_NbSamplesPhi; ++i) + //Compute the square root of each sample point + for (unsigned int i = 0; i < NB_SAMPLE_POINTS; ++i) squareRootOdf.put(i, 0, std::sqrt(std::max(0.0, odf[i]))); + //Convert from sample form to representation by coefficients squareRootCoef = m_SolveSHMatrix * squareRootOdf; + //Store the result in a vector instead of a matrix VectorType squareRootModelValue(m_VectorLength); for (unsigned int i = 0; i < m_VectorLength; ++i) squareRootModelValue[i] = squareRootCoef.get(i, 0); @@ -174,16 +387,26 @@ namespace anima return squareRootModelValue; } + ODFAverageImageFilter::OutputPixelType ODFAverageImageFilter::GetSquareODFCoef(const VectorType &odf) { - MatrixType squareOdf(m_NbSamplesTheta * m_NbSamplesPhi, 1); + /** + * Returns the square of an input ODF + * Output ODF is represented by coefficients + * Input ODF is in sample form + */ + + MatrixType squareOdf(NB_SAMPLE_POINTS, 1); MatrixType squareCoef(m_VectorLength, 1); - for (unsigned int i = 0; i < m_NbSamplesTheta * m_NbSamplesPhi; ++i) - squareOdf.put(i, 0, std::pow(odf[i], 2.0)); + //Compute the square of each sample point + for (unsigned int i = 0; i < NB_SAMPLE_POINTS; ++i) + squareOdf.put(i, 0, std::pow(odf[i], 2)); + //Convert from sample form to representation by coefficients squareCoef = m_SolveSHMatrix * squareOdf; + //Store the result in a vector instead of a matrix OutputPixelType squareModelValue(m_VectorLength); for (unsigned int i = 0; i < m_VectorLength; ++i) squareModelValue[i] = squareCoef.get(i, 0); @@ -191,16 +414,21 @@ namespace anima return squareModelValue; } - void ODFAverageImageFilter::GetAverageHisto(const HistoArrayType &coefs, const VectorType &weightValues, OutputPixelType &resCoef) + + ODFAverageImageFilter::VectorType ODFAverageImageFilter::ComputeAverage(const Vector2DType &coefs, const VectorType &weightValues, int nbODFs) { - unsigned int numImages = this->GetNumberOfIndexedInputs(); + /** + * Computes the weighted average of input ODFs using the method of Goh et al., 2O11 (Algorithm 1) + * coefs: Vector of size nbODFs, each element being a vector of size m_VectorLength representing an ODF + * weightValues: Vector of size nbODFs with the weights of the input ODFs + * Returns the average ODF represented by coefficients + */ VectorType mean, nextMean = coefs[0]; VectorType tangent(m_VectorLength), workValue(m_VectorLength); const unsigned int maxIter = 100; const double epsValue = 0.000035; - unsigned int nIter = 0; double normTan = 1.0; @@ -210,7 +438,7 @@ namespace anima std::fill(tangent.begin(), tangent.end(), 0.0); double weightSum = 0.0; - for (unsigned int i = 0; i < numImages; ++i) + for (unsigned int i = 0; i < nbODFs; ++i) { weightSum += weightValues[i]; std::fill(workValue.begin(), workValue.end(), 0.0); @@ -220,8 +448,7 @@ namespace anima } for (unsigned int j = 0; j < m_VectorLength; ++j) - tangent[j] /= weightSum; - + tangent[j] /= weightSum; //enables to make averaging work even if the sum of the weights is not equal to 1 normTan = anima::ComputeNorm(tangent); std::fill(nextMean.begin(), nextMean.end(), 0.0); @@ -230,8 +457,6 @@ namespace anima nIter++; } - for (unsigned int i = 0; i < m_VectorLength; ++i) - resCoef[i] = nextMean[i]; + return nextMean; } - -} // end namespace anima +} //end namespace anima diff --git a/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.h b/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.h index 3e272b51e..3f33d3a1c 100644 --- a/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.h +++ b/Anima/diffusion/odf/odf_average/animaODFAverageImageFilter.h @@ -1,19 +1,24 @@ #pragma once +#include "animaSphericalDesignPoints.h" #include #include #include + namespace anima { - class ODFAverageImageFilter : public anima::NumberedThreadImageToImageFilter, itk::VectorImage> + class ODFAverageImageFilter : public anima::NumberedThreadImageToImageFilter, itk::VectorImage> { public: + ODFAverageImageFilter(); + ~ODFAverageImageFilter() override; + using Self = ODFAverageImageFilter; - using InputImageType = itk::VectorImage; - using OutputImageType = itk::VectorImage; - using WeightImageType = itk::Image; + using InputImageType = itk::VectorImage; + using OutputImageType = itk::VectorImage; + using WeightImageType = itk::Image; using Superclass = anima::NumberedThreadImageToImageFilter; using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; @@ -24,12 +29,10 @@ namespace anima /** Run-time type information (and related methods) */ itkTypeMacro(ODFAverageImageFilter, anima::NumberedThreadImageToImageFilter); - /** Image typedef support */ using InputImagePointer = InputImageType::Pointer; using OutputImagePointer = OutputImageType::Pointer; using WeightImagePointer = WeightImageType::Pointer; - /** Superclass typedefs. */ using InputImageRegionType = Superclass::InputImageRegionType; using OutputImageRegionType = Superclass::OutputImageRegionType; using InputPixelType = InputImageType::PixelType; @@ -37,44 +40,38 @@ namespace anima /** Typedefs for computations. */ using VectorType = std::vector; - using HistoArrayType = std::vector; + using Vector2DType = std::vector; using MatrixType = vnl_matrix; void AddWeightImage(const unsigned int i, const WeightImagePointer &weightImage); protected: - ODFAverageImageFilter() : Superclass() - { - m_NbSamplesTheta = 10; - m_NbSamplesPhi = 2 * m_NbSamplesTheta; - m_WeightImages.clear(); - m_SolveSHMatrix.clear(); - m_SpherHarm.clear(); - } - - virtual ~ODFAverageImageFilter() {} - - void BeforeThreadedGenerateData() ITK_OVERRIDE; - void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE; - void AfterThreadedGenerateData() ITK_OVERRIDE; - - void DiscretizeODF(const InputPixelType &Coef, VectorType &resHisto); - void GetAverageHisto(const HistoArrayType &coefs, const VectorType &weightValues, OutputPixelType &resCoef); - VectorType GetSquareRootODFCoef(const VectorType &histo); - OutputPixelType GetSquareODFCoef(const VectorType &histo); + void BeforeThreadedGenerateData() override; + void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) override; + void AfterThreadedGenerateData() override; + + void ConvertODF(const InputPixelType &inputODF, VectorType &convertedODF); + int RemoveNullODFs(const Vector2DType &allODFs, const VectorType &allweights, int nbODFs, Vector2DType &selectedODFs, VectorType &selectedWeights); + int GetIndexODFtoNormalize(Vector2DType &odfs, int &numNotNullODFs); + void DiscretizeODF(const VectorType &coefODF, VectorType &sampledODF); + double NormalizeODF(const VectorType &inputODF, VectorType &normalizedODF); + void InverseNormalization(const double &normValue, VectorType &odf); + VectorType GetSquareRootODFCoef(const VectorType &odf); + OutputPixelType GetSquareODFCoef(const VectorType &odf); + VectorType ComputeAverage(const Vector2DType &coefs, const VectorType &weightValues, int nbODFs); private: ITK_DISALLOW_COPY_AND_ASSIGN(ODFAverageImageFilter); - unsigned int m_NbSamplesPhi; - unsigned int m_NbSamplesTheta; unsigned int m_VectorLength; + double m_EpsValueTestNormalized; //Will be used to test if an odf is normalized + double m_EpsValueTestNull; //Will be used to test if values are null (needs to be much smaller than m_EpsValueTestNormalized) std::vector m_WeightImages; - MatrixType m_SpherHarm; MatrixType m_SolveSHMatrix; - anima::ODFSphericalHarmonicBasis *m_ODFSHBasis; + + Array2DType m_SamplePoints; }; } // end of namespace anima diff --git a/Anima/diffusion/odf/odf_average/animaSphericalDesignPoints.h b/Anima/diffusion/odf/odf_average/animaSphericalDesignPoints.h new file mode 100644 index 000000000..e96529eac --- /dev/null +++ b/Anima/diffusion/odf/odf_average/animaSphericalDesignPoints.h @@ -0,0 +1,254 @@ +#pragma once + +#include +#include + +#define NB_SAMPLE_POINTS 240 + +using Array2DType = std::array, NB_SAMPLE_POINTS>; + +/*The following 3D points will be used to get a sampled representation of the ODF (namely a vector with the ODF's values on each of these points). +These points are called spherical design. They enable to compute sums on the unit sphere instead of integrals without any loss of information (see Hardin and Sloane, 1995).*/ +static const Array2DType g_SphericalDesignPoints ( +{{ + { 0.8926535357627230E+00, 0.4125340536573610E+00, -0.1816186104542530E+00}, + { 0.4125340534250320E+00, -0.1816186106417820E+00, 0.8926535358319380E+00}, + {-0.8926535358064071E+00, -0.4125340536278530E+00, -0.1816186103065750E+00}, + {-0.1816186106138490E+00, 0.8926535357404750E+00, 0.4125340536352400E+00}, + {-0.4125340534774350E+00, -0.1816186104226540E+00, -0.8926535358523040E+00}, + {-0.1816186104513840E+00, -0.8926535357628120E+00, -0.4125340536584320E+00}, + {-0.4125340533170900E+00, 0.1816186106118270E+00, 0.8926535358879180E+00}, + { 0.1816186104001360E+00, -0.8926535358123000E+00, 0.4125340535739110E+00}, + { 0.4125340533279960E+00, 0.1816186104204000E+00, -0.8926535359218249E+00}, + { 0.1816186105807890E+00, 0.8926535358109040E+00, -0.4125340534973990E+00}, + { 0.8926535358676440E+00, -0.4125340534725580E+00, 0.1816186103583390E+00}, + {-0.8926535358550640E+00, 0.4125340535351600E+00, 0.1816186102779710E+00}, + {-0.2920937425934330E+00, -0.2957670279931700E+00, 0.9095070701703470E+00}, + {-0.2957670280268870E+00, 0.9095070700892600E+00, -0.2920937428117760E+00}, + { 0.2920937424478640E+00, 0.2957670280397130E+00, 0.9095070702019620E+00}, + { 0.9095070701476120E+00, -0.2920937429267210E+00, -0.2957670277339340E+00}, + { 0.2957670281453960E+00, 0.9095070700844410E+00, 0.2920937427067830E+00}, + { 0.9095070701888540E+00, 0.2920937426892070E+00, 0.2957670278416750E+00}, + { 0.2957670279073110E+00, -0.9095070701484190E+00, -0.2920937427486510E+00}, + {-0.9095070701012210E+00, 0.2920937431592720E+00, -0.2957670276469270E+00}, + {-0.2957670278353330E+00, -0.9095070700472930E+00, 0.2920937431364140E+00}, + {-0.9095070702185910E+00, -0.2920937427217760E+00, 0.2957670277180690E+00}, + {-0.2920937425408960E+00, 0.2957670277931470E+00, -0.9095070702522660E+00}, + { 0.2920937428619380E+00, -0.2957670277476140E+00, -0.9095070701639690E+00}, + {-0.5752257180381920E+00, 0.2412057282507800E-01, 0.8176390225974030E+00}, + { 0.2412057278614400E-01, 0.8176390225112380E+00, -0.5752257181623009E+00}, + { 0.5752257181164780E+00, -0.2412057297921300E-01, 0.8176390225377810E+00}, + { 0.8176390225560030E+00, -0.5752257181034800E+00, 0.2412057267146900E-01}, + {-0.2412057304150300E-01, 0.8176390224407570E+00, 0.5752257182517770E+00}, + { 0.8176390224583791E+00, 0.5752257182291181E+00, -0.2412057298452600E-01}, + {-0.2412057281823900E-01, -0.8176390225812600E+00, -0.5752257180614240E+00}, + {-0.8176390225435780E+00, 0.5752257181238820E+00, 0.2412057260611100E-01}, + { 0.2412057271295000E-01, -0.8176390225272960E+00, 0.5752257181425460E+00}, + {-0.8176390226004950E+00, -0.5752257180351740E+00, -0.2412057279222800E-01}, + {-0.5752257179254689E+00, -0.2412057271105200E-01, -0.8176390226800700E+00}, + { 0.5752257179082300E+00, 0.2412057259415500E-01, -0.8176390226956460E+00}, + {-0.1288331617248000E+00, 0.5224764072024000E-01, 0.9902889479738530E+00}, + { 0.5224764069440900E-01, 0.9902889479588950E+00, -0.1288331618502510E+00}, + { 0.1288331618403250E+00, -0.5224764032003800E-01, 0.9902889479799381E+00}, + { 0.9902889479497170E+00, -0.1288331619247960E+00, 0.5224764068455800E-01}, + {-0.5224764038851000E-01, 0.9902889479675810E+00, 0.1288331619075380E+00}, + { 0.9902889479777300E+00, 0.1288331618780010E+00, -0.5224764026899200E-01}, + {-0.5224764039040900E-01, -0.9902889479621900E+00, -0.1288331619482090E+00}, + {-0.9902889479606261E+00, 0.1288331618966490E+00, 0.5224764054718700E-01}, + { 0.5224764052780800E-01, -0.9902889479532510E+00, 0.1288331619612000E+00}, + {-0.9902889479708680E+00, -0.1288331619362050E+00, -0.5224764025552600E-01}, + {-0.1288331617904780E+00, -0.5224764033764300E-01, -0.9902889479854941E+00}, + { 0.1288331618574160E+00, 0.5224764055154500E-01, -0.9902889479655000E+00}, + { 0.7180063860347500E+00, 0.6574468762559930E+00, -0.2285397875962860E+00}, + { 0.6574468762867370E+00, -0.2285397878319220E+00, 0.7180063859315960E+00}, + {-0.7180063861094420E+00, -0.6574468761714340E+00, -0.2285397876048770E+00}, + {-0.2285397877372190E+00, 0.7180063859474221E+00, 0.6574468763023740E+00}, + {-0.6574468762410210E+00, -0.2285397877138000E+00, -0.7180063860110540E+00}, + {-0.2285397876789970E+00, -0.7180063860313590E+00, -0.6574468762309450E+00}, + {-0.6574468763611850E+00, 0.2285397878605490E+00, 0.7180063858543150E+00}, + { 0.2285397877030650E+00, -0.7180063858573850E+00, 0.6574468764125770E+00}, + { 0.6574468763044540E+00, 0.2285397878740170E+00, -0.7180063859019750E+00}, + { 0.2285397877849670E+00, 0.7180063858138530E+00, -0.6574468764316480E+00}, + { 0.7180063858807600E+00, -0.6574468763634850E+00, 0.2285397877708510E+00}, + {-0.7180063858910180E+00, 0.6574468763715580E+00, 0.2285397877154010E+00}, + { 0.8631764731178030E+00, 0.4681818166531380E+00, 0.1890295289400010E+00}, + { 0.4681818164384860E+00, 0.1890295291974920E+00, 0.8631764731778400E+00}, + {-0.8631764731944460E+00, -0.4681818165764200E+00, 0.1890295287800330E+00}, + { 0.1890295291255270E+00, 0.8631764730643890E+00, 0.4681818166767080E+00}, + {-0.4681818163926710E+00, 0.1890295288974430E+00, -0.8631764732683980E+00}, + { 0.1890295287921740E+00, -0.8631764731436880E+00, -0.4681818166651000E+00}, + {-0.4681818164112130E+00, -0.1890295291281380E+00, 0.8631764732078210E+00}, + {-0.1890295288978520E+00, -0.8631764730897200E+00, 0.4681818167219310E+00}, + { 0.4681818165088670E+00, -0.1890295289305550E+00, -0.8631764731981230E+00}, + {-0.1890295290018230E+00, 0.8631764731066590E+00, -0.4681818166487220E+00}, + { 0.8631764731352290E+00, -0.4681818166486420E+00, -0.1890295288715610E+00}, + {-0.8631764731233340E+00, 0.4681818166987620E+00, -0.1890295288017440E+00}, + { 0.7726328568471330E+00, -0.5170594506955900E+00, 0.3683585114621520E+00}, + {-0.5170594505671320E+00, 0.3683585115855150E+00, 0.7726328568742860E+00}, + {-0.7726328568060810E+00, 0.5170594506473910E+00, 0.3683585116159150E+00}, + { 0.3683585116480010E+00, 0.7726328568060540E+00, -0.5170594506245730E+00}, + { 0.5170594504940070E+00, 0.3683585118165880E+00, -0.7726328568130560E+00}, + { 0.3683585117204960E+00, -0.7726328568024760E+00, 0.5170594505782730E+00}, + { 0.5170594505834450E+00, -0.3683585114871170E+00, 0.7726328569102801E+00}, + {-0.3683585115673300E+00, -0.7726328568594670E+00, -0.5170594506022290E+00}, + {-0.5170594505023690E+00, -0.3683585116659560E+00, -0.7726328568792750E+00}, + {-0.3683585114698030E+00, 0.7726328568556510E+00, 0.5170594506774120E+00}, + { 0.7726328569347490E+00, 0.5170594506919191E+00, -0.3683585112835310E+00}, + {-0.7726328569274850E+00, -0.5170594506337780E+00, -0.3683585113803780E+00}, + {-0.8478192319146480E+00, -0.6632577590016700E-01, -0.5261211281130020E+00}, + {-0.6632577591363099E-01, -0.5261211282576860E+00, -0.8478192318238090E+00}, + { 0.8478192318830180E+00, 0.6632577581985200E-01, -0.5261211281740970E+00}, + {-0.5261211283487620E+00, -0.8478192317669569E+00, -0.6632577591791000E-01}, + { 0.6632577584612000E-01, -0.5261211284070980E+00, 0.8478192317363720E+00}, + {-0.5261211284592400E+00, 0.8478192317090800E+00, 0.6632577578136600E-01}, + { 0.6632577594578500E-01, 0.5261211283443800E+00, -0.8478192317674960E+00}, + { 0.5261211284495320E+00, 0.8478192317006920E+00, -0.6632577596561300E-01}, + {-0.6632577587721100E-01, 0.5261211283063880E+00, 0.8478192317964360E+00}, + { 0.5261211285046690E+00, -0.8478192316652130E+00, 0.6632577598176000E-01}, + {-0.8478192318217250E+00, 0.6632577594100500E-01, 0.5261211282575941E+00}, + { 0.8478192318502640E+00, -0.6632577599665500E-01, 0.5261211282045900E+00}, + { 0.9805743229230000E-02, 0.9429838158425931E+00, 0.3326941094438920E+00}, + { 0.9429838158089230E+00, 0.3326941095397480E+00, 0.9805743214949999E-02}, + {-0.9805743379690000E-02, -0.9429838157872910E+00, 0.3326941095962070E+00}, + { 0.3326941092265540E+00, 0.9805743204272001E-02, 0.9429838159195320E+00}, + {-0.9429838157740400E+00, 0.3326941096356470E+00, -0.9805743315804000E-02}, + { 0.3326941093979960E+00, -0.9805743298910000E-02, -0.9429838158580620E+00}, + {-0.9429838157761140E+00, -0.3326941096300980E+00, 0.9805743304667000E-02}, + {-0.3326941093190270E+00, -0.9805743188507000E-02, 0.9429838158870700E+00}, + { 0.9429838157750819E+00, -0.3326941096351990E+00, -0.9805743230762999E-02}, + {-0.3326941094557650E+00, 0.9805743389761999E-02, -0.9429838158367350E+00}, + { 0.9805743301140001E-02, -0.9429838157525240E+00, -0.3326941096970650E+00}, + {-0.9805743287713000E-02, 0.9429838157913790E+00, -0.3326941095873310E+00}, + { 0.7855992483711520E+00, -0.4051569453122690E+00, -0.4676341204658960E+00}, + {-0.4051569449321250E+00, -0.4676341206498590E+00, 0.7855992484576980E+00}, + {-0.7855992482017899E+00, 0.4051569454340510E+00, -0.4676341206449040E+00}, + {-0.4676341206112420E+00, 0.7855992483346230E+00, -0.4051569452153390E+00}, + { 0.4051569451364230E+00, -0.4676341208682010E+00, -0.7855992482223660E+00}, + {-0.4676341208118040E+00, -0.7855992481456090E+00, 0.4051569453503470E+00}, + { 0.4051569448419850E+00, 0.4676341208613320E+00, 0.7855992483783050E+00}, + { 0.4676341207867260E+00, -0.7855992482498571E+00, -0.4051569451771560E+00}, + {-0.4051569449996430E+00, 0.4676341208710980E+00, -0.7855992482911820E+00}, + { 0.4676341208937130E+00, 0.7855992482342400E+00, 0.4051569450839530E+00}, + { 0.7855992483133410E+00, 0.4051569451171040E+00, 0.4676341207321060E+00}, + {-0.7855992482810999E+00, -0.4051569451973700E+00, 0.4676341207167270E+00}, + {-0.7373319991314921E+00, 0.6208515010137640E+00, -0.2662422519918900E+00}, + { 0.6208515009491860E+00, -0.2662422521548950E+00, -0.7373319991270100E+00}, + { 0.7373319990600610E+00, -0.6208515010887370E+00, -0.2662422520148830E+00}, + {-0.2662422519486310E+00, -0.7373319991032550E+00, 0.6208515010658500E+00}, + {-0.6208515010792210E+00, -0.2662422522338000E+00, 0.7373319989890250E+00}, + {-0.2662422520116240E+00, 0.7373319989962220E+00, -0.6208515011659510E+00}, + {-0.6208515010721239E+00, 0.2662422522225600E+00, -0.7373319989990600E+00}, + { 0.2662422521138640E+00, 0.7373319988329740E+00, 0.6208515013159830E+00}, + { 0.6208515011873870E+00, 0.2662422523283740E+00, 0.7373319988637970E+00}, + { 0.2662422519322500E+00, -0.7373319989389900E+00, -0.6208515012679590E+00}, + {-0.7373319989479430E+00, -0.6208515011832970E+00, 0.2662422521048790E+00}, + { 0.7373319988350070E+00, 0.6208515013057860E+00, 0.2662422521320100E+00}, + { 0.7268714691656590E+00, -0.2748828235042800E-01, -0.6862231864680610E+00}, + {-0.2748828218275500E-01, -0.6862231864483250E+00, 0.7268714691906329E+00}, + {-0.7268714691729310E+00, 0.2748828237188500E-01, -0.6862231864594990E+00}, + {-0.6862231864497120E+00, 0.7268714691854060E+00, -0.2748828228634100E-01}, + { 0.2748828235160700E-01, -0.6862231864911200E+00, -0.7268714691438450E+00}, + {-0.6862231865456220E+00, -0.7268714690897941E+00, 0.2748828242028100E-01}, + { 0.2748828226683600E-01, 0.6862231864703350E+00, 0.7268714691666740E+00}, + { 0.6862231866611830E+00, -0.7268714689834220E+00, -0.2748828234818500E-01}, + {-0.2748828225102900E-01, 0.6862231865230920E+00, -0.7268714691174650E+00}, + { 0.6862231866091120E+00, 0.7268714690334980E+00, 0.2748828232394800E-01}, + { 0.7268714690701070E+00, 0.2748828233555000E-01, 0.6862231865698690E+00}, + {-0.7268714690801830E+00, -0.2748828230971600E-01, 0.6862231865602320E+00}, + { 0.6653633857205150E+00, 0.5808602677392710E+00, 0.4689274083527160E+00}, + { 0.5808602675770870E+00, 0.4689274084886380E+00, 0.6653633857663080E+00}, + {-0.6653633856773800E+00, -0.5808602677195750E+00, 0.4689274084383180E+00}, + { 0.4689274083407830E+00, 0.6653633858218631E+00, 0.5808602676328130E+00}, + {-0.5808602675284530E+00, 0.4689274086788320E+00, -0.6653633856747230E+00}, + { 0.4689274083726140E+00, -0.6653633856988030E+00, -0.5808602677480780E+00}, + {-0.5808602676408769E+00, -0.4689274085527620E+00, 0.6653633856654270E+00}, + {-0.4689274084683360E+00, -0.6653633858479470E+00, 0.5808602674999610E+00}, + { 0.5808602673867520E+00, -0.4689274086545190E+00, -0.6653633858155630E+00}, + {-0.4689274083756990E+00, 0.6653633856513560E+00, -0.5808602677999380E+00}, + { 0.6653633856518190E+00, -0.5808602677912120E+00, -0.4689274083858500E+00}, + {-0.6653633857517340E+00, 0.5808602675480170E+00, -0.4689274085453260E+00}, + {-0.5801253673053040E+00, -0.7790995979244340E+00, 0.2376097109187070E+00}, + {-0.7790995980535180E+00, 0.2376097109099340E+00, -0.5801253671355390E+00}, + { 0.5801253671868080E+00, 0.7790995979777320E+00, 0.2376097110332580E+00}, + { 0.2376097106959320E+00, -0.5801253672761100E+00, -0.7790995980141140E+00}, + { 0.7790995980647319E+00, 0.2376097111473200E+00, 0.5801253670232500E+00}, + { 0.2376097108192850E+00, 0.5801253670474260E+00, 0.7790995981467740E+00}, + { 0.7790995981702240E+00, -0.2376097108496420E+00, -0.5801253670034990E+00}, + {-0.2376097108118020E+00, 0.5801253671572560E+00, -0.7790995980672760E+00}, + {-0.7790995980749610E+00, -0.2376097110451280E+00, 0.5801253670513690E+00}, + {-0.2376097106092530E+00, -0.5801253670223590E+00, 0.7790995982294950E+00}, + {-0.5801253670900940E+00, 0.7790995981519661E+00, -0.2376097106980860E+00}, + { 0.5801253672184110E+00, -0.7790995979667160E+00, -0.2376097109922150E+00}, + { 0.9586680253602000E+00, 0.1011136059005390E+00, -0.2659542363899560E+00}, + { 0.1011136058898930E+00, -0.2659542364771990E+00, 0.9586680253371200E+00}, + {-0.9586680253264100E+00, -0.1011136060954320E+00, -0.2659542364376600E+00}, + {-0.2659542366341790E+00, 0.9586680252945550E+00, 0.1011136058805580E+00}, + {-0.1011136060031710E+00, -0.2659542366563170E+00, -0.9586680252754820E+00}, + {-0.2659542367154550E+00, -0.9586680252461620E+00, -0.1011136061256020E+00}, + {-0.1011136058254380E+00, 0.2659542364146640E+00, 0.9586680253612671E+00}, + { 0.2659542362867390E+00, -0.9586680253935830E+00, 0.1011136058555220E+00}, + { 0.1011136058024440E+00, 0.2659542362606640E+00, -0.9586680254064150E+00}, + { 0.2659542365158540E+00, 0.9586680253225770E+00, -0.1011136059261060E+00}, + { 0.9586680254495000E+00, -0.1011136059091010E+00, 0.2659542360648080E+00}, + {-0.9586680254786000E+00, 0.1011136057894970E+00, 0.2659542360053860E+00}, + {-0.7844318144170850E+00, 0.2843190250072290E+00, 0.5512072392025160E+00}, + { 0.2843190248228480E+00, 0.5512072393207090E+00, -0.7844318144008620E+00}, + { 0.7844318144434220E+00, -0.2843190248881310E+00, 0.5512072392264670E+00}, + { 0.5512072394346770E+00, -0.7844318142918880E+00, 0.2843190249025560E+00}, + {-0.2843190246401610E+00, 0.5512072393475040E+00, 0.7844318144482491E+00}, + { 0.5512072394083570E+00, 0.7844318144009980E+00, -0.2843190246525460E+00}, + {-0.2843190247149400E+00, -0.5512072391601370E+00, -0.7844318145528040E+00}, + {-0.5512072394176490E+00, 0.7844318144267430E+00, 0.2843190245635030E+00}, + { 0.2843190244771060E+00, -0.5512072393940670E+00, 0.7844318144746290E+00}, + {-0.5512072392271640E+00, -0.7844318145108320E+00, -0.2843190247007970E+00}, + {-0.7844318146549000E+00, -0.2843190247577290E+00, -0.5512072389927720E+00}, + { 0.7844318145421389E+00, 0.2843190246898840E+00, -0.5512072391882400E+00}, + { 0.1666638785351180E+00, 0.9794687788666500E+00, 0.1134198519532850E+00}, + { 0.9794687788923619E+00, 0.1134198520112480E+00, 0.1666638783445640E+00}, + {-0.1666638783223350E+00, -0.9794687788772219E+00, 0.1134198521746590E+00}, + { 0.1134198518526030E+00, 0.1666638784650920E+00, 0.9794687788902240E+00}, + {-0.9794687789080510E+00, 0.1134198522332290E+00, -0.1666638781012970E+00}, + { 0.1134198520235320E+00, -0.1666638782131650E+00, -0.9794687789132980E+00}, + {-0.9794687788914180E+00, -0.1134198520887550E+00, 0.1666638782973680E+00}, + {-0.1134198519422990E+00, -0.1666638783837850E+00, 0.9794687788936730E+00}, + { 0.9794687788877920E+00, -0.1134198522526510E+00, -0.1666638782071420E+00}, + {-0.1134198518873330E+00, 0.1666638784200610E+00, -0.9794687788938650E+00}, + { 0.1666638785133120E+00, -0.9794687788588400E+00, -0.1134198520527750E+00}, + {-0.1666638785259920E+00, 0.9794687788524030E+00, -0.1134198520897270E+00}, + { 0.9035426353908700E+00, 0.9900269067959901E-01, 0.4169042735078650E+00}, + { 0.9900269051118001E-01, 0.4169042737536920E+00, 0.9035426352958970E+00}, + {-0.9035426353835330E+00, -0.9900269064792300E-01, 0.4169042735312880E+00}, + { 0.4169042739582500E+00, 0.9035426351937680E+00, 0.9900269058185000E-01}, + {-0.9900269041493300E-01, 0.4169042736997320E+00, -0.9035426353313410E+00}, + { 0.4169042738439640E+00, -0.9035426352375170E+00, -0.9900269066384500E-01}, + {-0.9900269046419199E-01, -0.4169042739372540E+00, 0.9035426352163480E+00}, + {-0.4169042742060360E+00, -0.9035426351101470E+00, 0.9900269030157501E-01}, + { 0.9900269012804400E-01, -0.4169042740643800E+00, -0.9035426351945230E+00}, + {-0.4169042741137440E+00, 0.9035426351313860E+00, -0.9900269049639200E-01}, + { 0.9035426352792750E+00, -0.9900269046710200E-01, -0.4169042738001830E+00}, + {-0.9035426352343990E+00, 0.9900269024582900E-01, -0.4169042739499880E+00}, + { 0.2787624045360920E+00, 0.3493121855370630E+00, -0.8945795206981750E+00}, + { 0.3493121855860560E+00, -0.8945795206085150E+00, 0.2787624047624310E+00}, + {-0.2787624045405250E+00, -0.3493121855034730E+00, -0.8945795207099100E+00}, + {-0.8945795207341440E+00, 0.2787624047279170E+00, 0.3493121852918660E+00}, + {-0.3493121854667010E+00, -0.8945795206777230E+00, -0.2787624046898960E+00}, + {-0.8945795207888640E+00, -0.2787624046586770E+00, -0.3493121852069840E+00}, + {-0.3493121855510410E+00, 0.8945795206827980E+00, 0.2787624045679230E+00}, + { 0.8945795207852190E+00, -0.2787624046804690E+00, 0.3493121851989290E+00}, + { 0.3493121855496230E+00, 0.8945795206792300E+00, -0.2787624045811490E+00}, + { 0.8945795207818050E+00, 0.2787624045559080E+00, -0.3493121853070750E+00}, + { 0.2787624044379500E+00, -0.3493121855065000E+00, 0.8945795207406920E+00}, + {-0.2787624044432590E+00, 0.3493121854287870E+00, 0.8945795207693820E+00}, + { 0.5558962301794150E+00, -0.6768332117366710E+00, 0.4825724658147600E+00}, + {-0.6768332116815670E+00, 0.4825724660401160E+00, 0.5558962300508760E+00}, + {-0.5558962303148920E+00, 0.6768332115229870E+00, 0.4825724659584010E+00}, + { 0.4825724659102830E+00, 0.5558962301646720E+00, -0.6768332116806730E+00}, + { 0.6768332114576920E+00, 0.4825724660928950E+00, -0.5558962302776390E+00}, + { 0.4825724659029810E+00, -0.5558962303679090E+00, 0.6768332115189570E+00}, + { 0.6768332116355920E+00, -0.4825724660719810E+00, 0.5558962300791910E+00}, + {-0.4825724661505860E+00, -0.5558962302300841E+00, -0.6768332114556160E+00}, + {-0.6768332114382860E+00, -0.4825724663277370E+00, -0.5558962300974000E+00}, + {-0.4825724659723730E+00, 0.5558962302677700E+00, 0.6768332115517270E+00}, + { 0.5558962301926910E+00, 0.6768332115894530E+00, -0.4825724660059490E+00}, + {-0.5558962301943380E+00, -0.6768332114555370E+00, -0.4825724661918750E+00} +}});