Sketch.cpp 11.3 KB
Newer Older
1
#include <random>
2
#include <utility>
3
4
5

#include "Sketch.h"
#include "Vertex.h"
6
7
#include "CircularArc.h"
#include "LineSegment.h"
8
9
10
11
12
13
14
15
16
17
18
19
#include "Curve.h"
#include "Constraint.h"
#include "Pattern.h"
#include "Constellation.h"
#include "Contour.h"
#include "Star.h"
#include "Branch.h"
#include "Variable.h"
#include "doublen.h"

#include <iostream>

20
double_t Sketch::solve() {
21
22
23
24
25
26
27
28
    // TODO: Use sparse matrix representation for Jacobian, (this is lower priority since typical uses will not produce very large matrices)
    // TODO: Try extracting sub matrix and solving (Certain constraints are known to eliminate variables (Radius, Fixation), if the matrix can be permuted so that a submatrix on the block diagonal is lower-triangular, then a subsystem of equations can be solved independently of the remaining equations. Moreover, if the jacobian can be permuted into a block lower triangular form, subsets of equations can be solved blockwise)
    // TODO: Use variable_feature_size() and equation_feature_size()

    // Matrix and vectors
    Eigen::MatrixXd J = Eigen::MatrixXd::Zero(NumEquations, NumVariables);
    Eigen::VectorXd r = Eigen::VectorXd::Zero(NumEquations);
    Eigen::VectorXd lfs = Eigen::VectorXd::Zero(NumEquations);
Jason Pries's avatar
Jason Pries committed
29

30
    // Tolerances and iteration limits
31
32
    double_t dbl_tol{std::sqrt(NumEquations + 1.0) * DBL_EPSILON};
    double_t flt_tol{std::sqrt(NumEquations + 1.0) * FLT_EPSILON};
Jason Pries's avatar
Jason Pries committed
33

34
    size_t double_bits = (size_t)std::log2(1/DBL_EPSILON);
Jason Pries's avatar
Jason Pries committed
35

36
37
    // Initial residual calculation
    update_linearization(J, r);
Jason Pries's avatar
Jason Pries committed
38

39
    // Setup iteration
40
41
    double_t res_norm{r.norm()};
    double_t res_tol{characteristic_length() * dbl_tol};
42
    long prev_rank(std::min(NumEquations, NumVariables));
43
    for (size_t j = 0; res_norm > res_tol && j < double_bits; ++j) {
44
45
        Eigen::FullPivHouseholderQR<Eigen::MatrixXd> QR = J.fullPivHouseholderQr();
        Eigen::VectorXd delta = QR.solve(r);
Jason Pries's avatar
Jason Pries committed
46

47
48
        if (((J * delta) - r).norm() > res_tol + r.norm() * flt_tol) {
            std::cerr << "Linearized sketch equation could not be solved with an average error per equation less than FLT_EPSILON. Some elements may be over-constrained." << std::endl;
Jason Pries's avatar
Jason Pries committed
49
50
        }

51
52
53
        long iter_rank = QR.rank();
        if (iter_rank < prev_rank) {
            // perturb(delta, characteristic_length()); // perturb update in case of stagnation at inflection point
Jason Pries's avatar
Jason Pries committed
54
        }
55
56
57
58
59
        prev_rank = iter_rank;

        update_variables(delta);
        update_linearization(J, r);

60
61
        double_t iter_norm{r.norm()};
        double_t iter_tol{res_tol};
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
        if (iter_norm > res_norm) { // Line search
            delta *= -0.5;

            size_t line_search_iter{0};
            while (iter_norm > res_norm && iter_norm > iter_tol && line_search_iter != double_bits) {
                update_variables(delta);
                update_linearization(J, r);

                iter_norm = r.norm();
                iter_tol = dbl_tol * characteristic_length();
                if (iter_norm > res_norm) {
                    delta *= 0.5;
                    line_search_iter += 1;

                    if (line_search_iter == double_bits) {
77
78
                        std::cerr << "Residual norm was not reduced for any step-size greater than DBL_EPSILON" << std::endl;
                        break;
79
80
81
                    }
                }
            }
Jason Pries's avatar
Jason Pries committed
82
        }
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
        res_norm = iter_norm;
        res_tol = iter_tol;
    }

    return res_norm;
}

void Sketch::update_linearization(Eigen::MatrixXd &J, Eigen::VectorXd &r) const {
    J.setZero();

    for (size_t i = 0; i != Verticies.size(); ++i) {
        Verticies[i]->update(J, r);
    }

    for (size_t i = 0; i != Curves.size(); ++i) {
        Curves[i]->update(J, r);
    }

    for (size_t i = 0; i != Constraints.size(); ++i) {
        Constraints[i]->update(J, r);
    }
}
Jason Pries's avatar
Jason Pries committed
105

106
107
108
109
110
void Sketch::update_variables(Eigen::VectorXd &delta) const {
    for (size_t i = 0; i != Variables.size(); ++i) {
        Variables[i]->update(delta);
    }
}
Jason Pries's avatar
Jason Pries committed
111

112
void Sketch::perturb(Eigen::VectorXd &delta, double_t scale) const {
113
114
115
116
117
118
119
120
121
122
    // Randomly perturb the delta vector to prevent stagnation at inflection point
    std::random_device rdev;
    std::mt19937 rgen(rdev());
    std::uniform_real_distribution<double> rdist = std::uniform_real_distribution<double>(-DBL_EPSILON * scale, DBL_EPSILON * scale);
    for (size_t i = 0; i != Variables.size(); ++i) {
        size_t j = Variables[i]->get_index();
        delta(j) += rdist(rgen) * scale;
    }
}

123
124
double_t Sketch::characteristic_length() const {
    double_t scale{0.0};
125
126
127
128
129
130

    if (Curves.size() > 0) {
        for (size_t i = 0; i != Curves.size(); ++i) {
            scale = std::max(scale, Curves[i]->length());
        }
    } else if (Verticies.size() > 0) {
131
132
133
134
        double_t xmax = DBL_MIN;
        double_t xmin = DBL_MAX;
        double_t ymax = DBL_MIN;
        double_t ymin = DBL_MAX;
135
136
137
138
139
140
141
142
143
144
145
        for (size_t i = 0; i != Verticies.size(); ++i) {
            xmax = std::max(xmax, Verticies[i]->x());
            xmin = std::min(xmin, Verticies[i]->x());

            ymax = std::max(ymax, Verticies[i]->y());
            ymin = std::max(ymin, Verticies[i]->y());
        }

        scale = std::max(xmax - xmin, ymax-ymin);
        if(scale == 0) {
            scale = std::max(xmax,ymax);
Jason Pries's avatar
Jason Pries committed
146
147
        }
    }
148
149

    return scale;
Jason Pries's avatar
Jason Pries committed
150
151
152
153
154
}

bool Sketch::build() {
    Constellation c = Constellation(this);

155
156
    Boundary = c.boundary();
    bool success = (Boundary->size() > 0);
Jason Pries's avatar
Jason Pries committed
157

158
159
    Contours = c.contours();
    success = success && (Contours.size() > 0);
Jason Pries's avatar
Jason Pries committed
160
161
162
163

    return success;
}

164
void Sketch::add_element(std::shared_ptr<Vertex> v) {
Jason Pries's avatar
Jason Pries committed
165
166
167
    add_element(Verticies, v);
}

168
void Sketch::add_element(std::shared_ptr<Curve> c) {
Jason Pries's avatar
Jason Pries committed
169
170
171
    add_element(Curves, c);
}

172
void Sketch::add_element(std::shared_ptr<Constraint> c) {
Jason Pries's avatar
Jason Pries committed
173
174
175
    add_element(Constraints, c);
}

176
void Sketch::add_element(std::shared_ptr<Pattern> p) {
Jason Pries's avatar
Jason Pries committed
177
    add_element(Patterns, p);
178
    p->register_elements(this);
Jason Pries's avatar
Jason Pries committed
179
180
181
}

template<>
182
void Sketch::save_as<SaveMethod::Rasterize>(std::string path, std::string file_name) const {
Jason Pries's avatar
Jason Pries committed
183
184
185
    /*
        This is a stub for visualization
    */
186

187
188
    if (!std::experimental::filesystem::exists(path)) {
        std::experimental::filesystem::create_directories(path);
189
190
    }

Jason Pries's avatar
Jason Pries committed
191
    std::fstream fs;
192
193
    fs.open(path + file_name + ".oesk", std::fstream::out);

194
    for (size_t i = 0; i != Curves.size(); ++i) {
195
        if (!Curves[i]->for_construction()) {
196
            for (size_t j = 0; j != 11; ++j) {
197
198
                double2 v = Curves[i]->point(j / 10.0);
                fs << v.X << ',' << v.Y << '\n';
Jason Pries's avatar
Jason Pries committed
199
            }
200
            fs << "NaN" << ',' << "NaN" << '\n';
Jason Pries's avatar
Jason Pries committed
201
202
        }
    }
203

Jason Pries's avatar
Jason Pries committed
204
    fs.close();
205
206
}

207
std::vector<PeriodicBoundaryPair> Sketch::select_periodic_boundary_pairs(std::shared_ptr<Vertex const> v0, double_t angle) const {
208
    std::vector<PeriodicBoundaryPair> cbp;
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229

    std::vector<std::shared_ptr<Curve const>> bc = Boundary->curves();

    while (bc.begin() != bc.end()) {
        auto c0 = bc.begin();
        auto c1 = --bc.end();

        MatchOrientation m{MatchOrientation::None};
        while (c0 != c1) {
            auto cc0 = (*c0).get();
            auto cc1 = (*c1).get();
            m = (*c1)->is_identical(*c0, v0, angle);

            if (m == MatchOrientation::None) {
                m = (*c0)->is_identical(*c1, v0, angle);
                if (m != MatchOrientation::None) {
                    std::swap(c0, c1);
                }
            }

            if (m == MatchOrientation::Forward) {
230
                cbp.emplace_back(*c0, *c1, true, v0, angle);
231
232
                break;
            } else if (m == MatchOrientation::Reverse) {
233
                cbp.emplace_back(*c0, *c1, false, v0, angle);
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
                break;
            } else {
                ++c0;
            }
        }

        if (m != MatchOrientation::None) { // if match found, erase c1
            bc.erase(c0);
            bc.erase(c1);
        } else {
            bc.erase(c1);
        }
    }

    return cbp;
}

251
std::vector<std::shared_ptr<Curve const>> Sketch::select_radial_boundary(std::shared_ptr<Vertex const> v0, double_t radius) const {
252
253
254
255
256
257
258
259
260
261
262
263
    std::vector<std::shared_ptr<Curve const>> rb;

    std::vector<std::shared_ptr<Curve const>> bc = Boundary->curves();

    auto test_circle = std::make_shared<CircularArc>(std::make_shared<Vertex>(), std::make_shared<Vertex>(), v0, radius);
    for (std::shared_ptr<Curve const> c : bc) {
        if (test_circle->is_coincident(c)) {
            rb.push_back(c);
        }
    }

    return rb;
264
265
}

266
std::shared_ptr<Vertex const> Sketch::select_periodic_vertex(std::shared_ptr<Vertex const> v, std::shared_ptr<Vertex const> origin, double_t angle) const {
267
268
269
270
271
272
273
274
    for(auto const &v_test : Verticies) {
        if (v_test->is_identical(v, origin, angle)) {
            return v_test;
        }
    }

    std::shared_ptr<Vertex const> nullv;
    return nullv;
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
}

std::vector<std::shared_ptr<Contour const>> Sketch::select_contours(BoundingBall bb) const {
    /*
     * Default behavior ensures a selected contour has all curve bounding balls within the specified bounding ball
     *
     * TODO: Options
     *          1) (Weakest) Select contours whose individual curves bounding ball centers are contained within the specified ball
     *          2) (Weaker) Select contours whose overall bounding ball center is contained within the specified ball, utilizing operator+(BoundingBall,BoundingBall) to construct contour bounding ball
     *          3) (Stronger) Select contours whose overall bounding ball is contained within the specified ball, contour bounding ball constructed as 2)
     */

    std::vector<std::shared_ptr<Contour const>> output;

    for (size_t i = 0; i != Contours.size(); ++i) {
        if (Contours[i]->is_contained_by(bb)) {
            output.push_back(Contours[i]);
        }
    }

    return output;
}

std::vector<std::shared_ptr<Contour const>> Sketch::select_contours(std::vector<BoundingBall> bb_in, std::vector<BoundingBall> bb_out) const {
    std::vector<std::shared_ptr<Contour const>> output;

    for (size_t i = 0; i != Contours.size(); ++i) {
        // Determine if Contours[i] is contained by any bounding ball in bb_in
        bool is_selected{false};
        size_t j = 0;
        while (!is_selected && j != bb_in.size()) {
            is_selected |= Contours[i]->is_contained_by(bb_in[j]);
            ++j;
        }

        if (is_selected) { // If it is, determine if Contours[i] is contained by any bounding ball in bb_out
            j = 0;
            while (is_selected && j != bb_out.size()) {
                is_selected &= !(Contours[i]->is_contained_by(bb_out[j]));
                ++j;
            }
        }

        if (is_selected) { // if (Contours[i] in any of bb_in) && !(Contours[i] in any of bb_out), then select Contour[i]
            output.push_back(Contours[i]);
        }
    }

    return output;
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
}

std::shared_ptr<Vertex> Sketch::select_nearest_vertex(double_t x, double_t y) const {
    std::shared_ptr<Vertex> nearest;
    double_t dmin = DBL_MAX;

    for (auto v : Verticies) {
        double_t d = hypot(x - v->x(), y - v->y());
        if (d < dmin) {
            nearest = v;
            dmin = d;
        }
    }

    return nearest;
Jason Pries's avatar
Jason Pries committed
339
}