Commit 06b36691 authored by Joachim Pouderoux's avatar Joachim Pouderoux Committed by Kitware Robot

Merge topic 'FixExtractFunctionalBagPlots'

582187de Fix the functional bag plot extraction filter.
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: Mathieu Westphal's avatarMathieu Westphal <mathieu.westphal@kitware.com>
Merge-request: !1110
parents 1749f98d 582187de
......@@ -1153,6 +1153,7 @@ vtkPlot * vtkChartXY::AddPlot(int type)
case FUNCTIONALBAG:
{
vtkPlotFunctionalBag *bag = vtkPlotFunctionalBag::New();
bag->GetPen()->SetColor(color.GetData());
bag->GetBrush()->SetColor(color.GetData());
plot = bag;
break;
......
......@@ -164,6 +164,7 @@ vtkStdString vtkPlotHistogram2D::GetTooltipLabel(const vtkVector2d &plotPos,
double bounds[4];
this->GetBounds(bounds);
int width = this->Input->GetExtent()[1] - this->Input->GetExtent()[0] + 1;
int height = this->Input->GetExtent()[3] - this->Input->GetExtent()[2] + 1;
int pointX = seriesIndex % width;
int pointY = seriesIndex / width;
......@@ -198,9 +199,12 @@ vtkStdString vtkPlotHistogram2D::GetTooltipLabel(const vtkVector2d &plotPos,
}
break;
case 'v':
tooltipLabel += this->GetNumber(this->Input->GetScalarComponentAsDouble(
pointX, pointY, 0, 0),
NULL);
if (pointX >= 0 && pointX < width && pointY >= 0 && pointY < height)
{
tooltipLabel +=
this->GetNumber(this->Input->GetScalarComponentAsDouble(
pointX, pointY, 0, 0), NULL);
}
break;
default: // If no match, insert the entire format tag
tooltipLabel += "%";
......
......@@ -230,6 +230,11 @@ bool vtkTransposeTableInternal::TransposeTable(vtkTable* inTable,
}
}
// Compute the number of chars needed to write the largest column id
std::stringstream ss;
ss << firstCol->GetNumberOfTuples();
std::size_t maxBLen = ss.str().length();
// Set id column on transposed table
firstCol = this->InTable->GetColumn(0);
for (int r = 0; r < firstCol->GetNumberOfComponents() *
......@@ -242,9 +247,12 @@ bool vtkTransposeTableInternal::TransposeTable(vtkTable* inTable,
}
else
{
std::stringstream ss;
ss << r;
destColumn->SetName(ss.str().c_str());
// Set the column name to the (padded) row id.
// We padd ids with 0 to avoid downstream dictionnary sort issues.
std::stringstream ss2;
ss2 << std::setw(maxBLen) << std::setfill('0');
ss2 << r;
destColumn->SetName(ss2.str().c_str());
}
}
......
......@@ -25,7 +25,11 @@
#include "vtkTestErrorObserver.h"
#include "vtkExecutive.h"
#include <algorithm>
#include <sstream>
#include <vector>
#include <vtksys/SystemTools.hxx>
#define CHECK_ERROR_MSG(observer, msg) \
{ \
......@@ -117,8 +121,37 @@ int TestExtractFunctionalBagPlot(int , char * [])
vtkNew<vtkTable> inTableDensity;
inTableDensity->AddColumn(density.GetPointer());
inTableDensity->AddColumn(varName.GetPointer());
std::vector<double> sortedDensities(
densities, densities + sizeof(densities) / sizeof(double));
double totalSumOfDensities = 0.;
for (std::size_t i = 0; i < sortedDensities.size(); i++)
{
totalSumOfDensities += sortedDensities[i];
}
std::sort(sortedDensities.begin(), sortedDensities.end());
double sumOfDensities = 0.;
double sumForP50 = totalSumOfDensities * 0.5;
double sumForP95 = totalSumOfDensities * ((100. - 95.) / 100.);
double p50 = 0.;
double p95 = 0.;
for (std::size_t i = 0; i < sortedDensities.size(); i++)
{
sumOfDensities += sortedDensities[i];
if (sumOfDensities >= sumForP50 && p50 == 0.)
{
p50 = sortedDensities[i];
}
if (sumOfDensities >= sumForP95 && p95 == 0.)
{
p95 = sortedDensities[i];
}
}
vtkNew<vtkExtractFunctionalBagPlot> ebp;
ebp->SetDensityForP50(p50);
ebp->SetDensityForPUser(p95);
ebp->SetPUser(95);
vtkNew<vtkTest::ErrorObserver> errorObserver1;
// First verify that absence of input does not cause trouble
......@@ -135,9 +168,17 @@ int TestExtractFunctionalBagPlot(int , char * [])
ebp->Update();
vtkTable* outBPTable = ebp->GetOutput();
vtkDoubleArray* q3Points =
vtkDoubleArray::SafeDownCast(outBPTable->GetColumnByName("Q3Points"));
vtkDoubleArray* q3Points = 0;
for (vtkIdType i = 0; i < outBPTable->GetNumberOfColumns(); i++)
{
const char* colName = outBPTable->GetColumnName(i);
if (vtksys::SystemTools::StringStartsWith(colName, "Q3Points"))
{
q3Points =
vtkDoubleArray::SafeDownCast(outBPTable->GetColumn(i));
break;
}
}
vtkDoubleArray* q2Points =
vtkDoubleArray::SafeDownCast(outBPTable->GetColumnByName("QMedPoints"));
......@@ -170,7 +211,7 @@ int TestExtractFunctionalBagPlot(int , char * [])
double q2v[2];
q2Points->GetTuple(19, q2v);
if (q3v[0] != 95 || q3v[1] != 304 || q2v[0] != 171 || q2v[1] != 209)
if (q3v[0] != 114 || q3v[1] != 285 || q2v[0] != 171 || q2v[1] != 209)
{
outBPTable->Dump();
cout << "## Failure: bad values found in Q3Points or QMedPoints" << endl;
......
......@@ -26,6 +26,7 @@
#include <algorithm>
#include <set>
#include <sstream>
#include <vector>
vtkStandardNewMacro(vtkExtractFunctionalBagPlot);
......@@ -34,6 +35,9 @@ vtkStandardNewMacro(vtkExtractFunctionalBagPlot);
vtkExtractFunctionalBagPlot::vtkExtractFunctionalBagPlot()
{
this->SetNumberOfInputPorts(2);
this->DensityForP50 = 0;
this->DensityForPUser = 0.;
this->PUser = 95;
}
//-----------------------------------------------------------------------------
......@@ -102,49 +106,35 @@ int vtkExtractFunctionalBagPlot::RequestData(vtkInformation* /*request*/,
vtkIdType nbPoints = varName->GetNumberOfValues();
// Fetch and sort arrays according their density
std::vector<DensityVal> varNames;
varNames.reserve(nbPoints);
for (int i = 0; i < nbPoints; i++)
{
varNames.push_back(DensityVal(density->GetValue(i),
inTable->GetColumnByName(varName->GetValue(i))));
}
std::sort(varNames.begin(), varNames.end());
std::vector<vtkAbstractArray*> medianLines;
std::vector<vtkAbstractArray*> q3Lines;
std::set<vtkAbstractArray*> outliersSeries;
// Compute total density sum
double densitySum = 0.0;
for (vtkIdType i = 0; i < nbPoints; i++)
{
densitySum += density->GetTuple1(i);
}
std::set<vtkIdType> outliersSeries;
double sum = 0.0;
for (vtkIdType i = 0; i < nbPoints; i++)
{
sum += varNames[i].Density;
if (sum < 0.5 * densitySum)
double d = density->GetValue(i);
vtkAbstractArray* c = inTable->GetColumnByName(varName->GetValue(i));
if (d < this->DensityForPUser)
{
medianLines.push_back(varNames[i].Array);
}
if (sum < 0.99 * densitySum)
{
q3Lines.push_back(varNames[i].Array);
outliersSeries.insert(i);
}
else
{
outliersSeries.insert(varNames[i].Array);
if (d > this->DensityForP50)
{
medianLines.push_back(c);
}
else
{
q3Lines.push_back(c);
}
}
}
vtkIdType nbRows = inTable->GetNumberOfRows();
vtkIdType nbCols = inTable->GetNumberOfColumns();
// Generate the median line
// Generate the median curve with median values for every sample.
vtkNew<vtkDoubleArray> qMedPoints;
qMedPoints->SetName("QMedianLine");
qMedPoints->SetNumberOfComponents(1);
......@@ -163,8 +153,10 @@ int vtkExtractFunctionalBagPlot::RequestData(vtkInformation* /*request*/,
}
// Generate the quad strip arrays
std::ostringstream ss;
ss << "Q3Points" << this->PUser;
vtkNew<vtkDoubleArray> q3Points;
q3Points->SetName("Q3Points");
q3Points->SetName(ss.str().c_str());
q3Points->SetNumberOfComponents(2);
q3Points->SetNumberOfTuples(nbRows);
......@@ -202,7 +194,7 @@ int vtkExtractFunctionalBagPlot::RequestData(vtkInformation* /*request*/,
for (vtkIdType i = 0; i < inNbColumns; i++)
{
vtkAbstractArray* arr = inTable->GetColumn(i);
if (outliersSeries.find(arr) != outliersSeries.end())
if (outliersSeries.find(i) != outliersSeries.end())
{
vtkAbstractArray* arrCopy = arr->NewInstance();
arrCopy->DeepCopy(arr);
......@@ -218,8 +210,14 @@ int vtkExtractFunctionalBagPlot::RequestData(vtkInformation* /*request*/,
}
// Then add the 2 "bag" columns into the output table
outTable->AddColumn(q3Points.GetPointer());
outTable->AddColumn(q2Points.GetPointer());
if (q3Lines.size() > 0)
{
outTable->AddColumn(q3Points.GetPointer());
}
if (medianLines.size() > 0)
{
outTable->AddColumn(q2Points.GetPointer());
}
outTable->AddColumn(qMedPoints.GetPointer());
return 1;
......
......@@ -39,6 +39,14 @@ public:
vtkTypeMacro(vtkExtractFunctionalBagPlot, vtkTableAlgorithm);
virtual void PrintSelf(ostream& os, vtkIndent indent);
// Density value for the median quartile.
vtkSetMacro(DensityForP50, double);
// Description:
// Density value for the user defined quartile.
vtkSetMacro(DensityForPUser, double);
vtkSetMacro(PUser, int);
protected:
vtkExtractFunctionalBagPlot();
virtual ~vtkExtractFunctionalBagPlot();
......@@ -47,6 +55,12 @@ protected:
vtkInformationVector**,
vtkInformationVector*);
char *P50String;
char *PUserString;
double DensityForP50;
double DensityForPUser;
int PUser;
private:
vtkExtractFunctionalBagPlot( const vtkExtractFunctionalBagPlot& ); // Not implemented.
void operator = ( const vtkExtractFunctionalBagPlot& ); // Not implemented.
......
......@@ -34,11 +34,9 @@ vtkStandardNewMacro(vtkHighestDensityRegionsStatistics);
// ----------------------------------------------------------------------
vtkHighestDensityRegionsStatistics::vtkHighestDensityRegionsStatistics()
{
this->SmoothHC1[0] = 0.;
// Initialize H smooth matrix to Identity.
this->SmoothHC1[0] = 1.0;
this->SmoothHC1[1] = 0.0;
this->SmoothHC2[0] = 0.0;
this->SmoothHC2[1] = 1.0;
this->SetSigma(1.0);
// At the construction, no columns pair are requested yet
this->NumberOfRequestedColumnsPair = 0;
......@@ -55,7 +53,7 @@ void vtkHighestDensityRegionsStatistics::PrintSelf(ostream& os,
{
this->Superclass::PrintSelf(os, indent);
os << indent << "Smooth matrix: " <<
os << indent << "Sigma matrix: " <<
this->SmoothHC1[0] << ", " <<
this->SmoothHC1[1] << ", " <<
this->SmoothHC2[0] << ", " <<
......@@ -63,23 +61,43 @@ void vtkHighestDensityRegionsStatistics::PrintSelf(ostream& os,
}
// ----------------------------------------------------------------------
void vtkHighestDensityRegionsStatistics::SetSigma(double sigma)
void vtkHighestDensityRegionsStatistics::SetSigmaMatrix(
double s11, double s12, double s21, double s22)
{
if (this->SmoothHC1[0] == sigma &&
this->SmoothHC1[1] == 0.0 &&
this->SmoothHC2[0] == 0.0 &&
this->SmoothHC2[1] == sigma)
if (this->SmoothHC1[0] == s11 && this->SmoothHC1[1] == s12 &&
this->SmoothHC2[0] == s21 && this->SmoothHC2[1] == s22)
{
return;
}
// Force H matrix to be equal to sigma * Identity.
this->SmoothHC1[0] = sigma;
this->SmoothHC1[1] = 0.0;
this->SmoothHC2[0] = 0.0;
this->SmoothHC2[1] = sigma;
this->SmoothHC1[0] = s11;
this->SmoothHC1[1] = s12;
this->SmoothHC2[0] = s21;
this->SmoothHC2[1] = s22;
this->Determinant =
vtkMath::Determinant2x2(this->SmoothHC1, this->SmoothHC2);
double invDet = 0.;
if (this->Determinant != 0.)
{
invDet = 1.0 / this->Determinant;
}
// We compute and store the inverse of the smoothing matrix
this->InvSigmaC1[0] = +invDet * this->SmoothHC2[1];
this->InvSigmaC1[1] = -invDet * this->SmoothHC1[1];
this->InvSigmaC2[0] = -invDet * this->SmoothHC2[0];
this->InvSigmaC2[1] = +invDet * this->SmoothHC1[0];
this->Modified();
}
// ----------------------------------------------------------------------
void vtkHighestDensityRegionsStatistics::SetSigma(double sigma)
{
this->SetSigmaMatrix(sigma * sigma, 0, 0, sigma * sigma);
}
// ----------------------------------------------------------------------
void vtkHighestDensityRegionsStatistics::Learn(vtkTable* inData,
vtkTable* vtkNotUsed(inParameters),
......@@ -229,16 +247,10 @@ double vtkHighestDensityRegionsStatistics
// Sum all gaussian kernel
for (vtkIdType j = 0; j < nbObservations; j++)
{
inObs->GetTuple(j, currentXj);
const double deltaX = currentXi[0] - currentXj[0];
const double deltaY = currentXi[1] - currentXj[1];
// Avoid case where point is compared to itself
if (deltaX == 0. && deltaY == 0.)
{
continue;
}
hdr += this->ComputeSmoothGaussianKernel(
inObs->GetNumberOfComponents(),
deltaX, deltaY);
......@@ -253,53 +265,12 @@ double vtkHighestDensityRegionsStatistics
// ----------------------------------------------------------------------
double vtkHighestDensityRegionsStatistics::ComputeSmoothGaussianKernel(
int dimension, double khx, double khy)
int vtkNotUsed(dimension), double khx, double khy)
{
double HDeterminant =
vtkMath::Determinant2x2(this->SmoothHC1, this->SmoothHC2);
if (HDeterminant > 0.0)
{
HDeterminant = 1.0 / sqrt(HDeterminant);
}
// We need to multiply the input vector by the smooth square root of
// H matrix parameter: sqrt(H) * [khx, khy] -> random vector of the
// standard gaussian input.
// If a H coefficient is equal to 0.0. we don't compute its sqrt to avoid
// domain error.
double SHC10 = 0.0;
double SHC11 = 0.0;
double SHC20 = 0.0;
double SHC21 = 0.0;
if (this->SmoothHC1[0] != 0.0)
{
SHC10 = 1.0 / sqrt(this->SmoothHC1[0]);
}
if (this->SmoothHC1[1] != 0.0)
{
SHC11 = 1.0 / sqrt(this->SmoothHC1[1]);
}
if (this->SmoothHC2[0] != 0.0)
{
SHC20 = 1.0 / sqrt(this->SmoothHC2[0]);
}
if (this->SmoothHC2[1] != 0.0)
{
SHC21 = 1.0 / sqrt(this->SmoothHC2[1]);
}
// Call the standard gaussian kernel with the new random vector.
return HDeterminant *
this->ComputeStandardGaussianKernel(dimension,
SHC10 * khx + SHC11 * khy,
SHC20 * khx + SHC21 * khy);
}
double d =
khx * (this->InvSigmaC1[0] * khx + this->InvSigmaC2[0] * khy) +
khy * (this->InvSigmaC1[1] * khx + this->InvSigmaC2[1] * khy);
// ----------------------------------------------------------------------
double vtkHighestDensityRegionsStatistics::ComputeStandardGaussianKernel(
int vtkNotUsed(dimension), double kx, double ky)
{
return exp(-(kx * kx + ky * ky) / 2.0) / (2.0 * vtkMath::Pi());
return (exp(-d * 0.5)) / (2.0 * vtkMath::Pi() * this->Determinant);
}
......@@ -56,18 +56,12 @@ public:
vtkMultiBlockDataSet*) { return; }
// Description:
// H is a positive matrix that defines the smooth direction.
// In a classical HDR, we don't set a specific smooth direction for the
// H matrix parameter (SmoothHC1, SmoothHC2). That mean H will be in a
// diagonal form and equal to sigma * Id.
// Set the width of the gaussian kernel.
void SetSigma(double sigma);
// Description:
// Get Smooth H matrix parameter of the HDR.
vtkGetVectorMacro(SmoothHC1, double, 2);
vtkSetVectorMacro(SmoothHC1, double, 2);
vtkGetVectorMacro(SmoothHC2, double, 2);
vtkSetVectorMacro(SmoothHC2, double, 2);
// Set the gaussian kernel matrix.
void SetSigmaMatrix(double s11, double s12, double s21, double s22);
// Description:
// Fill outDensity with density vector that is computed from
......@@ -85,6 +79,7 @@ public:
// Look ComputeSmoothGaussianKernel for KH kernel definition.
double ComputeHDR(vtkDataArray *inObs, vtkDataArray* inPOI,
vtkDataArray *outDensity);
protected:
vtkHighestDensityRegionsStatistics();
~vtkHighestDensityRegionsStatistics();
......@@ -125,6 +120,9 @@ protected:
// for the Gaussian kernel.
double SmoothHC1[2];
double SmoothHC2[2];
double InvSigmaC1[2];
double InvSigmaC2[2];
double Determinant;
// Description:
// Store the number of requested columns pair computed by learn method.
......@@ -138,12 +136,6 @@ private :
// Look ComputeStandardGaussianKernel for the K kernel definition.
double ComputeSmoothGaussianKernel(int dimension, double khx, double khy);
// Description:
// Helper that returns a standard gaussian kernel of a vector of dimension two,
// using its coordinates. For X = [kx, ky],
// K(X) = ( 1 / 2 * PI) * exp(-sqrt<X,X>).
double ComputeStandardGaussianKernel(int dimension, double kx, double ky);
private:
vtkHighestDensityRegionsStatistics(const vtkHighestDensityRegionsStatistics&); // Not implemented
void operator = (const vtkHighestDensityRegionsStatistics&); // Not implemented
......
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