git.sesse.net Git - kdenlive/blob - src/kis_cubic_curve.cpp

   1 /*
   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>
   5  *
   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.
  10  *
  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.
  15  *
  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.
  19  */
  20
  21 #include "kis_cubic_curve.h"
  22
  23 #include <QPointF>
  24 #include <QList>
  25 #include <QSharedData>
  26 #include <QStringList>
  27
  28 template <typename T>
  29 class KisTridiagonalSystem
  30 {
  31     /*
  32      * e.g.
  33      *      |b0 c0  0   0   0| |x0| |f0|
  34      *      |a0 b1 c1   0   0| |x1| |f1|
  35      *      |0  a1 b2  c2   0|*|x2|=|f2|
  36      *      |0   0 a2  b3  c3| |x3| |f3|
  37      *      |0   0  0  a3  b4| |x4| |f4|
  38      */
  39
  40 public:
  41
  42     /**
  43      * @return - vector that is storing x[]
  44      */
  45     static
  46     QVector<T> calculate(QList<T> &a,
  47                          QList<T> &b,
  48                          QList<T> &c,
  49                          QList<T> &f) {
  50         QVector<T> x;
  51         QVector<T> alpha;
  52         QVector<T> beta;
  53
  54         int i;
  55         int size = b.size();
  56
  57         Q_ASSERT(a.size() == size - 1 &&
  58                  c.size() == size - 1 &&
  59                  f.size() == size);
  60
  61         x.resize(size);
  62
  63         /**
  64          * Check for special case when
  65          * order of the matrix is equal to 1
  66          */
  67         if (size == 1) {
  68             x[0] = f[0] / b[0];
  69             return x;
  70         }
  71
  72         /**
  73          * Common case
  74          */
  75
  76         alpha.resize(size);
  77         beta.resize(size);
  78
  79
  80         alpha[1] = -c[0] / b[0];
  81         beta[1] =  f[0] / b[0];
  82
  83         for (i = 1; i < size - 1; i++) {
  84             alpha[i+1] = -c[i] /
  85                          (a[i-1] * alpha[i] + b[i]);
  86
  87             beta[i+1] = (f[i] - a[i-1] * beta[i])
  88                         /
  89                         (a[i-1] * alpha[i] + b[i]);
  90         }
  91
  92         x.last() = (f.last() - a.last() * beta.last())
  93                    /
  94                    (b.last() + a.last() * alpha.last());
  95
  96         for (i = size - 2; i >= 0; i--)
  97             x[i] = alpha[i+1] * x[i+1] + beta[i+1];
  98
  99         return x;
 100     }
 101 };
 102
 103 template <typename T_point, typename T>
 104 class KisCubicSpline
 105 {
 106     /**
 107      *  s[i](x)=a[i] +
 108      *          b[i] * (x-x[i]) +
 109      *    1/2 * c[i] * (x-x[i])^2 +
 110      *    1/6 * d[i] * (x-x[i])^3
 111      *
 112      *  h[i]=x[i+1]-x[i]
 113      *
 114      */
 115
 116 protected:
 117     QList<T> m_a;
 118     QVector<T> m_b;
 119     QVector<T> m_c;
 120     QVector<T> m_d;
 121
 122     QVector<T> m_h;
 123     T m_begin;
 124     T m_end;
 125     int m_intervals;
 126
 127 public:
 128     KisCubicSpline() {}
 129     KisCubicSpline(const QList<T_point> &a) {
 130         createSpline(a);
 131     }
 132
 133     /**
 134      * Create new spline and precalculate some values
 135      * for future
 136      *
 137      * @a - base points of the spline
 138      */
 139     void createSpline(const QList<T_point> &a) {
 140         int intervals = m_intervals = a.size() - 1;
 141         int i;
 142         m_begin = a.first().x();
 143         m_end = a.last().x();
 144
 145         m_a.clear();
 146         m_b.resize(intervals);
 147         m_c.clear();
 148         m_d.resize(intervals);
 149         m_h.resize(intervals);
 150
 151         for (i = 0; i < intervals; i++) {
 152             m_h[i] = a[i+1].x() - a[i].x();
 153             m_a.append(a[i].y());
 154         }
 155         m_a.append(a.last().y());
 156
 157
 158         QList<T> tri_b;
 159         QList<T> tri_f;
 160         QList<T> tri_a; /* equals to @tri_c */
 161
 162         for (i = 0; i < intervals - 1; i++) {
 163             tri_b.append(2.*(m_h[i] + m_h[i+1]));
 164
 165             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]));
 166         }
 167         for (i = 1; i < intervals - 1; i++)
 168             tri_a.append(m_h[i]);
 169
 170         if (intervals > 1) {
 171             KisTridiagonalSystem<T> tridia;
 172             m_c = tridia.calculate(tri_a, tri_b, tri_a, tri_f);
 173         }
 174         m_c.prepend(0);
 175         m_c.append(0);
 176
 177         for (i = 0; i < intervals; i++)
 178             m_d[i] = (m_c[i+1] - m_c[i]) / m_h[i];
 179
 180         for (i = 0; i < intervals; i++)
 181             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];
 182     }
 183
 184     /**
 185      * Get value of precalculated spline in the point @x
 186      */
 187     T getValue(T x) const {
 188         T x0;
 189         int i = findRegion(x, x0);
 190         /* TODO: check for asm equivalent */
 191         return m_a[i] +
 192                m_b[i] *(x - x0) +
 193                0.5 * m_c[i] *(x - x0) *(x - x0) +
 194                (1 / 6.0)* m_d[i] *(x - x0) *(x - x0) *(x - x0);
 195     }
 196
 197     T begin() const {
 198         return m_begin;
 199     }
 200
 201     T end() const {
 202         return m_end;
 203     }
 204
 205 protected:
 206
 207     /**
 208      * findRegion - Searches for the region containing @x
 209      * @x0 - out parameter, containing beginning of the region
 210      * @return - index of the region
 211      */
 212     int findRegion(T x, T &x0) const {
 213         int i;
 214         x0 = m_begin;
 215         for (i = 0; i < m_intervals; i++) {
 216             if (x >= x0 && x < x0 + m_h[i])
 217                 return i;
 218             x0 += m_h[i];
 219         }
 220         if (x >= x0) {
 221             x0 -= m_h[m_intervals-1];
 222             return m_intervals - 1;
 223         }
 224
 225         qDebug("X value: %f\n", x);
 226         qDebug("m_begin: %f\n", m_begin);
 227         qDebug("m_end  : %f\n", m_end);
 228         Q_ASSERT_X(0, "findRegion", "X value is outside regions");
 229         /* **never reached** */
 230         return -1;
 231     }
 232 };
 233
 234 static bool pointLessThan(const QPointF &a, const QPointF &b)
 235 {
 236     return a.x() < b.x();
 237 }
 238
 239 struct KisCubicCurve::Data : public QSharedData {
 240     Data() {
 241         init();
 242     }
 243     Data(const Data& data) : QSharedData() {
 244         init();
 245         points = data.points;
 246     }
 247     void init() {
 248         validSpline = false;
 249         validU16Transfer = false;
 250         validFTransfer = false;
 251     }
 252     ~Data() {
 253     }
 254     mutable KisCubicSpline<QPointF, qreal> spline;
 255     QList<QPointF> points;
 256     mutable bool validSpline;
 257     mutable QVector<quint16> u16Transfer;
 258     mutable bool validU16Transfer;
 259     mutable QVector<qreal> fTransfer;
 260     mutable bool validFTransfer;
 261     void updateSpline();
 262     void keepSorted();
 263     qreal value(qreal x);
 264     void invalidate();
 265     template<typename _T_, typename _T2_>
 266     void updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size);
 267 };
 268
 269 void KisCubicCurve::Data::updateSpline()
 270 {
 271     if (validSpline) return;
 272     validSpline = true;
 273     spline.createSpline(points);
 274 }
 275
 276 void KisCubicCurve::Data::invalidate()
 277 {
 278     validSpline = false;
 279     validFTransfer = false;
 280     validU16Transfer = false;
 281 }
 282
 283 void KisCubicCurve::Data::keepSorted()
 284 {
 285     qSort(points.begin(), points.end(), pointLessThan);
 286 }
 287
 288 qreal KisCubicCurve::Data::value(qreal x)
 289 {
 290     updateSpline();
 291     /* Automatically extend non-existing parts of the curve
 292      * (e.g. before the first point) and cut off big y-values
 293      */
 294     x = qBound(spline.begin(), x, spline.end());
 295     qreal y = spline.getValue(x);
 296     return qBound(0.0, y, 1.0);
 297 }
 298
 299 template<typename _T_, typename _T2_>
 300 void KisCubicCurve::Data::updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size)
 301 {
 302     if (!valid || transfer->size() != size) {
 303         if (transfer->size() != size) {
 304             transfer->resize(size);
 305         }
 306         qreal end = 1.0 / (size - 1);
 307         for (int i = 0; i < size; ++i) {
 308             /* Direct uncached version */
 309             _T2_ val = value(i * end) * max;
 310             val = qBound(min, val, max);
 311             (*transfer)[i] = val;
 312         }
 313         valid = true;
 314     }
 315 }
 316
 317 struct KisCubicCurve::Private {
 318     QSharedDataPointer<Data> data;
 319 };
 320
 321 KisCubicCurve::KisCubicCurve() : d(new Private)
 322 {
 323     d->data = new Data;
 324     QPointF p;
 325     p.rx() = 0.0;
 326     p.ry() = 0.0;
 327     d->data->points.append(p);
 328     p.rx() = 1.0;
 329     p.ry() = 1.0;
 330     d->data->points.append(p);
 331 }
 332
 333 KisCubicCurve::KisCubicCurve(const QList<QPointF>& points) : d(new Private)
 334 {
 335     d->data = new Data;
 336     d->data->points = points;
 337     d->data->keepSorted();
 338 }
 339
 340 KisCubicCurve::KisCubicCurve(const KisCubicCurve& curve) : d(new Private(*curve.d))
 341 {
 342 }
 343
 344 KisCubicCurve::~KisCubicCurve()
 345 {
 346     delete d;
 347 }
 348
 349 KisCubicCurve& KisCubicCurve::operator=(const KisCubicCurve & curve)
 350 {
 351     *d = *curve.d;
 352     return *this;
 353 }
 354
 355 bool KisCubicCurve::operator==(const KisCubicCurve& curve) const
 356 {
 357     if (d->data == curve.d->data) return true;
 358     return d->data->points == curve.d->data->points;
 359 }
 360
 361 qreal KisCubicCurve::value(qreal x) const
 362 {
 363     return d->data->value(x);
 364 }
 365
 366 QList<QPointF> KisCubicCurve::points() const
 367 {
 368     return d->data->points;
 369 }
 370
 371 void KisCubicCurve::setPoints(const QList<QPointF>& points)
 372 {
 373     d->data.detach();
 374     d->data->points = points;
 375     d->data->invalidate();
 376 }
 377
 378 void KisCubicCurve::setPoint(int idx, const QPointF& point)
 379 {
 380     d->data.detach();
 381     d->data->points[idx] = point;
 382     d->data->keepSorted();
 383     d->data->invalidate();
 384 }
 385
 386 int KisCubicCurve::addPoint(const QPointF& point)
 387 {
 388     d->data.detach();
 389     d->data->points.append(point);
 390     d->data->keepSorted();
 391     d->data->invalidate();
 392     return d->data->points.indexOf(point);
 393 }
 394
 395 void KisCubicCurve::removePoint(int idx)
 396 {
 397     d->data.detach();
 398     d->data->points.removeAt(idx);
 399     d->data->invalidate();
 400 }
 401
 402 QString KisCubicCurve::toString() const
 403 {
 404     QString sCurve;
 405     foreach(const QPointF & pair, d->data->points) {
 406         sCurve += QString::number(pair.x());
 407         sCurve += ',';
 408         sCurve += QString::number(pair.y());
 409         sCurve += ';';
 410     }
 411     return sCurve;
 412 }
 413
 414 void KisCubicCurve::fromString(const QString& string)
 415 {
 416     QStringList data = string.split(';');
 417
 418     QList<QPointF> points;
 419
 420     foreach(const QString & pair, data) {
 421         if (pair.indexOf(',') > -1) {
 422             QPointF p;
 423             p.rx() = pair.section(',', 0, 0).toDouble();
 424             p.ry() = pair.section(',', 1, 1).toDouble();
 425             points.append(p);
 426         }
 427     }
 428     setPoints(points);
 429 }
 430
 431 QVector<quint16> KisCubicCurve::uint16Transfer(int size) const
 432 {
 433     d->data->updateTransfer<quint16, int>(&d->data->u16Transfer, d->data->validU16Transfer, 0x0, 0xFFFF, size);
 434     return d->data->u16Transfer;
 435 }
 436
 437 QVector<qreal> KisCubicCurve::floatTransfer(int size) const
 438 {
 439     d->data->updateTransfer<qreal, qreal>(&d->data->fTransfer, d->data->validFTransfer, 0.0, 1.0, size);
 440     return d->data->fTransfer;
 441 }