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