curvedraw: 7f795f7fecd6 src/curve.cc

curvedraw

view src/curve.cc @ 16:7f795f7fecd6

readme and COPYING

author	John Tsiombikas <nuclear@member.fsf.org>
date	Sun, 20 Dec 2015 10:55:57 +0200
parents	4da693339d99
children

line source

1 /*

2 curvedraw - a simple program to draw curves

5 This program is free software: you can redistribute it and/or modify

6 it under the terms of the GNU General Public License as published by

7 the Free Software Foundation, either version 3 of the License, or

8 (at your option) any later version.

10 This program is distributed in the hope that it will be useful,

11 but WITHOUT ANY WARRANTY; without even the implied warranty of

12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

13 GNU General Public License for more details.

15 You should have received a copy of the GNU General Public License

16 along with this program. If not, see <http://www.gnu.org/licenses/>.

17 */

18 #include <float.h>

19 #include <assert.h>

20 #include <algorithm>

21 #include "curve.h"

23 Curve::Curve(CurveType type)

24 {

25 this->type = type;

26 bbvalid = true;

27 }

29 Curve::Curve(const Vector4 *cp, int numcp, CurveType type)

30 : Curve(type)

31 {

32 this->cp.resize(numcp);

33 for(int i=0; i<numcp; i++) {

34 this->cp[i] = cp[i];

35 }

36 }

38 Curve::Curve(const Vector3 *cp, int numcp, CurveType type)

39 : Curve(type)

40 {

41 this->cp.resize(numcp);

42 for(int i=0; i<numcp; i++) {

43 this->cp[i] = Vector4(cp[i].x, cp[i].y, cp[i].z, 1.0f);

44 }

45 }

47 Curve::Curve(const Vector2 *cp, int numcp, CurveType type)

48 : Curve(type)

49 {

50 this->cp.resize(numcp);

51 for(int i=0; i<numcp; i++) {

52 this->cp[i] = Vector4(cp[i].x, cp[i].y, 0.0f, 1.0f);

53 }

54 }

56 void Curve::set_type(CurveType type)

57 {

58 this->type = type;

59 }

61 CurveType Curve::get_type() const

62 {

63 return type;

64 }

66 void Curve::add_point(const Vector4 &p)

67 {

68 cp.push_back(p);

69 inval_bounds();

70 }

72 void Curve::add_point(const Vector3 &p, float weight)

73 {

74 add_point(Vector4(p.x, p.y, p.z, weight));

75 }

77 void Curve::add_point(const Vector2 &p, float weight)

78 {

79 add_point(Vector4(p.x, p.y, 0.0f, weight));

80 }

82 bool Curve::remove_point(int idx)

83 {

84 if(idx < 0 || idx >= (int)cp.size()) {

85 return false;

86 }

87 cp.erase(cp.begin() + idx);

88 inval_bounds();

89 return true;

90 }

92 void Curve::clear()

93 {

94 cp.clear();

95 inval_bounds();

96 }

98 bool Curve::empty() const

99 {

100 return cp.empty();

101 }

102

103 int Curve::size() const

104 {

105 return (int)cp.size();

106 }

107

108 Vector4 &Curve::operator [](int idx)

109 {

110 inval_bounds();

111 return cp[idx];

112 }

113

114 const Vector4 &Curve::operator [](int idx) const

115 {

116 return cp[idx];

117 }

118

119 const Vector4 &Curve::get_point(int idx) const

120 {

121 return cp[idx];

122 }

123

124 Vector3 Curve::get_point3(int idx) const

125 {

126 return Vector3(cp[idx].x, cp[idx].y, cp[idx].z);

127 }

128

129 Vector2 Curve::get_point2(int idx) const

130 {

131 return Vector2(cp[idx].x, cp[idx].y);

132 }

133

134 float Curve::get_weight(int idx) const

135 {

136 return cp[idx].w;

137 }

138

139 bool Curve::set_point(int idx, const Vector3 &p, float weight)

140 {

141 if(idx < 0 || idx >= (int)cp.size()) {

142 return false;

143 }

144 cp[idx] = Vector4(p.x, p.y, p.z, weight);

145 inval_bounds();

146 return true;

147 }

148

149 bool Curve::set_point(int idx, const Vector2 &p, float weight)

150 {

151 if(idx < 0 || idx >= (int)cp.size()) {

152 return false;

153 }

154 cp[idx] = Vector4(p.x, p.y, 0.0, weight);

155 inval_bounds();

156 return true;

157 }

158

159 bool Curve::set_weight(int idx, float weight)

160 {

161 if(idx < 0 || idx >= (int)cp.size()) {

162 return false;

163 }

164 cp[idx].w = weight;

165 return true;

166 }

167

168 bool Curve::move_point(int idx, const Vector3 &p)

169 {

170 if(idx < 0 || idx >= (int)cp.size()) {

171 return false;

172 }

173 cp[idx] = Vector4(p.x, p.y, p.z, cp[idx].w);

174 inval_bounds();

175 return true;

176 }

177

178 bool Curve::move_point(int idx, const Vector2 &p)

179 {

180 if(idx < 0 || idx >= (int)cp.size()) {

181 return false;

182 }

183 cp[idx] = Vector4(p.x, p.y, 0.0f, cp[idx].w);

184 inval_bounds();

185 return true;

186 }

187

188

189 int Curve::nearest_point(const Vector3 &p) const

190 {

191 int res = -1;

192 float bestsq = FLT_MAX;

193

194 for(size_t i=0; i<cp.size(); i++) {

195 float d = (get_point3(i) - p).length_sq();

196 if(d < bestsq) {

197 bestsq = d;

198 res = i;

199 }

200 }

201 return res;

202 }

203

204 int Curve::nearest_point(const Vector2 &p) const

205 {

206 int res = -1;

207 float bestsq = FLT_MAX;

208

209 for(size_t i=0; i<cp.size(); i++) {

210 float d = (get_point2(i) - p).length_sq();

211 if(d < bestsq) {

212 bestsq = d;

213 res = i;

214 }

215 }

216 return res;

217 }

218

219

220 void Curve::inval_bounds() const

221 {

222 bbvalid = false;

223 }

224

225 void Curve::calc_bounds() const

226 {

227 calc_bbox(&bbmin, &bbmax);

228 bbvalid = true;

229 }

230

231 void Curve::get_bbox(Vector3 *bbmin, Vector3 *bbmax) const

232 {

233 if(!bbvalid) {

234 calc_bounds();

235 }

236 *bbmin = this->bbmin;

237 *bbmax = this->bbmax;

238 }

239

240 void Curve::calc_bbox(Vector3 *bbmin, Vector3 *bbmax) const

241 {

242 if(empty()) {

243 *bbmin = *bbmax = Vector3(0, 0, 0);

244 return;

245 }

246

247 Vector3 bmin = cp[0];

248 Vector3 bmax = bmin;

249 for(size_t i=1; i<cp.size(); i++) {

250 const Vector4 &v = cp[i];

251 for(int j=0; j<3; j++) {

252 if(v[j] < bmin[j]) bmin[j] = v[j];

253 if(v[j] > bmax[j]) bmax[j] = v[j];

254 }

255 }

256 *bbmin = bmin;

257 *bbmax = bmax;

258 }

259

260 void Curve::normalize()

261 {

262 if(!bbvalid) {

263 calc_bounds();

264 }

265

266 Vector3 bsize = bbmax - bbmin;

267 Vector3 boffs = (bbmin + bbmax) * 0.5;

268

269 Vector3 bscale;

270 bscale.x = bsize.x == 0.0f ? 1.0f : 1.0f / bsize.x;

271 bscale.y = bsize.y == 0.0f ? 1.0f : 1.0f / bsize.y;

272 bscale.z = bsize.z == 0.0f ? 1.0f : 1.0f / bsize.z;

273

274 for(size_t i=0; i<cp.size(); i++) {

275 cp[i].x = (cp[i].x - boffs.x) * bscale.x;

276 cp[i].y = (cp[i].y - boffs.y) * bscale.y;

277 cp[i].z = (cp[i].z - boffs.z) * bscale.z;

278 }

279 inval_bounds();

280 }

281

282 Vector3 Curve::proj_point(const Vector3 &p, float refine_thres) const

283 {

284 // first step through the curve a few times and find the nearest of them

285 int num_cp = size();

286 int num_steps = num_cp * 5; // arbitrary number; sounds ok

287 float dt = 1.0f / (float)(num_steps - 1);

288

289 float best_distsq = FLT_MAX;

290 float best_t = 0.0f;

291 Vector3 best_pp;

292

293 float t = 0.0f;

294 for(int i=0; i<num_steps; i++) {

295 Vector3 pp = interpolate(t);

296 float distsq = (pp - p).length_sq();

297 if(distsq < best_distsq) {

298 best_distsq = distsq;

299 best_pp = pp;

300 best_t = t;

301 }

302 t += dt;

303 }

304

305 // refine by gradient descent

306 float dist = best_distsq;

307 t = best_t;

308 dt *= 0.05;

309 for(;;) {

310 float tn = t + dt;

311 float tp = t - dt;

312

313 float dn = (interpolate(tn) - p).length_sq();

314 float dp = (interpolate(tp) - p).length_sq();

315

316 if(fabs(dn - dp) < refine_thres * refine_thres) {

317 break;

318 }

319

320 if(dn < dist) {

321 t = tn;

322 dist = dn;

323 } else if(dp < dist) {

324 t = tp;

325 dist = dp;

326 } else {

327 break; // found the minimum

328 }

329 }

330

331 return interpolate(t);

332 }

333

334 float Curve::distance(const Vector3 &p) const

335 {

336 return (proj_point(p) - p).length();

337 }

338

339 float Curve::distance_sq(const Vector3 &p) const

340 {

341 return fabs((proj_point(p) - p).length_sq());

342 }

343

344

345 Vector3 Curve::interpolate_segment(int a, int b, float t) const

346 {

347 int num_cp = size();

348

349 if(t < 0.0) t = 0.0;

350 if(t > 1.0) t = 1.0;

351

352 Vector4 res;

353 if(type == CURVE_LINEAR || num_cp == 2) {

354 res = lerp(cp[a], cp[b], t);

355 } else {

356 int prev = a <= 0 ? a : a - 1;

357 int next = b >= num_cp - 1 ? b : b + 1;

358

359 if(type == CURVE_HERMITE) {

360 res = catmull_rom_spline(cp[prev], cp[a], cp[b], cp[next], t);

361 } else {

362 res = bspline(cp[prev], cp[a], cp[b], cp[next], t);

363 if(res.w != 0.0f) {

364 res.x /= res.w;

365 res.y /= res.w;

366 res.z /= res.w;

367 }

368 }

369 }

370

371 return Vector3(res.x, res.y, res.z);

372 }

373

374 Vector3 Curve::interpolate(float t) const

375 {

376 if(empty()) {

377 return Vector3(0, 0, 0);

378 }

379

380 int num_cp = (int)cp.size();

381 if(num_cp == 1) {

382 return Vector3(cp[0].x, cp[0].y, cp[0].z);

383 }

384

385 if(t < 0.0) t = 0.0;

386 if(t > 1.0) t = 1.0;

387

388 int idx0 = std::min((int)floor(t * (num_cp - 1)), num_cp - 2);

389 int idx1 = idx0 + 1;

390

391 float dt = 1.0 / (float)(num_cp - 1);

392 float t0 = (float)idx0 * dt;

393 float t1 = (float)idx1 * dt;

394

395 t = (t - t0) / (t1 - t0);

396

397 return interpolate_segment(idx0, idx1, t);

398 }

399

400 Vector2 Curve::interpolate2(float t) const

401 {

402 Vector3 res = interpolate(t);

403 return Vector2(res.x, res.y);

404 }

405

406 Vector3 Curve::operator ()(float t) const

407 {

408 return interpolate(t);

409 }