mirrors/JKQtPlotter
1
0
Fork 0
mirror of https://github.com/jkriege2/JKQtPlotter.git synced 2024-11-16 02:25:50 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity
028f192e7b
JKQtPlotter/lib/jkqtplotter/jkqtpgraphscontour.cpp

520 lines
17 KiB
C++
Raw Blame History

/*
Copyright (c) 2008-2019 Jan W. Krieger & Sebastian Isbaner (contour plot)
This software is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License (LGPL) as published by
the Free Software Foundation, either version 2.1 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 Lesser General Public License (LGPL) for more details.
You should have received a copy of the GNU Lesser General Public License (LGPL)
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "jkqtplotter/jkqtpgraphscontour.h"
#include "jkqtplotter/jkqtpbaseplotter.h"
#include "jkqtplottertools/jkqtpimagetools.h"
#include "jkqtcommon/jkqtptools.h"
#include "jkqtplottertools/jkqtpenhancedpainter.h"
#include "jkqtplotter/jkqtplotter.h"
#include <QDebug>
#include <QImageWriter>
#include <QFileDialog>
#include <QFileInfo>
#include <QApplication>
#include <QClipboard>
# include <QVector3D>
JKQTPContour::JKQTPContour(JKQTBasePlotter *parent) :
JKQTPMathImage(parent)
{
colorBarRightVisible=false;
ignoreOnPlane=false;
numberOfLevels=1;
colorFromPalette=true;
datatype=JKQTPMathImageBase::DoubleArray;
relativeLevels=false;
initLineStyle(parent, parentPlotStyle);
}
JKQTPContour::JKQTPContour(double x, double y, double width, double height, void* data, int Nx, int Ny, JKQTPMathImageColorPalette palette, DataType datatype, JKQTBasePlotter* parent) :
JKQTPMathImage( x, y, width, height, datatype, data, Nx, Ny, palette, parent)
{
colorBarRightVisible=false;
ignoreOnPlane=false;
numberOfLevels=1;
colorFromPalette=true;
relativeLevels=false;
initLineStyle(parent, parentPlotStyle);
}
JKQTPContour::JKQTPContour(JKQTPlotter *parent) :
JKQTPContour(parent->getPlotter())
{
}
JKQTPContour::JKQTPContour(double x, double y, double width, double height, void* data, int Nx, int Ny, JKQTPMathImageColorPalette palette, DataType datatype, JKQTPlotter* parent) :
JKQTPContour( x, y, width, height, data, Nx, Ny, palette, datatype, parent->getPlotter())
{
}
void JKQTPContour::draw(JKQTPEnhancedPainter &painter)
{
//qDebug()<<"JKQTPContourPlot::draw";
if(contourLevels.isEmpty()) createContourLevels(numberOfLevels);
else {
numberOfLevels=contourLevels.size();
// qSort(contourLevels);
}
if(contourLines.isEmpty()) { // contour lines are only calculated once
for(int i =0; i<contourLevels.size();++i) {
contourLines.append(QVector<QLineF> (0));
}
this->calcContourLines(contourLines);
}
// draw lines
painter.save(); auto __finalpaint=JKQTPFinally([&painter]() {painter.restore();});
QPen p=getLinePen(painter, parent);
painter.setPen(p);
QImage colorLevels = getPaletteImage(palette,numberOfLevels);
QVector<QLineF> contourLinesTransformedSingleLevel;
QLineF lineTranformed;
for(int i =0; i<numberOfLevels;++i) {
if(colorFromPalette) {
p.setColor(QColor(colorLevels.pixel(i,0)));
painter.setPen(p);
}
// transform into plot coordinates
for(QVector<QLineF >::const_iterator line =contourLines.at(i).begin(); line!=contourLines.at(i).end();++line ) {
lineTranformed.setP1(transform(x+line->p1().x()/double(Nx-1)*width, y+line->p1().y()/double(Ny-1)*height));
lineTranformed.setP2(transform(x+line->p2().x()/double(Nx-1)*width, y+line->p2().y()/double(Ny-1)*height));
contourLinesTransformedSingleLevel.append(lineTranformed);
}
painter.drawLines(contourLinesTransformedSingleLevel);
contourLinesTransformedSingleLevel.clear();
}
}
void JKQTPContour::createContourLevels(int nLevels)
{
if (!data) return;
if (nLevels<1) return;
double min,max;
getDataMinMax(min,max);
double delta=(max-min)/static_cast<double>(nLevels+1);
for(int i=1; i<=nLevels; ++i) {
contourLevels.append(min + i*delta);
}
relativeLevels=false;
}
void JKQTPContour::createContourLevelsLog(int nLevels, int m)
{
if (!data) return;
if (nLevels<1) return;
double min,max;
getDataMinMax(min,max);
if(min<=0) min=1; // FIXME get smallest number greater zero
int S=floor((log10(max)-log10(min))/log10(m));
if(S<2) S=1;
int P = floor(static_cast<double>(nLevels)/static_cast<double>(S));
if(P<1) P=1;
double delta=min;
contourLevels.append(2*delta);
for (long s=0; s<S; s++) {
for (long p=0; p<P; p++) {
{
contourLevels.append(contourLevels.last()+delta);
}
}
delta=delta*m;
}
if(nLevels!=contourLevels.size()) {
//qDebug()<<"nLevels="<<nLevels<<"contourLevels.size()="<<contourLevels.size();
//qDebug()<<"adapt m";
}
relativeLevels=false;
}
void JKQTPContour::setIgnoreOnPlane(bool __value)
{
this->ignoreOnPlane = __value;
}
bool JKQTPContour::getIgnoreOnPlane() const
{
return this->ignoreOnPlane;
}
void JKQTPContour::setNumberOfLevels(int __value)
{
this->numberOfLevels = __value;
}
int JKQTPContour::getNumberOfLevels() const
{
return this->numberOfLevels;
}
void JKQTPContour::setColorFromPalette(bool __value)
{
this->colorFromPalette = __value;
}
bool JKQTPContour::getColorFromPalette() const
{
return this->colorFromPalette;
}
void JKQTPContour::setContourLevels(const QList<double> &__value)
{
this->contourLevels = __value;
}
QList<double> JKQTPContour::getContourLevels() const
{
return this->contourLevels;
}
void JKQTPContour::setRelativeLevels(bool __value)
{
this->relativeLevels = __value;
}
bool JKQTPContour::getRelativeLevels() const
{
return this->relativeLevels;
}
void JKQTPContour::setImageColumn(size_t columnID)
{
datatype=JKQTPMathImageBase::DoubleArray;
data=parent->getDatastore()->getColumnPointer(columnID,0);
if (parent && columnID>=0 && parent->getDatastore()) {
setNx(parent->getDatastore()->getColumnImageWidth(columnID));
setNy(parent->getDatastore()->getColumnImageHeight(columnID));
}
}
void JKQTPContour::ensureImageData()
{
}
double JKQTPContour::value(int xIdx, int yIdx)
{
// row-major in datastore
ensureImageData();
if (!data) return 0;
switch(datatype) {
case JKQTPMathImageBase::DoubleArray:
return (static_cast<double*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::FloatArray:
return (static_cast<float*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::UInt8Array:
return (static_cast<uint8_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::UInt16Array:
return (static_cast<uint16_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::UInt32Array:
return (static_cast<uint32_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::UInt64Array:
return (static_cast<uint64_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::Int8Array:
return (static_cast<int8_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::Int16Array:
return (static_cast<int16_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::Int32Array:
return (static_cast<int32_t*>(data))[yIdx*getNx()+xIdx];
case JKQTPMathImageBase::Int64Array:
return (static_cast<int64_t*>(data))[yIdx*getNx()+xIdx];
default:
return 0;
}
}
bool JKQTPContour::intersect(QLineF &line, const QVector3D &vertex1,const QVector3D &vertex2,const QVector3D &vertex3,double level)
{
bool found = true;
// Are the vertices below (-1), on (0) or above (1) the plane ?
const int eq1 = compare2level(vertex1,level);
const int eq2 = compare2level(vertex2,level);
const int eq3 = compare2level(vertex3,level);
/*
(a) All the vertices lie below the contour level.
(b) Two vertices lie below and one on the contour level.
(c) Two vertices lie below and one above the contour level.
(d) One vertex lies below and two on the contour level.
(e) One vertex lies below, one on and one above the contour level.
(f) One vertex lies below and two above the contour level.
(g) Three vertices lie on the contour level.
(h) Two vertices lie on and one above the contour level.
(i) One vertex lies on and two above the contour level.
(j) All the vertices lie above the contour level.
*/
static const int caseLUT[3][3][3] =
{
// jump table to avoid nested case statements
{ { 0, 0, 8 }, { 0, 2, 5 }, { 7, 6, 9 } },
{ { 0, 3, 4 }, { 1, 10, 1 }, { 4, 3, 0 } },
{ { 9, 6, 7 }, { 5, 2, 0 }, { 8, 0, 0 } }
};
const int caseType = caseLUT[eq1+1][eq2+1][eq3+1];
switch (caseType)
{
case 1:
// d(0,0,-1), h(0,0,1)
line.setP1(vertex1.toPointF());
line.setP2(vertex2.toPointF());
break;
case 2:
// d(-1,0,0), h(1,0,0)
line.setP1(vertex2.toPointF());
line.setP2(vertex3.toPointF());
break;
case 3:
// d(0,-1,0), h(0,1,0)
line.setP1(vertex3.toPointF());
line.setP2(vertex1.toPointF());
break;
case 4:
// e(0,-1,1), e(0,1,-1)
line.setP1(vertex1.toPointF());
line.setP2(interpolatePoint(vertex2, vertex3, level));
break;
case 5:
// e(-1,0,1), e(1,0,-1)
line.setP1(vertex2.toPointF());
line.setP2(interpolatePoint(vertex3, vertex1, level));
break;
case 6:
// e(-1,1,0), e(1,0,-1)
line.setP1(vertex3.toPointF());
line.setP2(interpolatePoint(vertex1, vertex2, level));
break;
case 7:
// c(-1,1,-1), f(1,1,-1)
line.setP1(interpolatePoint(vertex1, vertex2, level));
line.setP2(interpolatePoint(vertex2, vertex3, level));
break;
case 8:
// c(-1,-1,1), f(1,1,-1)
line.setP1(interpolatePoint(vertex2, vertex3, level));
line.setP2(interpolatePoint(vertex3, vertex1, level));
break;
case 9:
// f(-1,1,1), c(1,-1,-1)
line.setP1(interpolatePoint(vertex3, vertex1, level));
line.setP2(interpolatePoint(vertex1, vertex2, level));
break;
case 10:
// g(0,0,0)
// The CONREC algorithm has no satisfying solution for
// what to do, when all vertices are on the plane.
if ( ignoreOnPlane )
found = false;
else
{
line.setP1(vertex3.toPointF());
line.setP2(vertex1.toPointF());
}
break;
default:
found = false;
}
// qDebug()<<caseType;
// qDebug()<<line;
return found;
}
int JKQTPContour::compare2level(const QVector3D &vertex, double level)
{
if (vertex.z() > level)
return 1;
if (vertex.z() < level)
return -1;
return 0;
}
void JKQTPContour::calcContourLines(QList<QVector<QLineF> > &ContourLines)
{
double scale=1; ///< scale of the contour levels;
if(relativeLevels) {
double min;
double max;
getDataMinMax(min,max);
scale=1/(max-min);
}
enum Position
{
// the positions of points of one box
// vertex 1 +-------------------+ vertex 2
// | \ / |
// | \ m=3 / |
// | \ / |
// | \ / |
// | m=2 X m=2 | the center is vertex 0
// | / \ |
// | / \ |
// | / m=1 \ |
// | / \ |
// vertex 4 +-------------------+ vertex 3
Center=0,
TopLeft=1,
TopRight=2,
BottomRight=3,
BottomLeft=4,
NumPositions=5
};
for ( int yp = 0; yp < (int64_t)getNy() - 1; ++yp ) { // go through image (pixel coordinates) in row major order
QVector<QVector3D> vertices(NumPositions);
for ( int xp = 0; xp < (int64_t)getNy() - 1; ++xp ) {
if ( xp == 0 )
{
vertices[TopRight].setX(xp); // will be used for TopLeft later
vertices[TopRight].setY(yp);
vertices[TopRight].setZ(
value( vertices[TopRight].x(), vertices[TopRight].y())*scale
);
vertices[BottomRight].setX(xp);
vertices[BottomRight].setY(yp+1);
vertices[BottomRight].setZ(
value(vertices[BottomRight].x(), vertices[BottomRight].y())*scale
);
}
vertices[TopLeft] = vertices[TopRight]; // use right vertices of the last box as new left vertices
vertices[BottomLeft] = vertices[BottomRight];
vertices[TopRight].setX(xp + 1);
vertices[TopRight].setY(yp); // <----
vertices[TopRight].setZ(
value(vertices[TopRight].x(), vertices[TopRight].y())*scale
);
vertices[BottomRight].setX(xp + 1);
vertices[BottomRight].setY(yp + 1);
vertices[BottomRight].setZ(
value(vertices[BottomRight].x(), vertices[BottomRight].y())*scale
);
double zMin = vertices[TopLeft].z();
double zMax = zMin;
double zSum = zMin;
for ( int i = TopRight; i <= BottomLeft; ++i ) {
const double z = vertices[i].z();
zSum += z;
if ( z < zMin )
zMin = z;
if ( z > zMax )
zMax = z;
}
if ( zMax >= contourLevels.first() && zMin <= contourLevels.last() ) {
vertices[Center].setX(xp + 0.5); // pseudo pixel coordinates
vertices[Center].setY(yp + 0.5);
vertices[Center].setZ(0.25 * zSum);
for (int levelIdx=0; levelIdx<contourLevels.size(); ++levelIdx) {
if ( contourLevels.at(levelIdx) >= zMin && contourLevels.at(levelIdx) <= zMax ) {
QLineF line;
QVector<QVector3D> triangle(3);
/* triangle[1]
X
/ \
/ \
/ m \
/ \
triangle[2] +-------------------+ triangle[0]
*/
for (int m = TopLeft; m < NumPositions; m++) { // construct triangles
triangle[0] = vertices[m];
triangle[1] = vertices[Center];
triangle[2] = vertices[(m!=BottomLeft)?(m + 1):TopLeft];
const bool intersects =intersect(line, triangle.at(0),triangle.at(1),triangle.at(2),
contourLevels.at(levelIdx));
if ( intersects ) {
ContourLines[levelIdx]<<line;
}
}
}
}
}
}
}
}
QPointF JKQTPContour::interpolatePoint(const QVector3D &point1, const QVector3D &point2,double level)
{
const double h1 = point1.z() - level; // height above contour level
const double h2 = point2.z() - level;
// // check if h1 or h2 is zero
// Division by zero is not possible (the intersect function is not called if h2-h1 is zero, !)
// if(h2-h1==0||h1==0||h2==0) {
// qDebug()<<h1<<h2;
// qDebug()<<"interpolate p1="<<point1<<", p2="<<point2<<" level="<<level;
// }
const double x = (h2 * point1.x() - h1 * point2.x()) / (h2 - h1); // linear interpolation in x.direction (independent of y)
const double y = (h2 * point1.y() - h1 * point2.y()) / (h2 - h1);
// const double alpha=(level-point1.z())/(point2.z()-point1.z());
// const double x=point1.x()+alpha*(point2.x()-point1.x());
// const double y=point1.y()+alpha*(point2.y()-point1.y());
return QPointF(x, y);
}
Reference in New Issue View Git Blame Copy Permalink