/*
* Copyright (c) 2002-2007 Samit Basu
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "HandleContour.hpp"
#include "HandleList.hpp"
#include "HandleObject.hpp"
#include "HandleAxis.hpp"
#include "IEEEFP.hpp"
#include <math.h>
HandleContour::HandleContour() {
ConstructProperties();
SetupDefaults();
}
HandleContour::~HandleContour() {
}
// Calculate the limits - should return a vector of the
// form...
// [xmin, xmax, ymin, ymax, zmin, zmax, cmin, cmax, amin, amax]
std::vector<double> HandleContour::GetLimits() {
// std::vector<double> xs(VectorPropertyLookup("xdata"));
// std::vector<double> ys(VectorPropertyLookup("ydata"));
// std::vector<double> zs(VectorPropertyLookup("zdata"));
UpdateState();
std::vector<double> limits;
bool x_initialized = false;
bool y_initialized = false;
bool z_initialized = false;
double xmin;
double xmax;
double ymin;
double ymax;
double zmin;
double zmax;
double cmin;
double cmax;
for (int i=0;i<pset.size();i++) {
// For each level.
lineset cset(pset[i]);
for (int j=0;j<cset.size();j++) {
cline aline(cset[j]);
for (int k=0;k<aline.size();k++) {
if (!x_initialized && IsFinite(aline[k].x)) {
xmin = xmax = aline[k].x;
x_initialized = true;
}
if (!y_initialized && IsFinite(aline[k].y)) {
ymin = ymax = aline[k].y;
y_initialized = true;
}
if (!z_initialized) {
zmin = zmax = zvals[0];
z_initialized = true;
}
if (x_initialized && IsFinite(aline[k].x)) {
xmin = qMin(xmin,aline[k].x);
xmax = qMax(xmax,aline[k].x);
}
if (y_initialized && IsFinite(aline[k].y)) {
ymin = qMin(ymin,aline[k].y);
ymax = qMax(ymax,aline[k].y);
}
if (z_initialized) {
zmin = qMin(zmin,zvals[i]);
zmax = qMax(zmax,zvals[i]);
}
}
}
}
if (!x_initialized) {
xmin = -1;
xmax = 1;
}
if (!y_initialized) {
ymin = -1;
ymax = 1;
}
if (!z_initialized) {
zmin = 0;
zmax = 0;
}
cmin = zmin;
cmax = zmax;
if (StringCheck("floating","off"))
zmin = zmax = 0;
limits.push_back(xmin);
limits.push_back(xmax);
limits.push_back(ymin);
limits.push_back(ymax);
limits.push_back(zmin);
limits.push_back(zmax);
limits.push_back(cmin);
limits.push_back(cmax);
limits.push_back(0);
limits.push_back(0);
return limits;
}
inline bool operator==(const contour_point& p1, const contour_point& p2) {
return ((p1.x == p2.x) && (p1.y == p2.y));
}
inline cline Reverse(const cline& src) {
cline ret;
for (int i=src.size()-1;i>=0;i--)
ret << src.at(i);
return ret;
}
inline bool Connected(const cline& current, const cline& test) {
return ((current.front() == test.front()) ||
(current.front() == test.back()) ||
(current.back() == test.front()) ||
(current.back() == test.back()));
}
inline void Join(cline& current, const cline& toadd) {
if (current.front() == toadd.front())
current = Reverse(toadd) + current;
else if (current.front() == toadd.back())
current = toadd + current;
else if (current.back() == toadd.front())
current += toadd;
else if (current.back() == toadd.back())
current += Reverse(toadd);
}
#define FOLD(x) MAP((x),row-1)
#define FNEW(x) MAP((x),row)
#define MAP(x,y) func[(y)+(x)*numy]
#define AINTER(a,b) ((val-(a))/((b)-(a)))
#define ALEFT(i,j) (((j)-1)+AINTER(FOLD((i)-1),FNEW((i)-1)))
#define TOP(i) (((i)-1)+AINTER(FNEW((i)-1),FNEW((i))))
#define BOT(i) (((i)-1)+AINTER(FOLD((i)-1),FOLD((i))))
#define RIGHT(i,j) (((j)-1)+AINTER(FOLD((i)),FNEW((i))))
#define DRAW(a,b,c,d) {allLines << (cline() << contour_point((double)a,(double)b) << contour_point((double)c,(double)d));}
lineset HandleContour::ContourCDriver(Array m, double val, Array x, Array y) {
lineset allLines;
lineset bundledLines;
m.promoteType(FM_DOUBLE);
const double *func = (const double *) m.getDataPointer();
int outcnt = 0;
int numy = m.rows();
int numx = m.columns();
for (int row=1;row<numy;row++)
for (int col=1;col<numx;col++) {
int l = 0;
if (FOLD(col) >= val) l = l + 1;
if (FOLD(col-1) >= val) l = l + 2;
if (FNEW(col) >= val) l = l + 4;
if (FNEW(col-1) >= val) l = l + 8;
switch (l) {
case 1:
case 14:
DRAW(BOT(col),row-1,col,RIGHT(col,row));
break;
case 2:
case 13:
DRAW(col-1,ALEFT(col,row),BOT(col),row-1);
break;
case 3:
case 12:
DRAW(col-1,ALEFT(col,row),col,RIGHT(col,row));
break;
case 4:
case 11:
DRAW(TOP(col),row,col,RIGHT(col,row));
break;
case 5:
case 10:
DRAW(BOT(col),row-1,TOP(col),row);
break;
case 6:
case 9:
{
double x0 = AINTER(FOLD(col-1),FOLD(col));
double x1 = AINTER(FNEW(col-1),FNEW(col));
double y0 = AINTER(FOLD(col-1),FNEW(col-1));
double y1 = AINTER(FOLD(col),FNEW(col));
double y = (x0*(y1-y0)+y0)/(1.0-(x1-x0)*(y1-y0));
double x = y*(x1-x0) + x0;
double fx1 = MAP(col-1,row-1)+x*(MAP(col,row-1)-MAP(col-1,row-1));
double fx2 = MAP(col-1,row)+x*(MAP(col,row)-MAP(col-1,row));
double f = fx1 + y*(fx2-fx1);
if (f==val) {
DRAW(BOT(col),row-1,TOP(col),row);
DRAW(col-1,ALEFT(col,row),col,RIGHT(col,row));
} else if (((f > val) && (FNEW(col) > val)) ||
((f < val) && (FNEW(col) < val))) {
DRAW(col-1,ALEFT(col,row),TOP(col),row);
DRAW(BOT(col),row-1,col,RIGHT(col,row));
} else {
DRAW(col-1,ALEFT(col,row),BOT(col),row-1);
DRAW(TOP(col),row,col,RIGHT(col,row));
}
}
break;
case 7:
case 8:
DRAW(col-1,ALEFT(col,row),TOP(col),row);
break;
}
}
// Now we link the line segments into longer lines.
int allcount = allLines.size();
while (!allLines.empty()) {
// Start a new line segment
cline current(allLines.takeAt(0));
bool lineGrown = true;
while (lineGrown) {
lineGrown = false;
int i = 0;
while (i<allLines.size()) {
if (Connected(current,allLines.at(i))) {
Join(current,allLines.takeAt(i));
lineGrown = true;
} else
i++;
}
}
bundledLines << current;
}
// Final step is to transform the lines into the
// given coordinates
const double *xp = (const double *) x.getDataPointer();
const double *yp = (const double *) y.getDataPointer();
#define X(a,b) xp[(b)+(a)*numy]
#define Y(a,b) yp[(b)+(a)*numy]
for (int i=0;i<bundledLines.size();i++)
for (int j=0;j<bundledLines[i].size();j++) {
double ndx_x = bundledLines[i][j].x;
double ndx_y = bundledLines[i][j].y;
// Compute a new x and y using these coordinates as
// bilinear interpolators into xp and yp
int indx_x = (int) ndx_x;
int indx_y = (int) ndx_y;
double eps_x = ndx_x - indx_x;
double eps_y = ndx_y - indx_y;
double xp_out_1 = X(indx_x,indx_y) +
eps_x*(X(indx_x+1,indx_y) - X(indx_x,indx_y));
double xp_out_2 = X(indx_x,indx_y+1) +
eps_x*(X(indx_x+1,indx_y+1) - X(indx_x,indx_y+1));
double xp_out = xp_out_1 + eps_y*(xp_out_2-xp_out_1);
double yp_out_1 = Y(indx_x,indx_y) +
eps_x*(Y(indx_x+1,indx_y) - Y(indx_x,indx_y));
double yp_out_2 = Y(indx_x,indx_y+1) +
eps_x*(Y(indx_x+1,indx_y+1) - Y(indx_x,indx_y+1));
double yp_out = yp_out_1 + eps_y*(yp_out_2-yp_out_1);
bundledLines[i][j].x = xp_out;
bundledLines[i][j].y = yp_out;
}
return bundledLines;
}
void HandleContour::RebuildContourMatrix() {
// Count how many total points there are
int pointcount = 0;
int linecount = 0;
for (int i=0;i<pset.size();i++) {
for (int j=0;j<pset[i].size();j++) {
linecount++;
pointcount += pset[i][j].size();
}
}
// Create the contour matrix
int outcount = pointcount+linecount;
Array out(Array::doubleMatrixConstructor(2,outcount));
double *output = (double *) out.getReadWriteDataPointer();
for (int i=0;i<pset.size();i++) {
for (int j=0;j<pset[i].size();j++) {
*output++ = zvals[i];
*output++ = pset[i][j].size();
cline bline(pset[i][j]);
for (int k=0;k<bline.size();k++) {
*output++ = bline[k].x;
*output++ = bline[k].y;
}
}
}
HPArray *hp = (HPArray*) LookupProperty("contourmatrix");
hp->Data(out);
}
Array HandleContour::GetCoordinateMatrix(std::string name, bool isXcoord) {
// Get the elevation data from the object
Array zdata(ArrayPropertyLookup("zdata"));
int zrows(zdata.rows());
int zcols(zdata.columns());
if (StringCheck(name+"mode","manual")) {
// not auto mode...
Array cdata(ArrayPropertyLookup(name));
if (cdata.isVector() &&
((isXcoord && (cdata.getLength() == zcols)) ||
(!isXcoord && (cdata.getLength() == zrows)))) {
cdata.promoteType(FM_DOUBLE);
const double *qp = (const double*) cdata.getDataPointer();
Array mat(Array::doubleMatrixConstructor(zrows,zcols));
double *dp = (double*) mat.getReadWriteDataPointer();
for (int i=0;i<zcols;i++)
for (int j=0;j<zrows;j++) {
if (isXcoord)
*dp = qp[i];
else
*dp = qp[j];
dp++;
}
return mat;
} else if (cdata.is2D() && (cdata.rows() == zrows) &&
(cdata.columns() == zcols)) {
return cdata;
}
}
// In auto mode, or the given data is bogus...
Array mat(Array::doubleMatrixConstructor(zrows,zcols));
double *dp = (double*) mat.getReadWriteDataPointer();
for (int i=0;i<zcols;i++)
for (int j=0;j<zrows;j++) {
if (isXcoord)
*dp = i+1;
else
*dp = j+1;
dp++;
}
return mat;
}
void HandleContour::UpdateState() {
if (HasChanged("levellist")) ToManual("levellistmode");
if (HasChanged("xdata")) ToManual("xdatamode");
if (HasChanged("ydata")) ToManual("ydatamode");
Array zdata(ArrayPropertyLookup("zdata"));
double zmin = ArrayMin(zdata);
double zmax = ArrayMax(zdata);
Array xdata(GetCoordinateMatrix("xdata",true));
Array ydata(GetCoordinateMatrix("ydata",false));
QList<double> levels;
if (StringCheck("levellistmode","auto")) {
levels = GetTicksInner(zmin,zmax,false,10);
if (levels.front() == zmin) levels.pop_front();
if (levels.back() == zmax) levels.pop_back();
std::vector<double> ulevels;
for (int i=0;i<levels.size();i++)
ulevels.push_back(levels[i]);
HPVector *hp = (HPVector*) LookupProperty("levellist");
hp->Data(ulevels);
hp->ClearModified();
} else {
std::vector<double> ulevels(VectorPropertyLookup("levellist"));
for (int i=0;i<ulevels.size();i++)
levels.push_back(ulevels[i]);
}
pset.clear();
zvals.clear();
for (int i=0;i<levels.size();i++) {
pset << ContourCDriver(zdata,levels[i],xdata,ydata);
zvals << levels[i];
}
RebuildContourMatrix();
}
void HandleContour::SelectColor(RenderEngine& gc, double zval) {
// Need to select a color for the contour
// Retrieve the colormap
std::vector<double> cmap(((HandleObject*)GetParentFigure())->VectorPropertyLookup("colormap"));
HandleAxis* ap(GetParentAxis());
std::vector<double> clim(((HandleObject*)ap)->VectorPropertyLookup("clim"));
double clim_min(qMin(clim[0],clim[1]));
double clim_max(qMax(clim[0],clim[1]));
// Calculate the colormap length
int cmaplen(cmap.size()/3);
int ndx = (int) ((zval-clim_min)/(clim_max-clim_min)*(cmaplen-1));
ndx = qMin(cmaplen-1,qMax(0,ndx));
std::vector<double> color;
color.push_back(cmap[3*ndx]);
color.push_back(cmap[3*ndx+1]);
color.push_back(cmap[3*ndx+2]);
gc.color(color);
}
void HandleContour::PaintMe(RenderEngine& gc) {
if (StringCheck("visible","off"))
return;
// Draw the line...
double width(ScalarPropertyLookup("linewidth"));
HPAutoColor *lc = (HPAutoColor*) LookupProperty("linecolor");
bool floatflag = StringCheck("floating","on");
if (!StringCheck("linecolor","none")) {
gc.setLineStyle(StringPropertyLookup("linestyle"));
gc.lineWidth(width);
HandleAxis *parent = (HandleAxis*) GetParentAxis();
for (int i=0;i<pset.size();i++) {
// For each level.
lineset cset(pset[i]);
if (StringCheck("linecolor","auto"))
SelectColor(gc,zvals[i]);
else
gc.color(lc->ColorSpec());
for (int j=0;j<cset.size();j++) {
cline aline(cset[j]);
std::vector<double> xs;
std::vector<double> ys;
std::vector<double> zs;
for (int k=0;k<aline.size();k++) {
xs.push_back(aline[k].x);
ys.push_back(aline[k].y);
if (floatflag)
zs.push_back(zvals[i]);
else
zs.push_back(0);
}
std::vector<double> mxs, mys, mzs;
parent->ReMap(xs,ys,zs,mxs,mys,mzs);
gc.lineSeries(mxs,mys,mzs);
}
}
}
}
void HandleContour::SetupDefaults() {
SetConstrainedStringDefault("linecolor","auto");
SetConstrainedStringDefault("levellistmode","auto");
SetConstrainedStringDefault("linestyle","-");
SetScalarDefault("linewidth",1.0);
SetConstrainedStringDefault("floating","off");
SetStringDefault("type","contour");
SetConstrainedStringDefault("visible","on");
SetConstrainedStringDefault("xdatamode","auto");
SetConstrainedStringDefault("ydatamode","auto");
}
void HandleContour::ConstructProperties() {
//!
//@Module COUNTOUR Contour Object Properties
//@@Section HANDLE
//@@Usage
//Below is a summary of the properties for a line series.
//\begin{itemize}
// \item @|contourmatrix| - @|array| - the matrix containing contour data
// for the plot. This is a @|2 x N| matrix containing x and y coordinates
// for points on the contours. In addition, each contour line starts with
// a column containing the number of points and the contour value.
// \item @|displayname| - @|string| - The name of this line series as it
// appears in a legend.
// \item @|floating| - @|{'on','off'}| - set to on to have floating (3D) contours
// \item @|levellist| - @|vector| - a vector of Z-values for the contours
// \item @|levellistmode| - @|{'auto','manual'}| - set to auto for
// automatic selection of Z-values of the contours.
// \item @|linecolor| - color of the contour lines.
// \item @|linestyle| - @|{'-','--',':','-.','none'}| - the line style to draw the contour in.
// \item @|linewidth| - @|scalar| - the width of the lines
// \item @|parent| - @|handle| - The axis that contains this object
// \item @|tag| - @|string| - A string that can be used to tag the object.
// \item @|type| - @|string| - Returns the string @|'contour'|.
// \item @|userdata| - @|array| - Available to store any variable you
// want in the handle object.
// \item @|visible| - @|{'on','off'}| - Controls visibility of the the line.
// \item @|xdata| - @|matrix| - Contains the x coordinates of the
// surface on which the zdata is defined. This can either be a monotonic
// vector of the same number of columns as @|zdata|, or a 2D matrix
// that is the same size as @|zdata|.
// \item @|xdatamode| - @|{'auto','manual'}| - When set to @|'auto'|
//FreeMat will autogenerate the x coordinates for the contours.
//These values will be @|1,..,N| where @|N| is the number of columns
//of @|zdata|.
// \item @|ydata| - @|matrix| - Contains the y coordinates of the
// surface on which the zdata is defined. This can either be a monotonic
// vector of the same number of rows as @|zdata| or a 2D matrix that is
// the same size as @|zdata|.
// \item @|ydatamode| - @|{'auto','manual'}| - When set to @|'auto'|
//FreeMat will autogenerate the y coordinates for the contour data.
// \item @|zdata| - @|matrix| - The matrix of z values that are to
// be contoured.
//\end{itemize}
//!
AddProperty(new HPArray,"contourmatrix"); //done
AddProperty(new HPHandles,"children"); //done
AddProperty(new HPString,"displayname"); //done
AddProperty(new HPOnOff,"floating"); //done
AddProperty(new HPVector,"levellist"); //done
AddProperty(new HPAutoManual,"levellistmode"); //done
AddProperty(new HPAutoColor,"linecolor"); //done
AddProperty(new HPLineStyle,"linestyle"); //done
AddProperty(new HPScalar,"linewidth"); //done
AddProperty(new HPHandles,"parent"); //done
AddProperty(new HPString,"tag"); //done
AddProperty(new HPString,"type"); //done
AddProperty(new HPArray,"userdata"); //done
AddProperty(new HPOnOff,"visible"); //done
AddProperty(new HPArray,"xdata"); //done
AddProperty(new HPAutoManual,"xdatamode"); //done
AddProperty(new HPArray,"ydata"); //done
AddProperty(new HPAutoManual,"ydatamode"); //done
AddProperty(new HPArray,"zdata"); //done
}
syntax highlighted by Code2HTML, v. 0.9.1