File indexing completed on 2024-05-19 04:26:10

0001 /*
0002  *  SPDX-FileCopyrightText: 2005 C. Boemann <cbo@boemann.dk>
0003  *  SPDX-FileCopyrightText: 2009 Dmitry Kazakov <dimula73@gmail.com>
0004  *  SPDX-FileCopyrightText: 2010 Cyrille Berger <cberger@cberger.net>
0005  *
0006  *  SPDX-License-Identifier: GPL-2.0-or-later
0007  */
0008 
0009 #include "kis_cubic_curve.h"
0010 
0011 #include <QPointF>
0012 #include <QList>
0013 #include <QSharedData>
0014 #include <QStringList>
0015 #include "kis_dom_utils.h"
0016 #include "kis_algebra_2d.h"
0017 
0018 template <typename T>
0019 class KisTridiagonalSystem
0020 {
0021     /*
0022      * e.g.
0023      *      |b0 c0  0   0   0| |x0| |f0|
0024      *      |a0 b1 c1   0   0| |x1| |f1|
0025      *      |0  a1 b2  c2   0|*|x2|=|f2|
0026      *      |0   0 a2  b3  c3| |x3| |f3|
0027      *      |0   0  0  a3  b4| |x4| |f4|
0028      */
0029 
0030 public:
0031 
0032     /**
0033      * @return - vector that is storing x[]
0034      */
0035     static
0036     QVector<T> calculate(QList<T> &a,
0037                          QList<T> &b,
0038                          QList<T> &c,
0039                          QList<T> &f) {
0040         QVector<T> x;
0041         QVector<T> alpha;
0042         QVector<T> beta;
0043 
0044         int i;
0045         int size = b.size();
0046 
0047         Q_ASSERT(a.size() == size - 1 &&
0048                  c.size() == size - 1 &&
0049                  f.size() == size);
0050 
0051         x.resize(size);
0052 
0053         /**
0054          * Check for special case when
0055          * order of the matrix is equal to 1
0056          */
0057         if (size == 1) {
0058             x[0] = f[0] / b[0];
0059             return x;
0060         }
0061 
0062         /**
0063          * Common case
0064          */
0065 
0066         alpha.resize(size);
0067         beta.resize(size);
0068 
0069 
0070         alpha[1] = -c[0] / b[0];
0071         beta[1] =  f[0] / b[0];
0072 
0073         for (i = 1; i < size - 1; i++) {
0074             alpha[i+1] = -c[i] /
0075                          (a[i-1] * alpha[i] + b[i]);
0076 
0077             beta[i+1] = (f[i] - a[i-1] * beta[i])
0078                         /
0079                         (a[i-1] * alpha[i] + b[i]);
0080         }
0081 
0082         x.last() = (f.last() - a.last() * beta.last())
0083                    /
0084                    (b.last() + a.last() * alpha.last());
0085 
0086         for (i = size - 2; i >= 0; i--)
0087             x[i] = alpha[i+1] * x[i+1] + beta[i+1];
0088 
0089         return x;
0090     }
0091 };
0092 
0093 template <typename T_point, typename T>
0094 class KisCubicSpline
0095 {
0096     /**
0097      *  s[i](x)=a[i] +
0098      *          b[i] * (x-x[i]) +
0099      *    1/2 * c[i] * (x-x[i])^2 +
0100      *    1/6 * d[i] * (x-x[i])^3
0101      *
0102      *  h[i]=x[i+1]-x[i]
0103      *
0104      */
0105 
0106 protected:
0107     QList<T> m_a;
0108     QVector<T> m_b;
0109     QVector<T> m_c;
0110     QVector<T> m_d;
0111 
0112     QVector<T> m_h;
0113     T m_begin;
0114     T m_end;
0115     int m_intervals {0};
0116 
0117 public:
0118     KisCubicSpline() {}
0119     KisCubicSpline(const QList<T_point> &a) {
0120         createSpline(a);
0121     }
0122 
0123     /**
0124      * Create new spline and precalculate some values
0125      * for future
0126      *
0127      * @a - base points of the spline
0128      */
0129     void createSpline(const QList<T_point> &a) {
0130         int intervals = m_intervals = a.size() - 1;
0131         int i;
0132         m_begin = a.first().x();
0133         m_end = a.last().x();
0134 
0135         m_a.clear();
0136         m_b.resize(intervals);
0137         m_c.clear();
0138         m_d.resize(intervals);
0139         m_h.resize(intervals);
0140 
0141         for (i = 0; i < intervals; i++) {
0142             m_h[i] = a[i+1].x() - a[i].x();
0143             m_a.append(a[i].y());
0144         }
0145         m_a.append(a.last().y());
0146 
0147 
0148         QList<T> tri_b;
0149         QList<T> tri_f;
0150         QList<T> tri_a; /* equals to @tri_c */
0151 
0152         for (i = 0; i < intervals - 1; i++) {
0153             tri_b.append(2.*(m_h[i] + m_h[i+1]));
0154 
0155             tri_f.append(6.*((m_a[i+2] - m_a[i+1]) / m_h[i+1] - (m_a[i+1] - m_a[i]) / m_h[i]));
0156         }
0157         for (i = 1; i < intervals - 1; i++)
0158             tri_a.append(m_h[i]);
0159 
0160         if (intervals > 1) {
0161             m_c = KisTridiagonalSystem<T>::calculate(tri_a, tri_b, tri_a, tri_f);
0162         }
0163         m_c.prepend(0);
0164         m_c.append(0);
0165 
0166         for (i = 0; i < intervals; i++)
0167             m_d[i] = (m_c[i+1] - m_c[i]) / m_h[i];
0168 
0169         for (i = 0; i < intervals; i++)
0170             m_b[i] = -0.5 * (m_c[i] * m_h[i])  - (1 / 6.0) * (m_d[i] * m_h[i] * m_h[i]) + (m_a[i+1] - m_a[i]) / m_h[i];
0171     }
0172 
0173     /**
0174      * Get value of precalculated spline in the point @x
0175      */
0176     T getValue(T x) const {
0177         T x0;
0178         int i = findRegion(x, x0);
0179         /* TODO: check for asm equivalent */
0180         return m_a[i] +
0181                m_b[i] *(x - x0) +
0182                0.5 * m_c[i] *(x - x0) *(x - x0) +
0183                (1 / 6.0)* m_d[i] *(x - x0) *(x - x0) *(x - x0);
0184     }
0185 
0186     T begin() const {
0187         return m_begin;
0188     }
0189 
0190     T end() const {
0191         return m_end;
0192     }
0193 
0194 protected:
0195 
0196     /**
0197      * findRegion - Searches for the region containing @x
0198      * @x0 - out parameter, containing beginning of the region
0199      * @return - index of the region
0200      */
0201     int findRegion(T x, T &x0) const {
0202         int i;
0203         x0 = m_begin;
0204         for (i = 0; i < m_intervals; i++) {
0205             if (x >= x0 && x < x0 + m_h[i])
0206                 return i;
0207             x0 += m_h[i];
0208         }
0209         if (x >= x0) {
0210             x0 -= m_h[m_intervals-1];
0211             return m_intervals - 1;
0212         }
0213 
0214         qDebug("X value: %f\n", x);
0215         qDebug("m_begin: %f\n", m_begin);
0216         qDebug("m_end  : %f\n", m_end);
0217         Q_ASSERT_X(0, "findRegion", "X value is outside regions");
0218         /* **never reached** */
0219         return -1;
0220     }
0221 };
0222 
0223 static bool pointLessThan(const QPointF &a, const QPointF &b)
0224 {
0225     return a.x() < b.x();
0226 }
0227 
0228 struct Q_DECL_HIDDEN KisCubicCurve::Data : public QSharedData {
0229     Data() {
0230     }
0231     Data(const Data& data) : QSharedData() {
0232         points = data.points;
0233         name = data.name;
0234     }
0235     ~Data() {
0236     }
0237 
0238     mutable QString name;
0239     mutable KisCubicSpline<QPointF, qreal> spline;
0240     QList<QPointF> points;
0241     mutable bool validSpline {false};
0242     mutable QVector<quint8> u8Transfer;
0243     mutable bool validU8Transfer {false};
0244     mutable QVector<quint16> u16Transfer;
0245     mutable bool validU16Transfer {false};
0246     mutable QVector<qreal> fTransfer;
0247     mutable bool validFTransfer {false};
0248 
0249     void updateSpline();
0250     void keepSorted();
0251     qreal value(qreal x);
0252     void invalidate();
0253 
0254     template<typename _T_, typename _T2_>
0255     void updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size);
0256 };
0257 
0258 void KisCubicCurve::Data::updateSpline()
0259 {
0260     if (validSpline) return;
0261     validSpline = true;
0262     spline.createSpline(points);
0263 }
0264 
0265 void KisCubicCurve::Data::invalidate()
0266 {
0267     validSpline = false;
0268     validFTransfer = false;
0269     validU16Transfer = false;
0270 }
0271 
0272 void KisCubicCurve::Data::keepSorted()
0273 {
0274     std::sort(points.begin(), points.end(), pointLessThan);
0275 }
0276 
0277 qreal KisCubicCurve::Data::value(qreal x)
0278 {
0279     updateSpline();
0280     /* Automatically extend non-existing parts of the curve
0281      * (e.g. before the first point) and cut off big y-values
0282      */
0283     x = qBound(spline.begin(), x, spline.end());
0284     qreal y = spline.getValue(x);
0285     return qBound(qreal(0.0), y, qreal(1.0));
0286 }
0287 
0288 template<typename _T_, typename _T2_>
0289 void KisCubicCurve::Data::updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size)
0290 {
0291     if (!valid || transfer->size() != size) {
0292         if (transfer->size() != size) {
0293             transfer->resize(size);
0294         }
0295         qreal end = 1.0 / (size - 1);
0296         for (int i = 0; i < size; ++i) {
0297             /* Direct uncached version */
0298             _T2_ val = value(i * end ) * max;
0299             val = qBound(min, val, max);
0300             (*transfer)[i] = val;
0301         }
0302         valid = true;
0303     }
0304 }
0305 
0306 struct Q_DECL_HIDDEN KisCubicCurve::Private {
0307     QSharedDataPointer<Data> data;
0308 };
0309 
0310 KisCubicCurve::KisCubicCurve()
0311     : d(new Private)
0312 {
0313     d->data = new Data;
0314     QPointF p;
0315     p.rx() = 0.0; p.ry() = 0.0;
0316     d->data->points.append(p);
0317     p.rx() = 1.0; p.ry() = 1.0;
0318     d->data->points.append(p);
0319 }
0320 
0321 KisCubicCurve::KisCubicCurve(const QList<QPointF>& points) : d(new Private)
0322 {
0323     d->data = new Data;
0324     d->data->points = points;
0325     d->data->keepSorted();
0326 }
0327 
0328 KisCubicCurve::KisCubicCurve(const QVector<QPointF> &points)
0329     : KisCubicCurve(points.toList())
0330 {
0331 }
0332 
0333 KisCubicCurve::KisCubicCurve(const KisCubicCurve& curve)
0334     : d(new Private(*curve.d))
0335 {
0336 }
0337 
0338 KisCubicCurve::KisCubicCurve(const QString &curveString)
0339     : d(new Private)
0340 {
0341     d->data = new Data;
0342 
0343     KIS_SAFE_ASSERT_RECOVER(!curveString.isEmpty()) {
0344         *this = KisCubicCurve();
0345         return;
0346     }
0347 
0348     const QStringList data = curveString.split(';');
0349 
0350     QList<QPointF> points;
0351     Q_FOREACH (const QString & pair, data) {
0352         if (pair.indexOf(',') > -1) {
0353             QPointF p;
0354             p.rx() = KisDomUtils::toDouble(pair.section(',', 0, 0));
0355             p.ry() = KisDomUtils::toDouble(pair.section(',', 1, 1));
0356             points.append(p);
0357         }
0358     }
0359 
0360     d->data->points = points;
0361     setPoints(points);
0362 }
0363 
0364 KisCubicCurve::~KisCubicCurve()
0365 {
0366     delete d;
0367 }
0368 
0369 KisCubicCurve& KisCubicCurve::operator=(const KisCubicCurve & curve)
0370 {
0371     if (&curve != this) {
0372         *d = *curve.d;
0373     }
0374     return *this;
0375 }
0376 
0377 bool KisCubicCurve::operator==(const KisCubicCurve& curve) const
0378 {
0379     if (d->data == curve.d->data) return true;
0380     return d->data->points == curve.d->data->points;
0381 }
0382 
0383 qreal KisCubicCurve::value(qreal x) const
0384 {
0385     qreal value = d->data->value(x);
0386     return value;
0387 }
0388 
0389 const QList<QPointF>& KisCubicCurve::points() const
0390 {
0391     return d->data->points;
0392 }
0393 
0394 void KisCubicCurve::setPoints(const QList<QPointF>& points)
0395 {
0396     d->data.detach();
0397     d->data->points = points;
0398     d->data->invalidate();
0399 }
0400 
0401 void KisCubicCurve::setPoint(int idx, const QPointF& point)
0402 {
0403     d->data.detach();
0404     d->data->points[idx] = point;
0405     d->data->keepSorted();
0406     d->data->invalidate();
0407 }
0408 
0409 int KisCubicCurve::addPoint(const QPointF& point)
0410 { 
0411     d->data.detach();
0412     d->data->points.append(point);
0413     d->data->keepSorted();
0414     d->data->invalidate();
0415 
0416     return d->data->points.indexOf(point);
0417 }
0418 
0419 void KisCubicCurve::removePoint(int idx)
0420 {
0421     d->data.detach();
0422     d->data->points.removeAt(idx);
0423     d->data->invalidate();
0424 }
0425 
0426 bool KisCubicCurve::isIdentity() const
0427 {
0428     const QList<QPointF> &points = d->data->points;
0429     const int size = points.size();
0430 
0431     if (points[0] != QPointF(0,0) || points[size-1] != QPointF(1,1)) {
0432         return false;
0433     }
0434 
0435     for (int i = 1; i < size-1; i++) {
0436         if (!qFuzzyCompare(points[i].x(), points[i].y())) {
0437             return false;
0438         }
0439     }
0440 
0441     return true;
0442 }
0443 
0444 bool KisCubicCurve::isConstant(qreal c) const
0445 {
0446     const QList<QPointF> &points = d->data->points;
0447 
0448     Q_FOREACH (const QPointF &pt, points) {
0449             if (!qFuzzyCompare(c, pt.y())) {
0450                 return false;
0451             }
0452         }
0453 
0454     return true;
0455 }
0456 
0457 const QString& KisCubicCurve::name() const
0458 {
0459     return d->data->name;
0460 }
0461 
0462 qreal KisCubicCurve::interpolateLinear(qreal normalizedValue, const QVector<qreal> &transfer)
0463 {
0464     const qreal maxValue = transfer.size() - 1;
0465 
0466     const qreal bilinearX = qBound(0.0, maxValue * normalizedValue, maxValue);
0467     const qreal xFloored = std::floor(bilinearX);
0468     const qreal xCeiled = std::ceil(bilinearX);
0469 
0470     const qreal t = bilinearX - xFloored;
0471 
0472     constexpr qreal eps = 1e-6;
0473 
0474     qreal newValue = normalizedValue;
0475 
0476     if (t < eps) {
0477         newValue = transfer[int(xFloored)];
0478     } else if (t > (1.0 - eps)) {
0479         newValue = transfer[int(xCeiled)];
0480     } else {
0481         qreal a = transfer[int(xFloored)];
0482         qreal b = transfer[int(xCeiled)];
0483 
0484         newValue = a + t * (b - a);
0485     }
0486 
0487     return KisAlgebra2D::copysign(newValue, normalizedValue);
0488 }
0489 
0490 void KisCubicCurve::setName(const QString& name)
0491 {
0492     d->data->name = name;
0493 }
0494 
0495 QString KisCubicCurve::toString() const
0496 {
0497     QString sCurve;
0498 
0499     if(d->data->points.count() < 1)
0500         return sCurve;
0501 
0502     Q_FOREACH (const QPointF & pair, d->data->points) {
0503         sCurve += QString::number(pair.x());
0504         sCurve += ',';
0505         sCurve += QString::number(pair.y());
0506         sCurve += ';';
0507     }
0508 
0509     return sCurve;
0510 }
0511 
0512 void KisCubicCurve::fromString(const QString& string)
0513 {
0514     *this = KisCubicCurve(string);
0515 }
0516 
0517 const QVector<quint16> KisCubicCurve::uint16Transfer(int size) const
0518 {
0519     d->data->updateTransfer<quint16, int>(&d->data->u16Transfer, d->data->validU16Transfer, 0x0, 0xFFFF, size);
0520     return d->data->u16Transfer;
0521 }
0522 
0523 const QVector<qreal> KisCubicCurve::floatTransfer(int size) const
0524 {
0525     d->data->updateTransfer<qreal, qreal>(&d->data->fTransfer, d->data->validFTransfer, 0.0, 1.0, size);
0526     return d->data->fTransfer;
0527 }