Commit 88e997f9 authored by hrchilds's avatar hrchilds

Update from November 21, 2005

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@574 18c085ea-50e0-402c-830e-de6fd14e8384
parent 080c1e46
#include <avtDistanceToBestFitLineFilter.h>
#include <vtkDataArray.h>
#include <avtDataTree.h>
#include <avtParallel.h>
#include <math.h>
#define N_SUM 0
#define X_SUM 1
#define Y_SUM 2
#define XY_SUM 3
#define X2_SUM 4
#define N_CALC_VALUES 5
// ****************************************************************************
// Method: avtDistanceToBestFitLineFilter::avtDistanceToBestFitLineFilter
//
// Purpose:
// Constructor for the avtDistanceToBestFitLineFilter class.
//
// Arguments:
// v : If true then difference vertically - otherwise calculate perpendicular
// distance.
//
// Programmer: Brad Whitlock
// Creation: Fri Nov 18 16:15:57 PST 2005
//
// Modifications:
//
// ****************************************************************************
avtDistanceToBestFitLineFilter::avtDistanceToBestFitLineFilter(bool v) :
avtBinaryMathFilter()
{
verticalDifference = v;
pass = 1;
sums[N_SUM] = 0.;
sums[X_SUM] = 0.;
sums[Y_SUM] = 0.;
sums[XY_SUM] = 0.;
sums[X2_SUM] = 0.;
}
// ****************************************************************************
// Method: avtDistanceToBestFitLineFilter::~avtDistanceToBestFitLineFilter
//
// Purpose:
// Destructor for the avtDistanceToBestFitLineFilter class.
//
// Programmer: Brad Whitlock
// Creation: Fri Nov 18 16:16:49 PST 2005
//
// Modifications:
//
// ****************************************************************************
avtDistanceToBestFitLineFilter::~avtDistanceToBestFitLineFilter()
{
}
// ****************************************************************************
// Method: avtDistanceToBestFitLineFilter::PreExecute
//
// Purpose:
// Sets all of the sums to zero.
//
// Programmer: Brad Whitlock
// Creation: Fri Nov 18 16:17:07 PST 2005
//
// Modifications:
//
// ****************************************************************************
void
avtDistanceToBestFitLineFilter::PreExecute(void)
{
sums[N_SUM] = 0.;
sums[X_SUM] = 0.;
sums[Y_SUM] = 0.;
sums[XY_SUM] = 0.;
sums[X2_SUM] = 0.;
}
// ****************************************************************************
// Method: avtDistanceToBestFitLineFilter::Execute
//
// Purpose:
// Performs a multi-pass execute. In pass 1 we calculate sums for the best
// fit line. In pass 2, we calculate the difference between the variables
// and the best fit line.
//
// Programmer: Brad Whitlock
// Creation: Fri Nov 18 16:17:24 PST 2005
//
// Modifications:
//
// ****************************************************************************
void
avtDistanceToBestFitLineFilter::Execute(void)
{
//
// Sum the values required to do the best fit line.
//
pass = 1;
avtDataTree_p tree = GetInputDataTree();
totalNodes = 2 * tree->GetNumberOfLeaves();
avtDataTree_p newTree = avtDataTreeStreamer::Execute(tree);
newTree = 0;
// Sum the array values over all processors, making sure each processor
// gets the results
double d[N_CALC_VALUES];
SumDoubleArrayAcrossAllProcessors(sums, d, N_CALC_VALUES);
for(int i = 0; i < N_CALC_VALUES; ++i)
sums[i] = d[i];
//
// Make it perform the expression.
//
pass = 2;
avtDataTree_p newTree2 = avtDataTreeStreamer::Execute(tree);
SetOutputDataTree(newTree2);
}
// ****************************************************************************
// Method: avtDistanceToBestFitLineFilter::DoOperation
//
// Purpose:
// Does the work of calculating the expression.
//
// Arguments:
// in1 : The data array containing the X coordinate.
// in2 : The data array containing the Y coordinate.
// out : The resulting data.
// ncomps : The number of components.
// ntuples : The number of tuples.
//
// Note: Note that we only populate the out data array in pass 2. The
// results of pass 1 get thrown away.
//
// Programmer: Brad Whitlock
// Creation: Fri Nov 18 16:18:09 PST 2005
//
// Modifications:
//
// ****************************************************************************
void
avtDistanceToBestFitLineFilter::DoOperation(vtkDataArray *in1,
vtkDataArray *in2, vtkDataArray *out, int ncomps, int ntuples)
{
if(pass == 1)
{
// Sum up the values required to calculate the best fit line.
sums[N_SUM] += double(ntuples);
for(vtkIdType i = 0; i < ntuples; ++i)
{
float x = in1->GetTuple1(i);
float y = in2->GetTuple1(i);
sums[X_SUM] += double(x);
sums[Y_SUM] += double(y);
sums[XY_SUM] += double(x * y);
sums[X2_SUM] += double(x * x);
}
}
else if(pass == 2)
{
double dY = (sums[N_SUM] * sums[XY_SUM] - sums[X_SUM] * sums[Y_SUM]);
double dX = (sums[N_SUM] * sums[X2_SUM] - sums[X_SUM] * sums[X_SUM]);
double m, b;
if(dX == 0.)
{
double xLine = sums[X_SUM] / sums[N_SUM];
for(vtkIdType i = 0; i < ntuples; ++i)
{
float x = in1->GetTuple1(i);
out->SetTuple1(i, x - xLine);
}
}
else
{
m = dY / dX;
b = (sums[Y_SUM] - m * sums[X_SUM]) / sums[N_SUM];
if(verticalDifference)
{
for(vtkIdType i = 0; i < ntuples; ++i)
{
float x = in1->GetTuple1(i);
float y = in2->GetTuple1(i);
float yLine = m * x + b;
out->SetTuple1(i, y - yLine);
}
}
else // perpendicular distance.
{
for(vtkIdType i = 0; i < ntuples; ++i)
{
float x0 = in1->GetTuple1(i);
float y0 = in2->GetTuple1(i);
// We know line eq. Y = MX + B
// Reorganized: ax + by + c = 0 where a=M, b=-1, c=B
//
// Dist from point to line:
// |a*x0 + b*y0 + c| / sqrt(a^2 + b^2)
//
// Subst a,b,c
// |M*x0 -1*y0 + B| / sqrt(M^2 + 1)
float d = fabs(m * x0 - y0 + b) / sqrt(m*m + 1);
out->SetTuple1(i, d);
}
}
}
}
}
// ************************************************************************* //
// avtDistanceToBestFitLineFilter.h //
// ************************************************************************* //
#ifndef AVT_DISTANCE_TO_BEST_FIT_LINE_FILTER_H
#define AVT_DISTANCE_TO_BEST_FIT_LINE_FILTER_H
#include <avtBinaryMathFilter.h>
class vtkDataArray;
// ****************************************************************************
// Class: avtDistanceToBestFitLineFilter
//
// Purpose:
// Computes the best fit line for a 2-tuple of input variables and then,
// as a second pass, calculates each data value's distance from that line.
//
// Programmer: Brad Whitlock
// Creation: Fri Nov 18 14:23:50 PST 2005
//
// Modifications:
//
// ****************************************************************************
class EXPRESSION_API avtDistanceToBestFitLineFilter
: public avtBinaryMathFilter
{
public:
avtDistanceToBestFitLineFilter(bool v);
virtual ~avtDistanceToBestFitLineFilter();
virtual const char *GetType(void)
{ return "avtDistanceToBestFitLineFilter"; }
virtual const char *GetDescription(void)
{ return "Distance to best fit line"; }
protected:
bool verticalDifference;
int pass;
double sums[6];
virtual void PreExecute(void);
virtual void Execute(void);
virtual void DoOperation(vtkDataArray *in1, vtkDataArray *in2,
vtkDataArray *out, int ncomps, int ntuples);
virtual avtVarType GetVariableType(void) { return AVT_SCALAR_VAR; };
};
#endif
......@@ -114,6 +114,9 @@
# Hank Childs, Wed Sep 21 17:16:30 PDT 2005
# Added avtExternalNodeExpression, avtSurfaceNormalExpression, avtEdgeLength.
#
# Brad Whitlock, Fri Nov 18 15:45:58 PST 2005
# Added avtDistanceToBestFitLineFilter.
#
##############################################################################
@SET_MAKE@
......@@ -164,6 +167,7 @@ General_src= \
General/avtDataIdFilter.C \
General/avtDegreeFilter.C \
General/avtDivergenceFilter.C \
General/avtDistanceToBestFitLineFilter.C \
General/avtExpressionComponentMacro.C \
General/avtExternalNodeExpression.C \
General/avtGradientFilter.C \
......
......@@ -14,6 +14,7 @@
#include <avtRoundFilter.h>
#include <avtSinFilter.h>
#include <avtCosFilter.h>
#include <avtDistanceToBestFitLineFilter.h>
#include <avtRandomFilter.h>
#include <avtArctanFilter.h>
#include <avtArcsinFilter.h>
......@@ -413,6 +414,9 @@ avtVectorExpr::CreateFilters(ExprPipelineState *state)
// Added external_node, surface normals, min_side_volume, max_side_volume,
// min_edge_length, and max_edge_length.
//
// Brad Whitlock, Fri Nov 18 15:59:27 PST 2005
// Added distance_to_best_fit_line.
//
// ****************************************************************************
void
avtFunctionExpr::CreateFilters(ExprPipelineState *state)
......@@ -715,6 +719,16 @@ avtFunctionExpr::CreateFilters(ExprPipelineState *state)
ecm->SetMacro("polar", 2);
f = ecm;
}
else if (functionName == "distance_to_best_fit_line")
f = new avtDistanceToBestFitLineFilter(true);
else if (functionName == "distance_to_best_fit_line2")
f = new avtDistanceToBestFitLineFilter(false);
else if (functionName == "polar_phi")
{
avtExpressionComponentMacro *ecm = new avtExpressionComponentMacro;
ecm->SetMacro("polar", 2);
f = ecm;
}
else
{
string error =
......
......@@ -110,6 +110,9 @@
# Kathleen Bonnell, Tue Nov 8 10:45:43 PST 2005
# Added TrajectoryByNode/Zone.
#
# Brad Whitlock, Thu Nov 17 10:15:40 PDT 2005
# Added Best Fit Line.
#
##############################################################################
@SET_MAKE@
......@@ -150,6 +153,7 @@ Queries_src= \
Queries/avtActualDataNumNodesQuery.C \
Queries/avtActualDataNumZonesQuery.C \
Queries/avtAreaBetweenCurvesQuery.C \
Queries/avtBestFitLineQuery.C \
Queries/avtCentroidQuery.C \
Queries/avtCompactnessQuery.C \
Queries/avtCycleQuery.C \
......
......@@ -7,6 +7,7 @@
#include <avtActualDataNumNodesQuery.h>
#include <avtActualDataNumZonesQuery.h>
#include <avtAreaBetweenCurvesQuery.h>
#include <avtBestFitLineQuery.h>
#include <avtCentroidQuery.h>
#include <avtCompactnessQuery.h>
#include <avtCycleQuery.h>
......@@ -151,6 +152,9 @@ avtQueryFactory::Instance()
// Kathleen Bonnell, Tue Nov 8 10:45:43 PST 2005
// Added TrajectoryByNode/Zone queries.
//
// Brad Whitlock, Thu Nov 17 10:18:41 PDT 2005
// Added Best Fit Line.
//
// ****************************************************************************
......@@ -325,6 +329,11 @@ avtQueryFactory::CreateQuery(const QueryAttributes *qa)
{
query = new avtTrajectoryByNode();
}
else if (qname == "Best Fit Line")
{
query = new avtBestFitLineQuery();
}
return query;
}
......
// ************************************************************************* //
// avtBestFitLineQuery.C //
// ************************************************************************* //
#include <avtBestFitLineQuery.h>
#include <vtkCellData.h>
#include <vtkDataSet.h>
#include <vtkIdList.h>
#include <vtkPointData.h>
#include <vtkUnsignedCharArray.h>
#include <avtCallback.h>
#include <avtParallel.h>
#include <DebugStream.h>
#include <InvalidVariableException.h>
#include <snprintf.h>
using std::string;
#define N_SUM 0
#define X_SUM 1
#define Y_SUM 2
#define XY_SUM 3
#define X2_SUM 4
#define Y2_SUM 5
// ****************************************************************************
// Method: avtBestFitLineQuery constructor
//
// Programmer: Brad Whitlock
// Creation: Wed Nov 16 14:38:19 PST 2005
//
// ****************************************************************************
avtBestFitLineQuery::avtBestFitLineQuery()
{
sums[N_SUM] = 0.;
sums[X_SUM] = 0.;
sums[Y_SUM] = 0.;
sums[XY_SUM] = 0.;
sums[X2_SUM] = 0.;
sums[Y2_SUM] = 0.;
strcpy(descriptionBuffer, "Fitting line to points");
}
// ****************************************************************************
// Method: avtBestFitLineQuery destructor
//
// Purpose:
// Defines the destructor. Note: this should not be inlined in the header
// because it causes problems for certain compilers.
//
// Programmer: Brad Whitlock
// Creation: Wed Nov 16 14:38:19 PST 2005
//
// ****************************************************************************
avtBestFitLineQuery::~avtBestFitLineQuery()
{
;
}
// ****************************************************************************
// Method: avtBestFitLineQuery::PreExecute
//
// Purpose:
// This is called before all of the domains are executed.
//
// Programmer: Brad Whitlock
// Creation: Wed Nov 16 14:38:19 PST 2005
//
// ****************************************************************************
void
avtBestFitLineQuery::PreExecute(void)
{
sums[N_SUM] = 0.;
sums[X_SUM] = 0.;
sums[Y_SUM] = 0.;
sums[XY_SUM] = 0.;
sums[X2_SUM] = 0.;
sums[Y2_SUM] = 0.;
}
// ****************************************************************************
// Method: avtBestFitLineQuery::PostExecute
//
// Purpose:
// This is called after all of the domains are executed.
//
// Programmer: Brad Whitlock
// Creation: Wed Nov 16 14:38:19 PST 2005
//
// Modifications:
//
// ****************************************************************************
void
avtBestFitLineQuery::PostExecute(void)
{
double d[6];
SumDoubleArrayAcrossAllProcessors(sums, d, 6);
debug4 << "avtBestFitLineQuery::PostExecute: \n"
<< "\tN_SUM=" << d[N_SUM] << endl
<< "\tX_SUM=" << d[X_SUM] << endl
<< "\tY_SUM=" << d[Y_SUM] << endl
<< "\tXY_SUM=" << d[XY_SUM] << endl
<< "\tX2_SUM=" << d[X2_SUM] << endl
<< "\tY2_SUM=" << d[Y2_SUM]
<< endl;
double dY = (d[N_SUM] * d[XY_SUM] - d[X_SUM] * d[Y_SUM]);
double dX = (d[N_SUM] * d[X2_SUM] - d[X_SUM] * d[X_SUM]);
double m, b, r;
std::string s;
if(dX == 0.)
{
double x = d[X_SUM] / d[N_SUM];
m = DBL_MAX;
r = 1.;
char buf[1024];
SNPRINTF(buf, 1024, "The best fit line is: X = %g with a "
"correlation coefficient of: %g", x, r);
s = std::string(buf);
}
else
{
m = dY / dX;
b = (d[Y_SUM] - m * d[X_SUM]) / d[N_SUM];
double rtop = d[N_SUM] * d[XY_SUM] - d[X_SUM] * d[Y_SUM];
double r0 = d[N_SUM] * d[X2_SUM] - d[X_SUM] * d[X_SUM];
double r1 = d[N_SUM] * d[Y2_SUM] - d[Y_SUM] * d[Y_SUM];
double rbottom = r0 * r1;
if (rbottom > 0.)
r = rtop / sqrt(rbottom);
// Create a return message.
char buf[1024];
SNPRINTF(buf, 1024, "The best fit line is: Y = %gX ", m);
s = std::string(buf);
if(b < 0.)
{
SNPRINTF(buf, 1024, "-%g ", b);
buf[1] = ' ';
s += buf;
}
else if(b > 0.)
{
SNPRINTF(buf, 1024, "+ %g ", b);
s += buf;
}
if (rbottom > 0.)
SNPRINTF(buf, 1024, "with a correlation coefficient of: %g", r);
else
SNPRINTF(buf, 1024, "with an undefined correlation coefficient.");
s += buf;
}
//
// Parent class uses this message to set the Results message
// in the Query Attributes that is sent back to the viewer.
// That is all that is required of this query.
//
SetResultMessage(s);
doubleVector retval;
retval.push_back(m);
retval.push_back(b);
retval.push_back(r);
SetResultValues(retval);
}
// ****************************************************************************
// Method: avtBestFitLineQuery::Execute
//
// Purpose:
// Processes a single domain.
//
// Programmer: Brad Whitlock
// Creation: Wed Nov 16 14:38:19 PST 2005
//
// Modifications:
//
// ****************************************************************************
void
avtBestFitLineQuery::Execute(vtkDataSet *ds, const int dom)
{
vtkIdType npts = ds->GetNumberOfPoints();
sums[N_SUM] += double(npts);
for(vtkIdType i = 0; i < npts; ++i)
{
float fptr[3];
ds->GetPoint(i, fptr);
sums[X_SUM] += double(fptr[0]);
sums[Y_SUM] += double(fptr[1]);
sums[XY_SUM] += double(fptr[0] * fptr[1]);
sums[X2_SUM] += double(fptr[0] * fptr[0]);
sums[Y2_SUM] += double(fptr[1] * fptr[1]);
}
}
// ************************************************************************* //
// avtBestFitLineQuery.h //
// ************************************************************************* //
#ifndef AVT_BEST_FIT_LINE_QUERY_H
#define AVT_BEST_FIT_LINE_QUERY_H
#include <query_exports.h>
#include <avtSummationQuery.h>
// ****************************************************************************
// Class: avtBestFitLineQuery
//
// Purpose:
// A query that will calculate the best fit line to a set of points.
//
// Programmer: Brad Whitlock
// Creation: Wed Nov 16 16:21:33 PST 2005
//
// Modifications:
//
// ****************************************************************************
class QUERY_API avtBestFitLineQuery : public avtSummationQuery
{
public:
avtBestFitLineQuery();
virtual ~avtBestFitLineQuery();
virtual const char *GetType(void)
{ return "avtBestFitLineQuery"; };
virtual void PreExecute(void);
virtual void Execute(vtkDataSet *, const int);
virtual void PostExecute(void);
protected:
double sums[6];
};
#endif
......@@ -23,6 +23,29 @@
#include <InvalidVariableException.h>
using std::string;
using std::find;
// ****************************************************************************
// Function: dune_getline
//
// Purpose:
// Implements the getline function that is available with gcc, but not
// with SGI's compiler.
//
// Programmer: Hank Childs
// Creation: November 21, 2005
//
// ****************************************************************************
bool
dune_getline(std::ifstream &ifile, std::string &str)
{
char buffer[1024];
ifile.getline(buffer, 1024, '\n');
str = buffer;
return !(ifile.eof());
}
// Note: particle and fragment are interchangable in comments/descriptions.
......@@ -50,6 +73,9 @@ using std::string;
// Combined LOADER and RESTART sections.
// Modest cleanups.
//
// Hank Childs, Mon Nov 21 11:25:48 PST 2005
// Replace getline with dune_getline.
//
// ****************************************************************************
avtDuneFileFormat::avtDuneFileFormat(const char *filename)
......@@ -72,7 +98,7 @@ avtDuneFileFormat::avtDuneFileFormat(const char *filename)
string buffer;
while (ftype == UNKNOWN && !ifile.eof()) {
getline(ifile, buffer);
dune_getline(ifile, buffer);
if (buffer.find(sinkedFragments) != string::npos) {
ftype = RESTART;
}
......@@ -81,7 +107,7 @@ avtDuneFileFormat::avtDuneFileFormat(const char *filename)
ifile.seekg((streampos) 0, ios::beg);
while (ftype == UNKNOWN) {
getline(ifile, buffer);
dune_getline(ifile, buffer);
if (buffer.find(title) != string::npos) {
ftype = TECPLOT;
}
......@@ -133,7 +159,7 @@ avtDuneFileFormat::avtDuneFileFormat(const char *filename)
switch (ftype) {
case TECPLOT:
while (getline(ifile, buffer)) {
while (dune_getline(ifile, buffer)) {
if (buffer.find(ztime) != string::npos) {
cycles.push_back(ntimes);
++ntimes;
......@@ -147,7 +173,7 @@ avtDuneFileFormat::avtDuneFileFormat(const char *filename)
int np = atoi(tokens[3].c_str());
num_particles.push_back(np);
for (int i = 0; i < np; i++) {
getline(ifile, buffer);
dune_getline(ifile, buffer);
get_tokens(buffer, empty, tokens);
ispecies = atoi(tokens[0].c_str());
max_species = MAX(max_species, ispecies