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() : m_begin(0), m_end(0), m_intervals(0) {}
 130     KisCubicSpline(const QList<T_point> &a) : m_begin(0), m_end(0),
 131       m_intervals(0) {
 132         createSpline(a);
 133     }
 134
 135     /**
 136      * Create new spline and precalculate some values
 137      * for future
 138      *
 139      * @a - base points of the spline
 140      */
 141     void createSpline(const QList<T_point> &a) {
 142         int intervals = m_intervals = a.size() - 1;
 143         int i;
 144         m_begin = a.first().x();
 145         m_end = a.last().x();
 146
 147         m_a.clear();
 148         m_b.resize(intervals);
 149         m_c.clear();
 150         m_d.resize(intervals);
 151         m_h.resize(intervals);
 152
 153         for (i = 0; i < intervals; ++i) {
 154             m_h[i] = a[i+1].x() - a[i].x();
 155             m_a.append(a[i].y());
 156         }
 157         m_a.append(a.last().y());
 158
 159
 160         QList<T> tri_b;
 161         QList<T> tri_f;
 162         QList<T> tri_a; /* equals to @tri_c */
 163
 164         for (i = 0; i < intervals - 1; ++i) {
 165             tri_b.append(2.*(m_h[i] + m_h[i+1]));
 166
 167             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         }
 169         for (i = 1; i < intervals - 1; ++i)
 170             tri_a.append(m_h[i]);
 171
 172         if (intervals > 1) {
 173             KisTridiagonalSystem<T> tridia;
 174             m_c = tridia.calculate(tri_a, tri_b, tri_a, tri_f);
 175         }
 176         m_c.prepend(0);
 177         m_c.append(0);
 178
 179         for (i = 0; i < intervals; ++i)
 180             m_d[i] = (m_c[i+1] - m_c[i]) / m_h[i];
 181
 182         for (i = 0; i < intervals; ++i)
 183             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];
 184     }
 185
 186     /**
 187      * Get value of precalculated spline in the point @x
 188      */
 189     T getValue(T x) const {
 190         T x0;
 191         int i = findRegion(x, x0);
 192         /* TODO: check for asm equivalent */
 193         return m_a[i] +
 194                m_b[i] *(x - x0) +
 195                0.5 * m_c[i] *(x - x0) *(x - x0) +
 196                (1 / 6.0)* m_d[i] *(x - x0) *(x - x0) *(x - x0);
 197     }
 198
 199     T begin() const {
 200         return m_begin;
 201     }
 202
 203     T end() const {
 204         return m_end;
 205     }
 206
 207 protected:
 208
 209     /**
 210      * findRegion - Searches for the region containing @x
 211      * @x0 - out parameter, containing beginning of the region
 212      * @return - index of the region
 213      */
 214     int findRegion(T x, T &x0) const {
 215         int i;
 216         x0 = m_begin;
 217         for (i = 0; i < m_intervals; ++i) {
 218             if (x >= x0 && x < x0 + m_h[i])
 219                 return i;
 220             x0 += m_h[i];
 221         }
 222         if (x >= x0) {
 223             x0 -= m_h[m_intervals-1];
 224             return m_intervals - 1;
 225         }
 226
 227         qDebug("X value: %f\n", x);
 228         qDebug("m_begin: %f\n", m_begin);
 229         qDebug("m_end  : %f\n", m_end);
 230         Q_ASSERT_X(0, "findRegion", "X value is outside regions");
 231         /* **never reached** */
 232         return -1;
 233     }
 234 };
 235
 236 static bool pointLessThan(const QPointF &a, const QPointF &b)
 237 {
 238     return a.x() < b.x();
 239 }
 240
 241 struct KisCubicCurve::Data : public QSharedData {
 242     Data() {
 243         init();
 244     }
 245     Data(const Data& data) : QSharedData() {
 246         init();
 247         points = data.points;
 248     }
 249     void init() {
 250         validSpline = false;
 251         validU16Transfer = false;
 252         validFTransfer = false;
 253     }
 254     ~Data() {
 255     }
 256     mutable KisCubicSpline<QPointF, qreal> spline;
 257     QList<QPointF> points;
 258     mutable bool validSpline;
 259     mutable QVector<quint16> u16Transfer;
 260     mutable bool validU16Transfer;
 261     mutable QVector<qreal> fTransfer;
 262     mutable bool validFTransfer;
 263     void updateSpline();
 264     void keepSorted();
 265     qreal value(qreal x);
 266     void invalidate();
 267     template<typename _T_, typename _T2_>
 268     void updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size);
 269 };
 270
 271 void KisCubicCurve::Data::updateSpline()
 272 {
 273     if (validSpline) return;
 274     validSpline = true;
 275     spline.createSpline(points);
 276 }
 277
 278 void KisCubicCurve::Data::invalidate()
 279 {
 280     validSpline = false;
 281     validFTransfer = false;
 282     validU16Transfer = false;
 283 }
 284
 285 void KisCubicCurve::Data::keepSorted()
 286 {
 287     qSort(points.begin(), points.end(), pointLessThan);
 288 }
 289
 290 qreal KisCubicCurve::Data::value(qreal x)
 291 {
 292     updateSpline();
 293     /* Automatically extend non-existing parts of the curve
 294      * (e.g. before the first point) and cut off big y-values
 295      */
 296     x = qBound(spline.begin(), x, spline.end());
 297     qreal y = spline.getValue(x);
 298     return qBound((qreal)0.0, y, (qreal)1.0);
 299 }
 300
 301 template<typename _T_, typename _T2_>
 302 void KisCubicCurve::Data::updateTransfer(QVector<_T_>* transfer, bool& valid, _T2_ min, _T2_ max, int size)
 303 {
 304     if (!valid || transfer->size() != size) {
 305         if (transfer->size() != size) {
 306             transfer->resize(size);
 307         }
 308         qreal end = 1.0 / (size - 1);
 309         for (int i = 0; i < size; ++i) {
 310             /* Direct uncached version */
 311             _T2_ val = value(i * end) * max;
 312             val = qBound(min, val, max);
 313             (*transfer)[i] = val;
 314         }
 315         valid = true;
 316     }
 317 }
 318
 319 struct KisCubicCurve::Private {
 320     QSharedDataPointer<Data> data;
 321 };
 322
 323 KisCubicCurve::KisCubicCurve() : d(new Private)
 324 {
 325     d->data = new Data;
 326     QPointF p;
 327     p.rx() = 0.0;
 328     p.ry() = 0.0;
 329     d->data->points.append(p);
 330     p.rx() = 1.0;
 331     p.ry() = 1.0;
 332     d->data->points.append(p);
 333 }
 334
 335 KisCubicCurve::KisCubicCurve(const QList<QPointF>& points) : d(new Private)
 336 {
 337     d->data = new Data;
 338     d->data->points = points;
 339     d->data->keepSorted();
 340 }
 341
 342 KisCubicCurve::KisCubicCurve(const KisCubicCurve& curve) : d(new Private(*curve.d))
 343 {
 344 }
 345
 346 KisCubicCurve::~KisCubicCurve()
 347 {
 348     delete d;
 349 }
 350
 351 KisCubicCurve& KisCubicCurve::operator=(const KisCubicCurve & curve)
 352 {
 353     *d = *curve.d;
 354     return *this;
 355 }
 356
 357 bool KisCubicCurve::operator==(const KisCubicCurve& curve) const
 358 {
 359     if (d->data == curve.d->data) return true;
 360     return d->data->points == curve.d->data->points;
 361 }
 362
 363 qreal KisCubicCurve::value(qreal x) const
 364 {
 365     return d->data->value(x);
 366 }
 367
 368 QList<QPointF> KisCubicCurve::points() const
 369 {
 370     return d->data->points;
 371 }
 372
 373 void KisCubicCurve::setPoints(const QList<QPointF>& points)
 374 {
 375     d->data.detach();
 376     d->data->points = points;
 377     d->data->invalidate();
 378 }
 379
 380 void KisCubicCurve::setPoint(int idx, const QPointF& point)
 381 {
 382     d->data.detach();
 383     d->data->points[idx] = point;
 384     d->data->keepSorted();
 385     d->data->invalidate();
 386 }
 387
 388 int KisCubicCurve::addPoint(const QPointF& point)
 389 {
 390     d->data.detach();
 391     d->data->points.append(point);
 392     d->data->keepSorted();
 393     d->data->invalidate();
 394     return d->data->points.indexOf(point);
 395 }
 396
 397 void KisCubicCurve::removePoint(int idx)
 398 {
 399     d->data.detach();
 400     d->data->points.removeAt(idx);
 401     d->data->invalidate();
 402 }
 403
 404 QString KisCubicCurve::toString() const
 405 {
 406     QString sCurve;
 407     QLocale locale;
 408     foreach(const QPointF & pair, d->data->points) {
 409         sCurve += locale.toString(pair.x());
 410         sCurve += '/';
 411         sCurve += QString::number(pair.y());
 412         sCurve += ';';
 413     }
 414     return sCurve;
 415 }
 416
 417 void KisCubicCurve::fromString(const QString& string)
 418 {
 419     QStringList data = string.split(';');
 420
 421     QList<QPointF> points;
 422
 423     foreach(const QString & pair, data) {
 424         if (pair.indexOf('/') > -1) {
 425             QPointF p;
 426             p.rx() = pair.section('/', 0, 0).toDouble();
 427             p.ry() = pair.section('/', 1, 1).toDouble();
 428             points.append(p);
 429         }
 430     }
 431     setPoints(points);
 432 }
 433
 434 QVector<quint16> KisCubicCurve::uint16Transfer(int size) const
 435 {
 436     d->data->updateTransfer<quint16, int>(&d->data->u16Transfer, d->data->validU16Transfer, 0x0, 0xFFFF, size);
 437     return d->data->u16Transfer;
 438 }
 439
 440 QVector<qreal> KisCubicCurve::floatTransfer(int size) const
 441 {
 442     d->data->updateTransfer<qreal, qreal>(&d->data->fTransfer, d->data->validFTransfer, 0.0, 1.0, size);
 443     return d->data->fTransfer;
 444 }