From dd72b368ce415d3c7d4393f60dd4c8147ab46ac5 Mon Sep 17 00:00:00 2001 From: Kyzentun Date: Sat, 3 Jan 2015 05:26:29 -0700 Subject: [PATCH] Fixed spline solving to be mathematically correct after Matt told me column operations are not actually allowed. Attempts to measure performance returned garbage results, hopefully the caching mechanism makes it too fast to measure? --- src/CubicSpline.cpp | 279 ++++++++++++++++++++++++++++---------------- 1 file changed, 179 insertions(+), 100 deletions(-) diff --git a/src/CubicSpline.cpp b/src/CubicSpline.cpp index 1c41cbacc7..c500060ec8 100644 --- a/src/CubicSpline.cpp +++ b/src/CubicSpline.cpp @@ -20,29 +20,39 @@ using std::list; struct SplineSolutionCache { - void solve_diagonals_straight(vector& diagonals); - void solve_diagonals_looped(vector& diagonals); + struct Entry + { + vector diagonals; + vector multiples; + }; + void solve_diagonals_straight(vector& diagonals, vector& multiples); + void solve_diagonals_looped(vector& diagonals, vector& multiples); private: void prep_inner(size_t last, vector& out); - bool find_in_cache(list >& cache, vector& out); - void add_to_cache(list >& cache, vector& out); - list > straight_diagonals; - list > looped_diagonals; + bool find_in_cache(list& cache, vector& outd, vector& outm); + void add_to_cache(list& cache, vector& outd, vector& outm); + list straight_diagonals; + list looped_diagonals; }; size_t const solution_cache_limit= 16; -bool SplineSolutionCache::find_in_cache(list >& cache, vector& out) +bool SplineSolutionCache::find_in_cache(list& cache, vector& outd, vector& outm) { - size_t out_size= out.size(); - for(list >::iterator entry= cache.begin(); + size_t out_size= outd.size(); + for(list::iterator entry= cache.begin(); entry != cache.end(); ++entry) { - if(out_size == entry->size()) + if(out_size == entry->diagonals.size()) { for(size_t i= 0; i < out_size; ++i) { - out[i]= (*entry)[i]; + outd[i]= entry->diagonals[i]; + } + outm.resize(entry->multiples.size()); + for(size_t i= 0; i < entry->multiples.size(); ++i) + { + outm[i]= entry->multiples[i]; } return true; } @@ -50,13 +60,15 @@ bool SplineSolutionCache::find_in_cache(list >& cache, vector >& cache, vector& out) +void SplineSolutionCache::add_to_cache(list& cache, vector& outd, vector& outm) { if(cache.size() >= solution_cache_limit) { cache.pop_back(); } - cache.push_front(out); + cache.push_front(Entry()); + cache.front().diagonals= outd; + cache.front().multiples= outm; } void SplineSolutionCache::prep_inner(size_t last, vector& out) @@ -67,64 +79,138 @@ void SplineSolutionCache::prep_inner(size_t last, vector& out) } } -void SplineSolutionCache::solve_diagonals_straight(vector& diagonals) +void SplineSolutionCache::solve_diagonals_straight(vector& diagonals, vector& multiples) { - if(find_in_cache(straight_diagonals, diagonals)) + if(find_in_cache(straight_diagonals, diagonals, multiples)) { return; } + + // Solution steps: + // Two stages: First, work downwards, zeroing the 1s below each diagonal. + // | 2 1 0 0 | -> | 2 1 0 0 | -> | 2 1 0 0 | -> | 2 1 0 0 | + // | 1 4 1 0 | -> | 0 a 1 0 | -> | 0 d 1 0 | -> | 0 a 1 0 | + // | 0 1 4 1 | -> | 0 1 4 1 | -> | 0 0 b 1 | -> | 0 0 b 1 | + // | 0 0 1 2 | -> | 0 0 1 2 | -> | 0 0 1 2 | -> | 0 0 0 c | + // Second stage: Work upwards, zeroing the 1s above each diagonal. + // V + // | 2 1 0 0 | -> | 2 1 0 0 | -> | 2 0 0 0 | + // | 0 a 1 0 | -> | 0 a 0 0 | -> | 0 a 0 0 | + // | 0 0 b 0 | -> | 0 0 b 0 | -> | 0 0 b 0 | + // | 0 0 0 c | -> | 0 0 0 c | -> | 0 0 0 c | + size_t last= diagonals.size(); diagonals[0]= 2.0f; prep_inner(last-1, diagonals); diagonals[last-1]= 2.0f; - // Operation: Add col[0] * -.5 to col[1] to zero [r0][c1]. + + // Stage one. + // Operation: Add row[0] * -.5 to row[1] to zero [r1][c0]. diagonals[1]-= .5f; + multiples.push_back(.5f); for(size_t i= 1; i < last-1; ++i) { - // Operation: Add col[i] / -[ri][ci] to col[i+1] to zero [ri][ci+1]. - diagonals[i+1]-= 1.0f / diagonals[i]; + // Operation: Add row[i] / -[ri][ci] to row[i+1] to zero [ri+1][ci]. + float const diag_recip= 1.0f / diagonals[i]; + diagonals[i+1]-= diag_recip; + multiples.push_back(diag_recip); + } + // Stage two. + for(size_t i= last-1; i > 0; --i) + { + // Operation: Add row [i] / -[ri][ci] to row[i-1] to zero [ri-1][ci]. + multiples.push_back(1.0f / diagonals[i]); } // Solving finished. - add_to_cache(straight_diagonals, diagonals); + add_to_cache(straight_diagonals, diagonals, multiples); } -void SplineSolutionCache::solve_diagonals_looped(vector& diagonals) +void SplineSolutionCache::solve_diagonals_looped(vector& diagonals, vector& multiples) { - if(find_in_cache(looped_diagonals, diagonals)) + if(find_in_cache(looped_diagonals, diagonals, multiples)) { return; } + + // The steps to solve the system of equations look like this: + // Stage one: Zero the 1s below the diagonals. + // | 4 1 0 0 1 | -> | 4 1 0 0 1 | -> | 4 1 0 0 1 | -> | 4 1 0 0 1 | + // | 1 4 1 0 0 | -> | 0 a 1 0 u | -> | 0 a 1 0 u | -> | 0 a 1 0 u | + // | 0 1 4 1 0 | -> | 0 1 4 1 0 | -> | 0 0 b 1 v | -> | 0 0 b 1 v | + // | 0 0 1 4 1 | -> | 0 0 1 4 1 | -> | 0 0 1 4 1 | -> | 0 0 0 c w | + // | 1 0 0 1 4 | -> | 1 0 0 1 4 | -> | 1 0 0 1 4 | -> | 1 0 0 1 4 | + // V + // | 4 1 0 0 1 | + // | 0 a 1 0 u | + // | 0 0 b 1 v | + // | 0 0 0 c w | + // | 1 0 0 0 d | + // The top of the right column is left unzeroed because it will be changed + // by stage two, nullifying the effect of zeroing it. + // V Stage two: Zero the 1s above the diagonals, starting with the second + // to last row to avoid carrying effects across the left column. + // | 4 1 0 0 1 | -> | 4 1 0 0 1 | -> | 4 0 0 0 z | -> | 4 0 0 0 z | + // | 0 a 1 0 u | -> | 0 a 0 0 y | -> | 0 a 0 0 y | -> | 0 a 0 0 y | + // | 0 0 b 0 x | -> | 0 0 b 0 x | -> | 0 0 b 0 x | -> | 0 0 b 0 x | + // | 0 0 0 c w | -> | 0 0 0 c w | -> | 0 0 0 c w | -> | 0 0 0 c w | + // | 1 0 0 0 d | -> | 1 0 0 0 d | -> | 1 0 0 0 d | -> | 0 0 0 0 f | + // V Stage three: Zero the right column. + // | 4 0 0 0 0 | -> | 4 0 0 0 0 | -> | 4 0 0 0 0 | -> | 4 0 0 0 0 | + // | 0 a 0 0 y | -> | 0 a 0 0 0 | -> | 0 a 0 0 0 | -> | 0 a 0 0 0 | + // | 0 0 b 0 x | -> | 0 0 b 0 x | -> | 0 0 b 0 0 | -> | 0 0 b 0 0 | + // | 0 0 0 c w | -> | 0 0 0 c w | -> | 0 0 0 c w | -> | 0 0 0 c 0 | + // | 0 0 0 0 f | -> | 0 0 0 0 f | -> | 0 0 0 0 f | -> | 0 0 0 0 f | + size_t last= diagonals.size(); diagonals[0]= 4.0f; prep_inner(last, diagonals); + // right_column is sized to not store the diagonal . + vector right_column(diagonals.size()-1, 0.0f); + right_column[0]= 1.0f; + right_column[last-2]= 1.0f; - size_t end= last-1; - size_t stop= end-1; - float cedge= 1.0f; // [ri][cl] - float redge= 1.0f; // [rl][ci] - // The loop stops before end because the case where [ri][cl] == [ri][ci+1] - // needs special handling. - for(size_t i= 0; i < stop; ++i) + // Stage one. + for(size_t i= 0; i < last-2; ++i) { - float next_cedge= 0.0f; // [ri+1][ce] - float next_redge= 0.0f; // [re][ci+1] - // Operation: Add col[i] / -[ri][ci] to col[i+1] to zero [ri][ci+1]. - float diag_recip= 1.0f / diagonals[i]; - diagonals[i+1]-= diag_recip; - next_redge-= redge * diag_recip; // Operation: Add row[i] / -[ri][ci] to row[i+1] to zero [ri+1][ci]. - next_cedge-= cedge * diag_recip; - // Operation: Add col[i] * -(cedge/[ri][ci]) to col[e] to zero cedge. - diagonals[end]-= redge * (cedge / diagonals[i]); - cedge= next_cedge; // Do not use cedge after this point in the loop. - // Operation: Add row[i] * -(redge/[ri][ci]) to row[e] to zero redge. - redge= next_redge; // Do not use redge after this point in the loop. + float const diag_recip= 1.0f / diagonals[i]; + diagonals[i+1]-= diag_recip; + right_column[i+1]-= right_column[i] * diag_recip; + multiples.push_back(diag_recip); } - // [rs][ce] is 1 - cedge, [re][cs] is 1 - redge - // Operation: Add col[s] * -([rs][ce] / [rs][cs]) to col[e] to zero redge. - diagonals[end]-= redge * ((1.0f - cedge) / diagonals[stop]); + // Last step of stage one needs special handling for right_column. + // Operation: Add row[l-2] / [rl-2][cl-2] to row[l-1] to zero [rl-1][cl-2]. + { + float const diag_recip= 1.0f / diagonals[last-2]; + diagonals[last-1]-= right_column[last-2] * diag_recip; + multiples.push_back(diag_recip); + } + // Stage two. + for(size_t i= last-2; i > 0; --i) + { + // Operation: Add row[i] / -[ri][ci] to row[i-1] to zero [ri-1][ci]. + float const diag_recip= 1.0f / diagonals[i]; + right_column[i-1]-= right_column[i] * diag_recip; + multiples.push_back(diag_recip); + } + // Last step of stage two. + { + // Operation: Add row[0] / [r0][c0] to row[l-1] to zero [rl-1][c0]. + float const diag_recip= 1.0f / diagonals[0]; + right_column[0]-= right_column[1] * diag_recip; + multiples.push_back(diag_recip); + } + // Stage three. + size_t const end= last-1; + for(size_t i= 0; i < end; ++i) + { + // Operation: Add row[e] * (right_column[i] / [re][ce]) to row[i] to + // zero right_column[i]. + multiples.push_back(right_column[i] / diagonals[end]); + } + // Solving finished. - add_to_cache(looped_diagonals, diagonals); + add_to_cache(looped_diagonals, diagonals, multiples); } SplineSolutionCache solution_cache; @@ -169,55 +255,46 @@ void CubicSpline::solve_looped() size_t last= m_points.size(); vector results(m_points.size()); vector diagonals(m_points.size()); - solution_cache.solve_diagonals_looped(diagonals); + vector multiples; + solution_cache.solve_diagonals_looped(diagonals, multiples); results[0]= 3 * loop_space_difference( m_points[1].a, m_points[last-1].a, m_spatial_extent); prep_inner(last, results); results[last-1]= 3 * loop_space_difference( m_points[0].a, m_points[last-2].a, m_spatial_extent); - // The steps to solve the system of equations look like this: - // | 4 1 0 0 1 | -> | 4 0 0 0 1 | -> | 4 0 0 0 0 | -> | 4 0 0 0 0 | - // | 1 4 1 0 0 | -> | 0 d 1 0 x | -> | 0 d 1 0 x | -> | 0 d 0 0 x | - // | 0 1 4 1 0 | -> | 0 1 4 1 0 | -> | 0 1 4 1 0 | -> | 0 0 d 1 x | - // | 0 0 1 4 1 | -> | 0 0 1 4 1 | -> | 0 0 1 4 1 | -> | 0 0 1 4 1 | - // | 1 0 0 1 4 | -> | 1 x 0 1 4 | -> | 0 x 0 1 q | -> | 0 x x 1 q | - // V - // | 4 0 0 0 0 | -> | 4 0 0 0 0 | -> | 4 0 0 0 0 | -> | 4 0 0 0 0 | - // | 0 d 0 0 0 | -> | 0 d 0 0 0 | -> | 0 d 0 0 0 | -> | 0 d 0 0 0 | - // | 0 0 d 1 x | -> | 0 0 d 0 x | -> | 0 0 d 0 0 | -> | 0 0 d 0 0 | - // | 0 0 1 d 1 | -> | 0 0 0 d n | -> | 0 0 0 d n | -> | 0 0 0 d 0 | - // | 0 0 x 1 r | -> | 0 0 x n r | -> | 0 0 0 n s | -> | 0 0 0 0 t | - // Each time through the loop performs two of these steps, 4 operations. - // All operations on diagonals are done by the solution cache, because the - // diagonals come out the same for all splines of a given size. - - size_t end= last-1; - size_t stop= end-1; - float cedge= 1.0f; // [ri][cl] - float redge= 1.0f; // [rl][ci] - // The loop stops before end because the case where [ri][cl] == [ri][ci+1] - // needs special handling. - for(size_t i= 0; i < stop; ++i) + // Steps explained in detail in SplineSolutionCache. + // Only the operations on the results column are performed here. + // Stage one. + // SplineSolutionCache's Stage one loop ends at last-2 because it has to + // handle right_column. This does not handle right_column, so the loop + // goes to last-1. + for(size_t i= 0; i < last-1; ++i) { - float next_cedge= 0.0f; // [ri+1][ce] - float next_redge= 0.0f; // [re][ci+1] - // Operation: Add col[i] / -[ri][ci] to col[i+1] to zero [ri][ci+1]. - float diag_recip= 1.0f / diagonals[i]; - next_redge-= redge * diag_recip; - // Operation: Add row[i] / -[ri][ci] to row[i+1] to zero [ri+1][ci]. - results[i+1]-= results[i] * diag_recip; - next_cedge-= cedge * diag_recip; - // Operation: Add col[i] * -(cedge/[ri][ci]) to col[e] to zero cedge. - cedge= next_cedge; // Do not use cedge after this point in the loop. - // Operation: Add row[i] * -(redge/[ri][ci]) to row[e] to zero redge. - results[end]-= results[i] * (redge / diagonals[i]); - redge= next_redge; // Do not use redge after this point in the loop. + // Operation: Add row[i] * -multiples[i] to row[i+1]. + results[i+1]-= results[i] * multiples[i]; } - // [rs][ce] is 1 - cedge, [re][cs] is 1 - redge - // Operation: Add row[s] * -([re][cs] / [rs][cs]) to row[e] to zero redge. - results[end]-= results[stop] * ((1.0f - redge) / diagonals[stop]); - + size_t next_mult= last-1; + // Stage two. + for(size_t i= last-2; i > 0; --i) + { + // Operation: Add row[i] * -multiples[nm] to row[i-1]. + results[i-1]-= results[i] * multiples[next_mult]; + ++next_mult; + } + // Last step of stage two. + // Operation: Add row[0] * -multiples[nm] to row[l-1]. + results[last-1]-= results[0] * multiples[next_mult]; + ++next_mult; + // Stage three. + size_t const end= last-1; + for(size_t i= 0; i < end; ++i) + { + // Operation: Add row[e] * -multiples[nm] to row[i]. + results[i]-= results[end] * multiples[next_mult]; + ++next_mult; + } + // Solving finished. set_results(last, diagonals, results); } @@ -227,28 +304,30 @@ void CubicSpline::solve_straight() size_t last= m_points.size(); vector results(m_points.size()); vector diagonals(m_points.size()); - solution_cache.solve_diagonals_straight(diagonals); + vector multiples; + solution_cache.solve_diagonals_straight(diagonals, multiples); results[0]= 3 * (m_points[1].a - m_points[0].a); prep_inner(last, results); results[last-1]= 3 * loop_space_difference( m_points[last-1].a, m_points[last-2].a, m_spatial_extent); - // The system of equations to be solved looks like this: - // | 2 1 0 0 | = | results[0] | - // | 1 4 1 0 | = | results[1] | - // | 0 1 4 1 | = | results[2] | - // | 0 0 1 2 | = | results[3] | - // Operations are carefully chosen to only modify the values in the - // diagonals and the results, leaving the 1s unchanged. - // All operations on diagonals are done by the solution cache, because the - // diagonals come out the same for all splines of a given size. - // Operation: Add row[0] * -.5 to row[1] to zero [r1][c0]. - results[1]-= results[0] * .5f; - for(size_t i= 1; i < last - 1; ++i) + // Steps explained in detail in SplineSolutionCache. + // Only the operations on the results column are performed here. + // Stage one. + for(size_t i= 0; i < last-1; ++i) { - // Operation: Add row[i] / -[ri][ci] to row[i+1] to zero [ri+1][ci]; - results[i]-= results[i-1] * (1.0f / diagonals[i]); + // Operation: Add row[i] * -multiples[i] to row[i+1]. + results[i+1]-= results[i] * multiples[i]; } + size_t next_mult= last-1; + // Stage two. + for(size_t i= last-1; i > 0; --i) + { + // Operation: Add row[i] * -multiples[nm] to row [i-1]. + results[i-1]-= results[i] * multiples[next_mult]; + ++next_mult; + } + // Solving finished. set_results(last, diagonals, results); }