curvedraw

annotate src/curve.cc @ 16:7f795f7fecd6

find changesets by author, revision, files, or words in the commit message
readme and COPYING
author John Tsiombikas <nuclear@member.fsf.org>
date Sun, 20 Dec 2015 10:55:57 +0200
parents 4da693339d99
children
rev   line source
nuclear@16 1 /*
nuclear@16 2 curvedraw - a simple program to draw curves
nuclear@16 3 Copyright (C) 2015 John Tsiombikas <nuclear@member.fsf.org>
nuclear@16 4
nuclear@16 5 This program is free software: you can redistribute it and/or modify
nuclear@16 6 it under the terms of the GNU General Public License as published by
nuclear@16 7 the Free Software Foundation, either version 3 of the License, or
nuclear@16 8 (at your option) any later version.
nuclear@16 9
nuclear@16 10 This program is distributed in the hope that it will be useful,
nuclear@16 11 but WITHOUT ANY WARRANTY; without even the implied warranty of
nuclear@16 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
nuclear@16 13 GNU General Public License for more details.
nuclear@16 14
nuclear@16 15 You should have received a copy of the GNU General Public License
nuclear@16 16 along with this program. If not, see <http://www.gnu.org/licenses/>.
nuclear@16 17 */
nuclear@0 18 #include <float.h>
nuclear@12 19 #include <assert.h>
nuclear@5 20 #include <algorithm>
nuclear@0 21 #include "curve.h"
nuclear@0 22
nuclear@0 23 Curve::Curve(CurveType type)
nuclear@0 24 {
nuclear@0 25 this->type = type;
nuclear@12 26 bbvalid = true;
nuclear@12 27 }
nuclear@12 28
nuclear@12 29 Curve::Curve(const Vector4 *cp, int numcp, CurveType type)
nuclear@12 30 : Curve(type)
nuclear@12 31 {
nuclear@12 32 this->cp.resize(numcp);
nuclear@12 33 for(int i=0; i<numcp; i++) {
nuclear@12 34 this->cp[i] = cp[i];
nuclear@12 35 }
nuclear@12 36 }
nuclear@12 37
nuclear@12 38 Curve::Curve(const Vector3 *cp, int numcp, CurveType type)
nuclear@12 39 : Curve(type)
nuclear@12 40 {
nuclear@12 41 this->cp.resize(numcp);
nuclear@12 42 for(int i=0; i<numcp; i++) {
nuclear@12 43 this->cp[i] = Vector4(cp[i].x, cp[i].y, cp[i].z, 1.0f);
nuclear@12 44 }
nuclear@12 45 }
nuclear@12 46
nuclear@12 47 Curve::Curve(const Vector2 *cp, int numcp, CurveType type)
nuclear@12 48 : Curve(type)
nuclear@12 49 {
nuclear@12 50 this->cp.resize(numcp);
nuclear@12 51 for(int i=0; i<numcp; i++) {
nuclear@12 52 this->cp[i] = Vector4(cp[i].x, cp[i].y, 0.0f, 1.0f);
nuclear@12 53 }
nuclear@0 54 }
nuclear@0 55
nuclear@0 56 void Curve::set_type(CurveType type)
nuclear@0 57 {
nuclear@0 58 this->type = type;
nuclear@0 59 }
nuclear@0 60
nuclear@0 61 CurveType Curve::get_type() const
nuclear@0 62 {
nuclear@0 63 return type;
nuclear@0 64 }
nuclear@0 65
nuclear@12 66 void Curve::add_point(const Vector4 &p)
nuclear@12 67 {
nuclear@12 68 cp.push_back(p);
nuclear@12 69 inval_bounds();
nuclear@12 70 }
nuclear@12 71
nuclear@12 72 void Curve::add_point(const Vector3 &p, float weight)
nuclear@12 73 {
nuclear@12 74 add_point(Vector4(p.x, p.y, p.z, weight));
nuclear@12 75 }
nuclear@12 76
nuclear@0 77 void Curve::add_point(const Vector2 &p, float weight)
nuclear@0 78 {
nuclear@12 79 add_point(Vector4(p.x, p.y, 0.0f, weight));
nuclear@0 80 }
nuclear@0 81
nuclear@0 82 bool Curve::remove_point(int idx)
nuclear@0 83 {
nuclear@0 84 if(idx < 0 || idx >= (int)cp.size()) {
nuclear@0 85 return false;
nuclear@0 86 }
nuclear@0 87 cp.erase(cp.begin() + idx);
nuclear@12 88 inval_bounds();
nuclear@0 89 return true;
nuclear@0 90 }
nuclear@0 91
nuclear@12 92 void Curve::clear()
nuclear@0 93 {
nuclear@12 94 cp.clear();
nuclear@12 95 inval_bounds();
nuclear@0 96 }
nuclear@0 97
nuclear@2 98 bool Curve::empty() const
nuclear@2 99 {
nuclear@2 100 return cp.empty();
nuclear@2 101 }
nuclear@2 102
nuclear@2 103 int Curve::size() const
nuclear@0 104 {
nuclear@0 105 return (int)cp.size();
nuclear@0 106 }
nuclear@0 107
nuclear@12 108 Vector4 &Curve::operator [](int idx)
nuclear@12 109 {
nuclear@12 110 inval_bounds();
nuclear@12 111 return cp[idx];
nuclear@12 112 }
nuclear@12 113
nuclear@12 114 const Vector4 &Curve::operator [](int idx) const
nuclear@2 115 {
nuclear@2 116 return cp[idx];
nuclear@2 117 }
nuclear@2 118
nuclear@12 119 const Vector4 &Curve::get_point(int idx) const
nuclear@2 120 {
nuclear@2 121 return cp[idx];
nuclear@2 122 }
nuclear@2 123
nuclear@12 124 Vector3 Curve::get_point3(int idx) const
nuclear@0 125 {
nuclear@12 126 return Vector3(cp[idx].x, cp[idx].y, cp[idx].z);
nuclear@0 127 }
nuclear@0 128
nuclear@12 129 Vector2 Curve::get_point2(int idx) const
nuclear@0 130 {
nuclear@0 131 return Vector2(cp[idx].x, cp[idx].y);
nuclear@0 132 }
nuclear@0 133
nuclear@0 134 float Curve::get_weight(int idx) const
nuclear@0 135 {
nuclear@12 136 return cp[idx].w;
nuclear@12 137 }
nuclear@12 138
nuclear@12 139 bool Curve::set_point(int idx, const Vector3 &p, float weight)
nuclear@12 140 {
nuclear@12 141 if(idx < 0 || idx >= (int)cp.size()) {
nuclear@12 142 return false;
nuclear@12 143 }
nuclear@12 144 cp[idx] = Vector4(p.x, p.y, p.z, weight);
nuclear@12 145 inval_bounds();
nuclear@12 146 return true;
nuclear@0 147 }
nuclear@0 148
nuclear@0 149 bool Curve::set_point(int idx, const Vector2 &p, float weight)
nuclear@0 150 {
nuclear@0 151 if(idx < 0 || idx >= (int)cp.size()) {
nuclear@0 152 return false;
nuclear@0 153 }
nuclear@12 154 cp[idx] = Vector4(p.x, p.y, 0.0, weight);
nuclear@12 155 inval_bounds();
nuclear@0 156 return true;
nuclear@0 157 }
nuclear@0 158
nuclear@0 159 bool Curve::set_weight(int idx, float weight)
nuclear@0 160 {
nuclear@0 161 if(idx < 0 || idx >= (int)cp.size()) {
nuclear@0 162 return false;
nuclear@0 163 }
nuclear@12 164 cp[idx].w = weight;
nuclear@12 165 return true;
nuclear@12 166 }
nuclear@12 167
nuclear@12 168 bool Curve::move_point(int idx, const Vector3 &p)
nuclear@12 169 {
nuclear@12 170 if(idx < 0 || idx >= (int)cp.size()) {
nuclear@12 171 return false;
nuclear@12 172 }
nuclear@12 173 cp[idx] = Vector4(p.x, p.y, p.z, cp[idx].w);
nuclear@12 174 inval_bounds();
nuclear@0 175 return true;
nuclear@0 176 }
nuclear@0 177
nuclear@2 178 bool Curve::move_point(int idx, const Vector2 &p)
nuclear@2 179 {
nuclear@2 180 if(idx < 0 || idx >= (int)cp.size()) {
nuclear@2 181 return false;
nuclear@2 182 }
nuclear@12 183 cp[idx] = Vector4(p.x, p.y, 0.0f, cp[idx].w);
nuclear@12 184 inval_bounds();
nuclear@2 185 return true;
nuclear@2 186 }
nuclear@2 187
nuclear@12 188
nuclear@12 189 int Curve::nearest_point(const Vector3 &p) const
nuclear@0 190 {
nuclear@12 191 int res = -1;
nuclear@12 192 float bestsq = FLT_MAX;
nuclear@12 193
nuclear@12 194 for(size_t i=0; i<cp.size(); i++) {
nuclear@12 195 float d = (get_point3(i) - p).length_sq();
nuclear@12 196 if(d < bestsq) {
nuclear@12 197 bestsq = d;
nuclear@12 198 res = i;
nuclear@12 199 }
nuclear@12 200 }
nuclear@12 201 return res;
nuclear@12 202 }
nuclear@12 203
nuclear@12 204 int Curve::nearest_point(const Vector2 &p) const
nuclear@12 205 {
nuclear@12 206 int res = -1;
nuclear@12 207 float bestsq = FLT_MAX;
nuclear@12 208
nuclear@12 209 for(size_t i=0; i<cp.size(); i++) {
nuclear@12 210 float d = (get_point2(i) - p).length_sq();
nuclear@12 211 if(d < bestsq) {
nuclear@12 212 bestsq = d;
nuclear@12 213 res = i;
nuclear@12 214 }
nuclear@12 215 }
nuclear@12 216 return res;
nuclear@12 217 }
nuclear@12 218
nuclear@12 219
nuclear@12 220 void Curve::inval_bounds() const
nuclear@12 221 {
nuclear@12 222 bbvalid = false;
nuclear@12 223 }
nuclear@12 224
nuclear@12 225 void Curve::calc_bounds() const
nuclear@12 226 {
nuclear@12 227 calc_bbox(&bbmin, &bbmax);
nuclear@12 228 bbvalid = true;
nuclear@12 229 }
nuclear@12 230
nuclear@12 231 void Curve::get_bbox(Vector3 *bbmin, Vector3 *bbmax) const
nuclear@12 232 {
nuclear@12 233 if(!bbvalid) {
nuclear@12 234 calc_bounds();
nuclear@12 235 }
nuclear@12 236 *bbmin = this->bbmin;
nuclear@12 237 *bbmax = this->bbmax;
nuclear@12 238 }
nuclear@12 239
nuclear@12 240 void Curve::calc_bbox(Vector3 *bbmin, Vector3 *bbmax) const
nuclear@12 241 {
nuclear@12 242 if(empty()) {
nuclear@12 243 *bbmin = *bbmax = Vector3(0, 0, 0);
nuclear@12 244 return;
nuclear@12 245 }
nuclear@12 246
nuclear@12 247 Vector3 bmin = cp[0];
nuclear@12 248 Vector3 bmax = bmin;
nuclear@12 249 for(size_t i=1; i<cp.size(); i++) {
nuclear@12 250 const Vector4 &v = cp[i];
nuclear@12 251 for(int j=0; j<3; j++) {
nuclear@12 252 if(v[j] < bmin[j]) bmin[j] = v[j];
nuclear@12 253 if(v[j] > bmax[j]) bmax[j] = v[j];
nuclear@12 254 }
nuclear@12 255 }
nuclear@12 256 *bbmin = bmin;
nuclear@12 257 *bbmax = bmax;
nuclear@12 258 }
nuclear@12 259
nuclear@12 260 void Curve::normalize()
nuclear@12 261 {
nuclear@12 262 if(!bbvalid) {
nuclear@12 263 calc_bounds();
nuclear@12 264 }
nuclear@12 265
nuclear@12 266 Vector3 bsize = bbmax - bbmin;
nuclear@12 267 Vector3 boffs = (bbmin + bbmax) * 0.5;
nuclear@12 268
nuclear@12 269 Vector3 bscale;
nuclear@12 270 bscale.x = bsize.x == 0.0f ? 1.0f : 1.0f / bsize.x;
nuclear@12 271 bscale.y = bsize.y == 0.0f ? 1.0f : 1.0f / bsize.y;
nuclear@12 272 bscale.z = bsize.z == 0.0f ? 1.0f : 1.0f / bsize.z;
nuclear@12 273
nuclear@12 274 for(size_t i=0; i<cp.size(); i++) {
nuclear@12 275 cp[i].x = (cp[i].x - boffs.x) * bscale.x;
nuclear@12 276 cp[i].y = (cp[i].y - boffs.y) * bscale.y;
nuclear@12 277 cp[i].z = (cp[i].z - boffs.z) * bscale.z;
nuclear@12 278 }
nuclear@12 279 inval_bounds();
nuclear@12 280 }
nuclear@12 281
nuclear@13 282 Vector3 Curve::proj_point(const Vector3 &p, float refine_thres) const
nuclear@12 283 {
nuclear@13 284 // first step through the curve a few times and find the nearest of them
nuclear@12 285 int num_cp = size();
nuclear@13 286 int num_steps = num_cp * 5; // arbitrary number; sounds ok
nuclear@13 287 float dt = 1.0f / (float)(num_steps - 1);
nuclear@12 288
nuclear@13 289 float best_distsq = FLT_MAX;
nuclear@13 290 float best_t = 0.0f;
nuclear@13 291 Vector3 best_pp;
nuclear@12 292
nuclear@13 293 float t = 0.0f;
nuclear@13 294 for(int i=0; i<num_steps; i++) {
nuclear@13 295 Vector3 pp = interpolate(t);
nuclear@13 296 float distsq = (pp - p).length_sq();
nuclear@13 297 if(distsq < best_distsq) {
nuclear@13 298 best_distsq = distsq;
nuclear@13 299 best_pp = pp;
nuclear@13 300 best_t = t;
nuclear@13 301 }
nuclear@13 302 t += dt;
nuclear@13 303 }
nuclear@12 304
nuclear@13 305 // refine by gradient descent
nuclear@13 306 float dist = best_distsq;
nuclear@13 307 t = best_t;
nuclear@13 308 dt *= 0.05;
nuclear@13 309 for(;;) {
nuclear@13 310 float tn = t + dt;
nuclear@13 311 float tp = t - dt;
nuclear@12 312
nuclear@13 313 float dn = (interpolate(tn) - p).length_sq();
nuclear@13 314 float dp = (interpolate(tp) - p).length_sq();
nuclear@12 315
nuclear@13 316 if(fabs(dn - dp) < refine_thres * refine_thres) {
nuclear@12 317 break;
nuclear@12 318 }
nuclear@12 319
nuclear@13 320 if(dn < dist) {
nuclear@13 321 t = tn;
nuclear@13 322 dist = dn;
nuclear@13 323 } else if(dp < dist) {
nuclear@13 324 t = tp;
nuclear@13 325 dist = dp;
nuclear@12 326 } else {
nuclear@13 327 break; // found the minimum
nuclear@12 328 }
nuclear@12 329 }
nuclear@13 330
nuclear@13 331 return interpolate(t);
nuclear@12 332 }
nuclear@12 333
nuclear@13 334 float Curve::distance(const Vector3 &p) const
nuclear@13 335 {
nuclear@13 336 return (proj_point(p) - p).length();
nuclear@13 337 }
nuclear@13 338
nuclear@13 339 float Curve::distance_sq(const Vector3 &p) const
nuclear@13 340 {
nuclear@13 341 return fabs((proj_point(p) - p).length_sq());
nuclear@13 342 }
nuclear@13 343
nuclear@13 344
nuclear@12 345 Vector3 Curve::interpolate_segment(int a, int b, float t) const
nuclear@12 346 {
nuclear@12 347 int num_cp = size();
nuclear@12 348
nuclear@12 349 if(t < 0.0) t = 0.0;
nuclear@12 350 if(t > 1.0) t = 1.0;
nuclear@12 351
nuclear@12 352 Vector4 res;
nuclear@12 353 if(type == CURVE_LINEAR || num_cp == 2) {
nuclear@12 354 res = lerp(cp[a], cp[b], t);
nuclear@12 355 } else {
nuclear@12 356 int prev = a <= 0 ? a : a - 1;
nuclear@12 357 int next = b >= num_cp - 1 ? b : b + 1;
nuclear@12 358
nuclear@12 359 if(type == CURVE_HERMITE) {
nuclear@12 360 res = catmull_rom_spline(cp[prev], cp[a], cp[b], cp[next], t);
nuclear@12 361 } else {
nuclear@12 362 res = bspline(cp[prev], cp[a], cp[b], cp[next], t);
nuclear@12 363 if(res.w != 0.0f) {
nuclear@12 364 res.x /= res.w;
nuclear@12 365 res.y /= res.w;
nuclear@13 366 res.z /= res.w;
nuclear@12 367 }
nuclear@12 368 }
nuclear@12 369 }
nuclear@12 370
nuclear@12 371 return Vector3(res.x, res.y, res.z);
nuclear@12 372 }
nuclear@12 373
nuclear@12 374 Vector3 Curve::interpolate(float t) const
nuclear@12 375 {
nuclear@12 376 if(empty()) {
nuclear@12 377 return Vector3(0, 0, 0);
nuclear@0 378 }
nuclear@2 379
nuclear@2 380 int num_cp = (int)cp.size();
nuclear@0 381 if(num_cp == 1) {
nuclear@12 382 return Vector3(cp[0].x, cp[0].y, cp[0].z);
nuclear@0 383 }
nuclear@0 384
nuclear@13 385 if(t < 0.0) t = 0.0;
nuclear@13 386 if(t > 1.0) t = 1.0;
nuclear@13 387
nuclear@2 388 int idx0 = std::min((int)floor(t * (num_cp - 1)), num_cp - 2);
nuclear@0 389 int idx1 = idx0 + 1;
nuclear@0 390
nuclear@0 391 float dt = 1.0 / (float)(num_cp - 1);
nuclear@0 392 float t0 = (float)idx0 * dt;
nuclear@0 393 float t1 = (float)idx1 * dt;
nuclear@0 394
nuclear@0 395 t = (t - t0) / (t1 - t0);
nuclear@0 396
nuclear@12 397 return interpolate_segment(idx0, idx1, t);
nuclear@12 398 }
nuclear@0 399
nuclear@12 400 Vector2 Curve::interpolate2(float t) const
nuclear@12 401 {
nuclear@12 402 Vector3 res = interpolate(t);
nuclear@0 403 return Vector2(res.x, res.y);
nuclear@0 404 }
nuclear@0 405
nuclear@12 406 Vector3 Curve::operator ()(float t) const
nuclear@0 407 {
nuclear@0 408 return interpolate(t);
nuclear@0 409 }