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 }