2 * Copyright (c) 2005 Casper Boemann <cbr@boemann.dk>
3 * Copyright (c) 2009 Dmitry Kazakov <dimula73@gmail.com>
4 * Copyright (c) 2010 Cyrille Berger <cberger@cberger.net>
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21 #include "kis_cubic_curve.h"
25 #include <QSharedData>
26 #include <QStringList>
30 class KisTridiagonalSystem
34 * |b0 c0 0 0 0| |x0| |f0|
35 * |a0 b1 c1 0 0| |x1| |f1|
36 * |0 a1 b2 c2 0|*|x2|=|f2|
37 * |0 0 a2 b3 c3| |x3| |f3|
38 * |0 0 0 a3 b4| |x4| |f4|
44 * @return - vector that is storing x[]
47 QVector<T> calculate(QList<T> &a,
58 Q_ASSERT(a.size() == size - 1 &&
59 c.size() == size - 1 &&
65 * Check for special case when
66 * order of the matrix is equal to 1
81 alpha[1] = -c[0] / b[0];
82 beta[1] = f[0] / b[0];
84 for (i = 1; i < size - 1; i++) {
86 (a[i-1] * alpha[i] + b[i]);
88 beta[i+1] = (f[i] - a[i-1] * beta[i])
90 (a[i-1] * alpha[i] + b[i]);
93 x.last() = (f.last() - a.last() * beta.last())
95 (b.last() + a.last() * alpha.last());
97 for (i = size - 2; i >= 0; i--)
98 x[i] = alpha[i+1] * x[i+1] + beta[i+1];
104 template <typename T_point, typename T>
110 * 1/2 * c[i] * (x-x[i])^2 +
111 * 1/6 * d[i] * (x-x[i])^3
130 KisCubicSpline(const QList<T_point> &a) {
135 * Create new spline and precalculate some values
138 * @a - base points of the spline
140 void createSpline(const QList<T_point> &a) {
141 int intervals = m_intervals = a.size() - 1;
143 m_begin = a.first().x();
144 m_end = a.last().x();
147 m_b.resize(intervals);
149 m_d.resize(intervals);
150 m_h.resize(intervals);
152 for (i = 0; i < intervals; i++) {
153 m_h[i] = a[i+1].x() - a[i].x();
154 m_a.append(a[i].y());
156 m_a.append(a.last().y());
161 QList<T> tri_a; /* equals to @tri_c */
163 for (i = 0; i < intervals - 1; i++) {
164 tri_b.append(2.*(m_h[i] + m_h[i+1]));
166 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]));
168 for (i = 1; i < intervals - 1; i++)
169 tri_a.append(m_h[i]);
172 KisTridiagonalSystem<T> tridia;
173 m_c = tridia.calculate(tri_a, tri_b, tri_a, tri_f);
178 for (i = 0; i < intervals; i++)
179 m_d[i] = (m_c[i+1] - m_c[i]) / m_h[i];
181 for (i = 0; i < intervals; i++)
182 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];
186 * Get value of precalculated spline in the point @x
188 T getValue(T x) const {
190 int i = findRegion(x, x0);
191 /* TODO: check for asm equivalent */
194 0.5 * m_c[i] *(x - x0) *(x - x0) +
195 (1 / 6.0)* m_d[i] *(x - x0) *(x - x0) *(x - x0);
209 * findRegion - Searches for the region containing @x
210 * @x0 - out parameter, containing beginning of the region
211 * @return - index of the region
213 int findRegion(T x, T &x0) const {
216 for (i = 0; i < m_intervals; i++) {
217 if (x >= x0 && x < x0 + m_h[i])
222 x0 -= m_h[m_intervals-1];
223 return m_intervals - 1;
226 qDebug("X value: %f\n", x);
227 qDebug("m_begin: %f\n", m_begin);
228 qDebug("m_end : %f\n", m_end);
229 Q_ASSERT_X(0, "findRegion", "X value is outside regions");
230 /* **never reached** */
235 static bool pointLessThan(const QPointF &a, const QPointF &b)
237 return a.x() < b.x();
240 struct KisCubicCurve::Data : public QSharedData {
244 Data(const Data& data) : QSharedData() {
246 points = data.points;
250 validU16Transfer = false;
251 validFTransfer = false;
255 mutable KisCubicSpline<QPointF, qreal> spline;
256 QList<QPointF> points;
257 mutable bool validSpline;
258 mutable QVector<quint16> u16Transfer;
259 mutable bool validU16Transfer;
260 mutable QVector<qreal> fTransfer;
261 mutable bool validFTransfer;
264 qreal value(qreal x);
266 template<typename _T_, typename _T2_>
267 void updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size);
270 void KisCubicCurve::Data::updateSpline()
272 if (validSpline) return;
274 spline.createSpline(points);
277 void KisCubicCurve::Data::invalidate()
280 validFTransfer = false;
281 validU16Transfer = false;
284 void KisCubicCurve::Data::keepSorted()
286 qSort(points.begin(), points.end(), pointLessThan);
289 qreal KisCubicCurve::Data::value(qreal x)
292 /* Automatically extend non-existing parts of the curve
293 * (e.g. before the first point) and cut off big y-values
295 x = qBound(spline.begin(), x, spline.end());
296 qreal y = spline.getValue(x);
297 return qBound((qreal)0.0, y, (qreal)1.0);
300 template<typename _T_, typename _T2_>
301 void KisCubicCurve::Data::updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size)
303 if (!valid || transfer->size() != size) {
304 if (transfer->size() != size) {
305 transfer->resize(size);
307 qreal end = 1.0 / (size - 1);
308 for (int i = 0; i < size; ++i) {
309 /* Direct uncached version */
310 _T2_ val = value(i * end) * max;
311 val = qBound(min, val, max);
312 (*transfer)[i] = val;
318 struct KisCubicCurve::Private {
319 QSharedDataPointer<Data> data;
322 KisCubicCurve::KisCubicCurve() : d(new Private)
328 d->data->points.append(p);
331 d->data->points.append(p);
334 KisCubicCurve::KisCubicCurve(const QList<QPointF>& points) : d(new Private)
337 d->data->points = points;
338 d->data->keepSorted();
341 KisCubicCurve::KisCubicCurve(const KisCubicCurve& curve) : d(new Private(*curve.d))
345 KisCubicCurve::~KisCubicCurve()
350 KisCubicCurve& KisCubicCurve::operator=(const KisCubicCurve & curve)
356 bool KisCubicCurve::operator==(const KisCubicCurve& curve) const
358 if (d->data == curve.d->data) return true;
359 return d->data->points == curve.d->data->points;
362 qreal KisCubicCurve::value(qreal x) const
364 return d->data->value(x);
367 QList<QPointF> KisCubicCurve::points() const
369 return d->data->points;
372 void KisCubicCurve::setPoints(const QList<QPointF>& points)
375 d->data->points = points;
376 d->data->invalidate();
379 void KisCubicCurve::setPoint(int idx, const QPointF& point)
382 d->data->points[idx] = point;
383 d->data->keepSorted();
384 d->data->invalidate();
387 int KisCubicCurve::addPoint(const QPointF& point)
390 d->data->points.append(point);
391 d->data->keepSorted();
392 d->data->invalidate();
393 return d->data->points.indexOf(point);
396 void KisCubicCurve::removePoint(int idx)
399 d->data->points.removeAt(idx);
400 d->data->invalidate();
403 QString KisCubicCurve::toString() const
407 foreach(const QPointF & pair, d->data->points) {
408 sCurve += locale.toString(pair.x());
410 sCurve += QString::number(pair.y());
416 void KisCubicCurve::fromString(const QString& string)
418 QStringList data = string.split(';');
420 QList<QPointF> points;
422 foreach(const QString & pair, data) {
423 if (pair.indexOf('/') > -1) {
425 p.rx() = pair.section('/', 0, 0).toDouble();
426 p.ry() = pair.section('/', 1, 1).toDouble();
433 QVector<quint16> KisCubicCurve::uint16Transfer(int size) const
435 d->data->updateTransfer<quint16, int>(&d->data->u16Transfer, d->data->validU16Transfer, 0x0, 0xFFFF, size);
436 return d->data->u16Transfer;
439 QVector<qreal> KisCubicCurve::floatTransfer(int size) const
441 d->data->updateTransfer<qreal, qreal>(&d->data->fTransfer, d->data->validFTransfer, 0.0, 1.0, size);
442 return d->data->fTransfer;