mslang.cpp 28.6 KB
Newer Older
1 2 3
#include "Algorithms/Vektor.h"
#include "Algorithms/Quaternion.h"
#include "AppTypes/Tenzor.h"
kraus's avatar
kraus committed
4 5
#include "Physics/Physics.h"

6 7
#include <boost/regex.hpp>

8 9
#include <gsl/gsl_rng.h>

10 11
#include <iostream>
#include <string>
12 13
#include <fstream>
#include <streambuf>
14
#include <cstdlib>
15
#include <cmath>
16 17 18 19 20 21

std::string UDouble("([0-9]+\\.?[0-9]*([Ee][+-]?[0-9]+)?)");
std::string Double("(-?[0-9]+\\.?[0-9]*([Ee][+-]?[0-9]+)?)");
std::string UInt("([0-9]+)");
std::string FCall("([a-z]*)\\((.*)");

kraus's avatar
kraus committed
22 23
typedef std::string::iterator iterator;

24
struct BoundingBox {
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
    Vector_t center_m;
    double width_m;
    double height_m;

    BoundingBox():
        center_m(0.0),
        width_m(0.0),
        height_m(0.0)
    { }

    BoundingBox(const Vector_t &llc,
                const Vector_t &urc):
        center_m(0.5 * (llc + urc)),
        width_m(urc[0] - llc[0]),
        height_m(urc[1] - llc[1])
    { }

42
    bool isInside(const Vector_t &X) const {
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
        if (2 * std::abs(X[0] - center_m[0]) <= width_m &&
            2 * std::abs(X[1] - center_m[1]) <= height_m)
            return true;

        return false;
    }

    virtual void writeGnuplot(std::ofstream &out) const {
        std::vector<Vector_t> pts({Vector_t(center_m[0] + 0.5 * width_m, center_m[1] + 0.5 * height_m, 0),
                                   Vector_t(center_m[0] - 0.5 * width_m, center_m[1] + 0.5 * height_m, 0),
                                   Vector_t(center_m[0] - 0.5 * width_m, center_m[1] - 0.5 * height_m, 0),
                                   Vector_t(center_m[0] + 0.5 * width_m, center_m[1] - 0.5 * height_m, 0)});
        unsigned int width = out.precision() + 8;
        for (unsigned int i = 0; i < 5; ++ i) {
            Vector_t & pt = pts[i % 4];

            out << std::setw(width) << pt[0]
                << std::setw(width) << pt[1]
                << std::endl;
        }
        out << std::endl;
    }
};

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
struct AffineTransformation: public Tenzor<double, 3> {
    AffineTransformation(const Vector_t& row0,
                         const Vector_t& row1):
        Tenzor(row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], 0.0, 0.0, 1.0) {
    }

    AffineTransformation():
        Tenzor(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0) { }

    AffineTransformation getInverse() const {
        AffineTransformation Ret;
        double det = (*this)(0, 0) * (*this)(1, 1) - (*this)(1, 0) * (*this)(0, 1);

        Ret(0, 0) = (*this)(1, 1) / det;
        Ret(1, 0) = -(*this)(1, 0) / det;
        Ret(0, 1) = -(*this)(0, 1) / det;
        Ret(1, 1) = (*this)(0, 0) / det;

        Ret(0, 2) = - Ret(0, 0) * (*this)(0, 2) - Ret(0, 1) * (*this)(1, 2);
        Ret(1, 2) = - Ret(1, 0) * (*this)(0, 2) - Ret(1, 1) * (*this)(1, 2);
        Ret(2, 2) = 1.0;

        return Ret;
    }

    Vector_t getOrigin() const {
        return Vector_t(-(*this)(0, 2), -(*this)(1, 2), 0.0);
    }

    double getAngle() const {
        return atan2((*this)(1, 0), (*this)(0, 0));
    }

    Vector_t transformTo(const Vector_t &v) const {
        const Tenzor<double, 3> &A = *static_cast<const Tenzor<double, 3>* >(this);
        Vector_t b(v[0], v[1], 1.0);
103 104 105
        Vector_t w = dot(A, b);

        return Vector_t(w[0], w[1], 0.0);
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
    }

    Vector_t transformFrom(const Vector_t &v) const {
        AffineTransformation inv = getInverse();
        return inv.transformTo(v);
    }

    AffineTransformation mult(const AffineTransformation &B) {
        AffineTransformation Ret;
        const Tenzor<double, 3> &A = *static_cast<const Tenzor<double, 3> *>(this);
        const Tenzor<double, 3> &BTenz = *static_cast<const Tenzor<double, 3> *>(&B);
        Tenzor<double, 3> &C = *static_cast<Tenzor<double, 3> *>(&Ret);

        C = dot(A, BTenz);

        return Ret;
    }
};

kraus's avatar
kraus committed
125
struct Base;
kraus's avatar
kraus committed
126

kraus's avatar
kraus committed
127 128
struct Function {
    virtual ~Function() {};
kraus's avatar
kraus committed
129

130
    virtual void print(int indent) = 0;
kraus's avatar
kraus committed
131
    virtual void apply(std::vector<Base*> &bfuncs) = 0;
132 133
};

kraus's avatar
kraus committed
134 135 136
bool parse(iterator &it, const iterator &end, Function* &fun);

struct Base: public Function {
137
    AffineTransformation trafo_m;
138
    BoundingBox bb_m;
139

kraus's avatar
kraus committed
140
    Base():
141
        trafo_m()
kraus's avatar
kraus committed
142 143
    { }

kraus's avatar
kraus committed
144
    virtual Base* clone() = 0;
145
    virtual void writeGnuplot(std::ofstream &out) const = 0;
146
    virtual void computeBoundingBox() = 0;
147
    virtual bool isInside(const Vector_t &R) const = 0;
kraus's avatar
kraus committed
148 149
};

kraus's avatar
kraus committed
150 151 152
struct Rectangle: public Base {
    double width_m;
    double height_m;
153

kraus's avatar
kraus committed
154 155 156 157
    Rectangle():
        Base(),
        width_m(0.0),
        height_m(0.0)
kraus's avatar
kraus committed
158 159
    { }

kraus's avatar
kraus committed
160
    virtual ~Rectangle() { }
kraus's avatar
kraus committed
161

162 163 164
    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
kraus's avatar
kraus committed
165
        Vector_t origin = trafo_m.getOrigin();
166
        double angle = trafo_m.getAngle() * Physics::rad2deg;
167
        std::cout << indent << "rectangle, \n"
kraus's avatar
kraus committed
168 169
                  << indent2 << "w: " << width_m << ", \n"
                  << indent2 << "h: " << height_m << ", \n"
kraus's avatar
kraus committed
170
                  << indent2 << "origin: " << origin[0] << ", " << origin[1] << ",\n"
171 172 173 174
                  << indent2 << "angle: " << angle << "\n"
                  << indent2 << trafo_m(0, 0) << "\t" << trafo_m(0, 1) << "\t" << trafo_m(0, 2) << "\n"
                  << indent2 << trafo_m(1, 0) << "\t" << trafo_m(1, 1) << "\t" << trafo_m(1, 2) << "\n"
                  << indent2 << trafo_m(2, 0) << "\t" << trafo_m(2, 1) << "\t" << trafo_m(2, 2) << std::endl;
kraus's avatar
kraus committed
175 176
    }

177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
    virtual void computeBoundingBox() {
        std::vector<Vector_t> corners({Vector_t(0.5 * width_m, 0.5 * height_m, 0),
                    Vector_t(-0.5 * width_m, 0.5 * height_m, 0),
                    Vector_t(-0.5 * width_m, -0.5 * height_m, 0),
                    Vector_t(0.5 * width_m, -0.5 * height_m, 0)});

        for (Vector_t &v: corners) {
            v = trafo_m.transformFrom(v);
        }

        Vector_t llc = corners[0], urc = corners[0];
        for (unsigned int i = 1; i < 4; ++ i) {
            if (corners[i][0] < llc[0]) llc[0] = corners[i][0];
            else if (corners[i][0] > urc[0]) urc[0] = corners[i][0];

            if (corners[i][1] < llc[1]) llc[1] = corners[i][1];
            else if (corners[i][1] > urc[1]) urc[1] = corners[i][1];
        }

        bb_m = BoundingBox(llc, urc);
    }

199 200 201 202 203 204 205 206 207 208
    virtual bool isInside(const Vector_t &R) const {
        if (!bb_m.isInside(R)) return false;

        Vector_t X = trafo_m.transformTo(R);
        if (2 * std::abs(X[0]) <= width_m &&
            2 * std::abs(X[1]) <= height_m) return true;

        return false;
    }

209
    virtual void writeGnuplot(std::ofstream &out) const {
kraus's avatar
kraus committed
210 211 212 213
        std::vector<Vector_t> pts({Vector_t(0.5 * width_m, 0.5 * height_m, 0),
                                   Vector_t(-0.5 * width_m, 0.5 * height_m, 0),
                                   Vector_t(-0.5 * width_m, -0.5 * height_m, 0),
                                   Vector_t(0.5 * width_m, -0.5 * height_m, 0)});
214 215 216
        unsigned int width = out.precision() + 8;
        for (unsigned int i = 0; i < 5; ++ i) {
            Vector_t pt = pts[i % 4];
kraus's avatar
kraus committed
217
            pt = trafo_m.transformFrom(pt);
218 219 220 221 222 223

            out << std::setw(width) << pt[0]
                << std::setw(width) << pt[1]
                << std::endl;
        }
        out << std::endl;
224 225

        bb_m.writeGnuplot(out);
226 227
    }

kraus's avatar
kraus committed
228
    virtual void apply(std::vector<Base*> &bfuncs) {
kraus's avatar
kraus committed
229 230 231
        bfuncs.push_back(this->clone());
    }

kraus's avatar
kraus committed
232 233 234 235 236
    virtual Base* clone() {
        Rectangle *rect = new Rectangle;
        rect->width_m = width_m;
        rect->height_m = height_m;
        rect->trafo_m = trafo_m;
kraus's avatar
kraus committed
237 238

        return rect;
239 240 241
    }

    static
kraus's avatar
kraus committed
242
    bool parse_detail(iterator &it, const iterator &end, Function* fun) {
243 244 245 246 247 248
        std::string str(it, end);
        boost::regex argumentList(UDouble + "," + UDouble + "(\\).*)");
        boost::smatch what;

        if (!boost::regex_match(str, what, argumentList)) return false;

kraus's avatar
kraus committed
249 250 251
        Rectangle *rect = static_cast<Rectangle*>(fun);
        rect->width_m = atof(std::string(what[1]).c_str());
        rect->height_m = atof(std::string(what[3]).c_str());
252 253 254 255 256 257 258 259 260 261

        std::string fullMatch = what[0];
        std::string rest = what[5];
        it += (fullMatch.size() - rest.size() + 1);

        return true;

    }
};

kraus's avatar
kraus committed
262 263 264
struct Ellipse: public Base {
    double width_m;
    double height_m;
265

kraus's avatar
kraus committed
266 267 268
    Ellipse():
        width_m(0.0),
        height_m(0.0)
kraus's avatar
kraus committed
269 270
    { }

kraus's avatar
kraus committed
271
    virtual ~Ellipse() { }
kraus's avatar
kraus committed
272

273 274 275
    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
kraus's avatar
kraus committed
276
        Vector_t origin = trafo_m.getOrigin();
277
        double angle = trafo_m.getAngle() * Physics::rad2deg;
278
        std::cout << indent << "ellipse, \n"
kraus's avatar
kraus committed
279 280
                  << indent2 << "w: " << width_m << ", \n"
                  << indent2 << "h: " << height_m << ", \n"
kraus's avatar
kraus committed
281
                  << indent2 << "origin: " << origin[0] << ", " << origin[1] << ",\n"
282
                  << indent2 << "angle: " << angle << "\n"
283 284 285
                  << indent2 << std::setw(14) << trafo_m(0, 0) << std::setw(14) << trafo_m(0, 1) << std::setw(14) << trafo_m(0, 2) << "\n"
                  << indent2 << std::setw(14) << trafo_m(1, 0) << std::setw(14) << trafo_m(1, 1) << std::setw(14) << trafo_m(1, 2) << "\n"
                  << indent2 << std::setw(14) << trafo_m(2, 0) << std::setw(14) << trafo_m(2, 1) << std::setw(14) << trafo_m(2, 2)
286
                  << std::endl;
kraus's avatar
kraus committed
287 288
    }

289 290 291 292 293 294 295 296
    virtual void writeGnuplot(std::ofstream &out) const {
        const unsigned int N = 101;
        const double dp = Physics::two_pi / (N - 1);
        const unsigned int colwidth = out.precision() + 8;

        double phi = 0;
        for (unsigned int i = 0; i < N; ++ i, phi += dp) {
            Vector_t pt(0.0);
kraus's avatar
kraus committed
297 298 299
            pt[0] = std::copysign(sqrt(std::pow(height_m * width_m * 0.25, 2) /
                                       (std::pow(height_m * 0.5, 2) +
                                        std::pow(width_m * 0.5 * tan(phi), 2))),
300 301
                                  cos(phi));
            pt[1] = pt[0] * tan(phi);
kraus's avatar
kraus committed
302
            pt = trafo_m.transformFrom(pt);
303 304 305 306 307 308 309

            out << std::setw(colwidth) << pt[0]
                << std::setw(colwidth) << pt[1]
                << std::endl;
        }
        out << std::endl;

310
        bb_m.writeGnuplot(out);
311 312
    }

kraus's avatar
kraus committed
313
    virtual void apply(std::vector<Base*> &bfuncs) {
kraus's avatar
kraus committed
314 315 316
        bfuncs.push_back(this->clone());
    }

kraus's avatar
kraus committed
317 318 319 320 321
    virtual Base* clone() {
        Ellipse *elps = new Ellipse;
        elps->width_m = width_m;
        elps->height_m = height_m;
        elps->trafo_m = trafo_m;
kraus's avatar
kraus committed
322 323

        return elps;
324 325
    }

326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353
    virtual void computeBoundingBox() {
        Vector_t llc(0.0), urc(0.0);
        const Vector_t e_x(1.0, 0.0, 0.0), e_y(0.0, 1.0, 0.0);
        const Vector_t center = trafo_m.transformFrom(Vector_t(0.0));
        const Vector_t e_xp = trafo_m.transformFrom(e_x) - center;
        const Vector_t e_yp = trafo_m.transformFrom(e_y) - center;
        const double &M11 = e_xp[0];
        const double &M12 = e_yp[0];
        const double &M21 = e_xp[1];
        const double &M22 = e_yp[1];

        double t = atan2(height_m * M12, width_m * M11);
        double halfwidth = 0.5 * (M11 * width_m * cos(t) +
                                  M12 * height_m * sin(t));
        llc[0] = center[0] - std::abs(halfwidth);
        urc[0] = center[0] + std::abs(halfwidth);

        t = atan2(height_m * M22, width_m * M21);

        double halfheight = 0.5 * (M21 * width_m * cos(t) +
                                   M22 * height_m * sin(t));

        llc[1] = center[1] - std::abs(halfheight);
        urc[1] = center[1] + std::abs(halfheight);

        bb_m = BoundingBox(llc, urc);
    }

354 355 356 357 358 359 360 361 362 363
    virtual bool isInside(const Vector_t &R) const {
        if (!bb_m.isInside(R)) return false;

        Vector_t X = trafo_m.transformTo(R);
        if (4 * (std::pow(X[0] / width_m, 2) + std::pow(X[1] / height_m, 2)) <= 1)
            return true;

        return false;
    }

364
    static
kraus's avatar
kraus committed
365
    bool parse_detail(iterator &it, const iterator &end, Function* fun) {
366 367 368 369 370 371
        std::string str(it, end);
        boost::regex argumentList(UDouble + "," + UDouble + "(\\).*)");
        boost::smatch what;

        if (!boost::regex_match(str, what, argumentList)) return false;

kraus's avatar
kraus committed
372 373 374
        Ellipse *elps = static_cast<Ellipse*>(fun);
        elps->width_m = atof(std::string(what[1]).c_str());
        elps->height_m = atof(std::string(what[3]).c_str());
375 376 377 378 379 380 381 382 383

        std::string fullMatch = what[0];
        std::string rest = what[5];
        it += (fullMatch.size() - rest.size() + 1);

        return true;
    }
};

kraus's avatar
kraus committed
384 385 386 387 388
struct Repeat: public Function {
    Function* func_m;
    unsigned int N_m;
    double shiftx_m;
    double shifty_m;
389

kraus's avatar
kraus committed
390 391
    virtual ~Repeat() {
        delete func_m;
kraus's avatar
kraus committed
392 393
    }

394 395 396 397
    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
        std::cout << indent << "repeat, " << std::endl;
kraus's avatar
kraus committed
398
        func_m->print(indentwidth + 8);
399
        std::cout << ",\n"
kraus's avatar
kraus committed
400 401 402
                  << indent2 << "N: " << N_m << ", \n"
                  << indent2 << "dx: " << shiftx_m << ", \n"
                  << indent2 << "dy: " << shifty_m;
kraus's avatar
kraus committed
403 404
    }

kraus's avatar
kraus committed
405
    virtual void apply(std::vector<Base*> &bfuncs) {
406 407
        AffineTransformation shift(Vector_t(1.0, 0.0, -shiftx_m),
                                   Vector_t(0.0, 1.0, -shifty_m));
408

kraus's avatar
kraus committed
409
        func_m->apply(bfuncs);
kraus's avatar
kraus committed
410 411
        const unsigned int size = bfuncs.size();

412
        AffineTransformation current_shift = shift;
kraus's avatar
kraus committed
413
        for (unsigned int i = 0; i < N_m; ++ i) {
kraus's avatar
kraus committed
414
            for (unsigned int j = 0; j < size; ++ j) {
kraus's avatar
kraus committed
415
                Base *obj = bfuncs[j]->clone();
416
                obj->trafo_m = obj->trafo_m.mult(current_shift);
kraus's avatar
kraus committed
417 418 419
                bfuncs.push_back(obj);
            }

420
            current_shift = current_shift.mult(shift);
kraus's avatar
kraus committed
421
        }
422 423 424
    }

    static
kraus's avatar
kraus committed
425 426 427
    bool parse_detail(iterator &it, const iterator &end, Function* &fun) {
        Repeat *rep = static_cast<Repeat*>(fun);
        if (!parse(it, end, rep->func_m)) return false;
428 429 430 431 432 433 434

        boost::regex argumentList("," + UInt + "," + Double + "," + Double + "\\)(.*)");
        boost::smatch what;

        std::string str(it, end);
        if (!boost::regex_match(str, what, argumentList)) return false;

kraus's avatar
kraus committed
435 436 437
        rep->N_m = atof(std::string(what[1]).c_str());
        rep->shiftx_m = atof(std::string(what[2]).c_str());
        rep->shifty_m = atof(std::string(what[4]).c_str());
438 439 440 441 442 443 444 445 446 447

        std::string fullMatch = what[0];
        std::string rest = what[6];

        it += (fullMatch.size() - rest.size());

        return true;
    }
};

kraus's avatar
kraus committed
448 449 450 451
struct Translate: public Function {
    Function* func_m;
    double shiftx_m;
    double shifty_m;
452

kraus's avatar
kraus committed
453 454
    virtual ~Translate() {
        delete func_m;
kraus's avatar
kraus committed
455 456
    }

457 458 459 460
    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
        std::cout << indent << "translate, " << std::endl;
kraus's avatar
kraus committed
461
        func_m->print(indentwidth + 8);
462
        std::cout << ",\n"
kraus's avatar
kraus committed
463 464
                  << indent2 << "dx: " << shiftx_m << ", \n"
                  << indent2 << "dy: " << shifty_m;
kraus's avatar
kraus committed
465 466
    }

kraus's avatar
kraus committed
467
    virtual void apply(std::vector<Base*> &bfuncs) {
468 469 470
        AffineTransformation shift(Vector_t(1.0, 0.0, -shiftx_m),
                                   Vector_t(0.0, 1.0, -shifty_m));

kraus's avatar
kraus committed
471
        func_m->apply(bfuncs);
kraus's avatar
kraus committed
472 473 474
        const unsigned int size = bfuncs.size();

        for (unsigned int j = 0; j < size; ++ j) {
kraus's avatar
kraus committed
475
            Base *obj = bfuncs[j];
476
            obj->trafo_m = obj->trafo_m.mult(shift);
kraus's avatar
kraus committed
477
        }
478 479 480
    }

    static
kraus's avatar
kraus committed
481 482 483
    bool parse_detail(iterator &it, const iterator &end, Function* &fun) {
        Translate *trans = static_cast<Translate*>(fun);
        if (!parse(it, end, trans->func_m)) return false;
484 485 486 487 488 489 490

        boost::regex argumentList("," + Double + "," + Double + "\\)(.*)");
        boost::smatch what;

        std::string str(it, end);
        if (!boost::regex_match(str, what, argumentList)) return false;

kraus's avatar
kraus committed
491 492
        trans->shiftx_m = atof(std::string(what[1]).c_str());
        trans->shifty_m = atof(std::string(what[3]).c_str());
493 494 495 496 497 498 499 500 501 502 503

        std::string fullMatch = what[0];
        std::string rest = what[5];

        it += (fullMatch.size() - rest.size());

        return true;
    }
};


kraus's avatar
kraus committed
504 505 506
struct Rotate: public Function {
    Function* func_m;
    double angle_m;
507

kraus's avatar
kraus committed
508 509
    virtual ~Rotate() {
        delete func_m;
kraus's avatar
kraus committed
510 511
    }

512 513 514 515
    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
        std::cout << indent << "rotate, " << std::endl;
kraus's avatar
kraus committed
516
        func_m->print(indentwidth + 8);
517
        std::cout << ",\n"
kraus's avatar
kraus committed
518
                  << indent2 << "angle: " << angle_m;
kraus's avatar
kraus committed
519 520
    }

kraus's avatar
kraus committed
521
    virtual void apply(std::vector<Base*> &bfuncs) {
522 523 524
        AffineTransformation rotation(Vector_t(cos(angle_m), sin(angle_m), 0.0),
                                      Vector_t(-sin(angle_m), cos(angle_m), 0.0));

kraus's avatar
kraus committed
525
        func_m->apply(bfuncs);
kraus's avatar
kraus committed
526 527 528
        const unsigned int size = bfuncs.size();

        for (unsigned int j = 0; j < size; ++ j) {
kraus's avatar
kraus committed
529
            Base *obj = bfuncs[j];
530
            obj->trafo_m = obj->trafo_m.mult(rotation);
kraus's avatar
kraus committed
531
        }
532 533 534
    }

    static
kraus's avatar
kraus committed
535 536 537
    bool parse_detail(iterator &it, const iterator &end, Function* &fun) {
        Rotate *rot = static_cast<Rotate*>(fun);
        if (!parse(it, end, rot->func_m)) return false;
538

kraus's avatar
kraus committed
539
        boost::regex argumentList("," + Double + "\\)(.*)");
540 541 542 543 544
        boost::smatch what;

        std::string str(it, end);
        if (!boost::regex_match(str, what, argumentList)) return false;

kraus's avatar
kraus committed
545
        rot->angle_m = atof(std::string(what[1]).c_str());
546 547 548 549 550 551 552 553 554 555

        std::string fullMatch = what[0];
        std::string rest = what[3];

        it += (fullMatch.size() - rest.size());

        return true;
    }
};

kraus's avatar
kraus committed
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
struct Shear: public Function {
    Function* func_m;
    double angleX_m;
    double angleY_m;

    virtual ~Shear() {
        delete func_m;
    }

    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
        std::cout << indent << "shear, " << std::endl;
        func_m->print(indentwidth + 8);
        if (std::abs(angleX_m) > 0.0) {
            std::cout << ",\n"
                      << indent2 << "angle X: " << angleX_m;
        } else {
            std::cout << ",\n"
                      << indent2 << "angle Y: " << angleY_m;
        }
    }

    virtual void apply(std::vector<Base*> &bfuncs) {
        AffineTransformation shear(Vector_t(1.0, tan(angleX_m), 0.0),
                                   Vector_t(-tan(angleY_m), 1.0, 0.0));

        func_m->apply(bfuncs);
        const unsigned int size = bfuncs.size();

        for (unsigned int j = 0; j < size; ++ j) {
            Base *obj = bfuncs[j];
            obj->trafo_m = obj->trafo_m.mult(shear);
        }
    }

    static
    bool parse_detail(iterator &it, const iterator &end, Function* &fun) {
        Shear *shr = static_cast<Shear*>(fun);
        if (!parse(it, end, shr->func_m)) return false;

        boost::regex argumentList("," + Double + "," + Double + "\\)(.*)");
        boost::smatch what;

        std::string str(it, end);
        if (!boost::regex_match(str, what, argumentList)) return false;

        shr->angleX_m = atof(std::string(what[1]).c_str());
        shr->angleY_m = atof(std::string(what[3]).c_str());

        if (std::abs(shr->angleX_m) > 0.0 && std::abs(shr->angleY_m) > 0.0)
            return false;

        std::string fullMatch = what[0];
        std::string rest = what[5];

        it += (fullMatch.size() - rest.size());

        return true;
    }
};

kraus's avatar
kraus committed
618 619
struct Union: public Function {
    std::vector<Function*> funcs_m;
620

kraus's avatar
kraus committed
621 622
    virtual ~Union () {
        for (Function* func: funcs_m) {
kraus's avatar
kraus committed
623 624 625 626
            delete func;
        }
    }

627 628 629 630 631 632
    virtual void print(int indentwidth) {
        std::string indent(indentwidth, ' ');
        std::string indent2(indentwidth + 8, ' ');
        std::string indent3(indentwidth + 16, ' ');
        std::cout << indent << "union, " << std::endl;
        std::cout << indent2 << "funcs: {\n";
kraus's avatar
kraus committed
633 634
        funcs_m.front()->print(indentwidth + 16);
        for (unsigned int i = 1; i < funcs_m.size(); ++ i) {
kraus's avatar
kraus committed
635 636
            std::cout << "\n"
                      << indent3 << "," << std::endl;
kraus's avatar
kraus committed
637
            funcs_m[i]->print(indentwidth + 16);
638 639
        }
        std::cout << "\n"
kraus's avatar
kraus committed
640 641 642
                  << indent2 << "} ";
    }

kraus's avatar
kraus committed
643 644 645 646
    virtual void apply(std::vector<Base*> &bfuncs) {
        for (unsigned int i = 0; i < funcs_m.size(); ++ i) {
            std::vector<Base*> children;
            Function *func = funcs_m[i];
647
            func->apply(children);
kraus's avatar
kraus committed
648
            bfuncs.insert(bfuncs.end(), children.begin(), children.end());
kraus's avatar
kraus committed
649
        }
650 651 652
    }

    static
kraus's avatar
kraus committed
653 654 655 656
    bool parse_detail(iterator &it, const iterator &end, Function* &fun) {
        Union *unin = static_cast<Union*>(fun);
        unin->funcs_m.push_back(NULL);
        if (!parse(it, end, unin->funcs_m.back())) return false;
657 658 659 660 661 662 663 664

        boost::regex argumentList("(,[a-z]+\\(.*)");
        boost::regex endParenthesis("\\)(.*)");
        boost::smatch what;

        std::string str(it, end);
        while (boost::regex_match(str, what, argumentList)) {
            iterator it2 = it + 1;
kraus's avatar
kraus committed
665
            unin->funcs_m.push_back(NULL);
666

kraus's avatar
kraus committed
667
            if (!parse(it2, end, unin->funcs_m.back())) return false;
668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684

            it = it2;
            str = std::string(it, end);
        }

        str = std::string(it, end);
        if (!boost::regex_match(str, what, endParenthesis)) return false;

        std::string fullMatch = what[0];
        std::string rest = what[1];

        it += (fullMatch.size() - rest.size());

        return true;
    }
};

kraus's avatar
kraus committed
685 686

bool parse(std::string str, Function* &fun) {
687
    str = boost::regex_replace(str, boost::regex("//.*?\\n"), std::string(""), boost::match_default | boost::format_all);
688 689 690 691 692 693 694 695 696 697 698
    str = boost::regex_replace(str, boost::regex("\\s"), std::string(""), boost::match_default | boost::format_all);
    iterator it = str.begin();
    iterator end = str.end();
    if (!parse(it, end, fun)) {
        std::cout << "parsing failed here:" << std::string(it, end) << std::endl;
        return false;
    }

    return true;
}

kraus's avatar
kraus committed
699
bool parse(iterator &it, const iterator &end, Function* &fun) {
700 701 702 703 704 705 706 707 708 709 710
    boost::regex functionCall(FCall);
    boost::smatch what;

    std::string str(it, end);
    if( !boost::regex_match(str , what, functionCall ) ) return false;

    std::string identifier = what[1];
    std::string arguments = what[2];
    unsigned int shift = identifier.size() + 1;

    if (identifier == "rectangle") {
kraus's avatar
kraus committed
711
        fun = new Rectangle;
kraus's avatar
kraus committed
712 713
        /*iterator it2 = */it += shift;
        if (!Rectangle::parse_detail(it, end, fun)) return false;
714

kraus's avatar
kraus committed
715
        // it = it2;
716 717
        return true;
    } else if (identifier == "ellipse") {
kraus's avatar
kraus committed
718
        fun = new Ellipse;
kraus's avatar
kraus committed
719 720
        /*iterator it2 = */it += shift;
        if (!Ellipse::parse_detail(it, end, fun)) return false;
721

kraus's avatar
kraus committed
722
        // it = it2;
723 724
        return true;
    } else if (identifier == "repeat") {
kraus's avatar
kraus committed
725
        fun = new Repeat;
kraus's avatar
kraus committed
726 727
        /*iterator it2 = */it += shift;
        if (!Repeat::parse_detail(it, end, fun)) return false;
728

kraus's avatar
kraus committed
729
        // it = it2;
730 731 732

        return true;
    } else if (identifier == "rotate") {
kraus's avatar
kraus committed
733
        fun = new Rotate;
kraus's avatar
kraus committed
734 735
        // iterator it2 =
            it += shift;
kraus's avatar
kraus committed
736
        if (!Rotate::parse_detail(it, end, fun)) return false;
737

kraus's avatar
kraus committed
738
        // // it = it2;
739 740 741

        return true;
    } else if (identifier == "translate") {
kraus's avatar
kraus committed
742
        fun = new Translate;
kraus's avatar
kraus committed
743 744
        /*iterator it2 = */it += shift;
        if (!Translate::parse_detail(it, end, fun)) return false;
745

kraus's avatar
kraus committed
746 747 748 749 750 751 752 753 754
        // it = it2;

        return true;
    } else if (identifier == "shear") {
        fun = new Shear;
        /*iterator it2 = */it += shift;
        if (!Shear::parse_detail(it, end, fun)) return false;

        // it = it2;
755 756 757

        return true;
    } else if (identifier == "union") {
kraus's avatar
kraus committed
758
        fun = new Union;
kraus's avatar
kraus committed
759 760
        /*iterator it2 = */it += shift;
        if (!Union::parse_detail(it, end, fun)) return false;
761

kraus's avatar
kraus committed
762
        // it = it2;
763 764 765 766 767 768 769 770 771

        return true;
    }


    return (it == end);

}

772
int main(int argc, char *argv[])
773
{
774 775 776 777
    if (argc < 2) {
        std::cout << "please provide the name of the file that contains your code" << std::endl;
        return 1;
    }
kraus's avatar
kraus committed
778
    Function *fun;
779 780 781 782 783 784

    std::ifstream t(argv[1]);
    std::string str((std::istreambuf_iterator<char>(t)),
                    std::istreambuf_iterator<char>());

    // std::string str("repeat( translate(union(rectangle(0.1, 0.1), ellipse(0.1, 0.1)), -0.01, -0.02), 2, 0.1, 0.2)");
785 786 787

    if (parse(str, fun)) {
        fun->print(0);
kraus's avatar
kraus committed
788 789
        std::cout << "\n" << std::endl;

kraus's avatar
kraus committed
790
        std::vector<Base*> baseBlocks;
kraus's avatar
kraus committed
791 792
        fun->apply(baseBlocks);

793
        std::ofstream out("test.gpl");
kraus's avatar
kraus committed
794
        for (Base* bfun: baseBlocks) {
kraus's avatar
kraus committed
795 796
            bfun->print(0);
            std::cout << std::endl;
797
            bfun->computeBoundingBox();
798
            bfun->writeGnuplot(out);
kraus's avatar
kraus committed
799
        }
800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819

        if (baseBlocks.size() > 0) {
            Vector_t llc, urc;
            Base* first = baseBlocks.front();
            const BoundingBox &bb = first->bb_m;
            llc = Vector_t(bb.center_m[0] - 0.5 * bb.width_m,
                           bb.center_m[1] - 0.5 * bb.height_m,
                           0.0);
            urc = Vector_t(bb.center_m[0] + 0.5 * bb.width_m,
                           bb.center_m[1] + 0.5 * bb.height_m,
                           0.0);

            for (unsigned int i = 1; i < baseBlocks.size(); ++ i) {
                const BoundingBox &bb = baseBlocks[i]->bb_m;
                llc[0] = std::min(llc[0], bb.center_m[0] - 0.5 * bb.width_m);
                llc[1] = std::min(llc[1], bb.center_m[1] - 0.5 * bb.height_m);
                urc[0] = std::max(urc[0], bb.center_m[0] + 0.5 * bb.width_m);
                urc[1] = std::max(urc[1], bb.center_m[1] + 0.5 * bb.height_m);
            }

820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868
            struct CoordinateComp {
                double delta_m;
                CoordinateComp(double delta):
                    delta_m(delta)
                { }

                bool operator()(double a, double b) const {
                    return a < b && std::abs(a - b) > delta_m;
                }
            };

            std::set<double, CoordinateComp> xCoordinates(CoordinateComp(1e-6 * (urc[0] - llc[0])));
            std::set<double, CoordinateComp> yCoordinates(CoordinateComp(1e-6 * (urc[1] - llc[1])));
            for (Base *func: baseBlocks) {
                Vector_t center = func->trafo_m.transformFrom(Vector_t(0.0));
                xCoordinates.insert(center[0]);
                yCoordinates.insert(center[1]);
            }
            std::set<double> dualXCoordinates;
            auto it = xCoordinates.begin();
            auto next = std::next(it);
            dualXCoordinates.insert(std::min(llc[0], 1.5 * (*it) - 0.5 * (*next)));
            for (; next != xCoordinates.end(); ++ it, ++ next) {
                dualXCoordinates.insert(0.5 * ((*next) + (*it)));
            }
            dualXCoordinates.insert(std::max(urc[0], 1.5 * (*it) - 0.5 * (*std::prev(it))));

            std::set<double> dualYCoordinates;
            it = yCoordinates.begin();
            next = std::next(it);
            dualYCoordinates.insert(std::min(llc[1], 1.5 * (*it) - 0.5 * (*next)));
            for (; next != yCoordinates.end(); ++ it, ++ next) {
                dualYCoordinates.insert(0.5 * ((*next) + (*it)));
            }
            dualXCoordinates.insert(std::max(urc[0], 1.5 * (*it) - 0.5 * (*std::prev(it))));

            for (auto it = dualXCoordinates.begin(); it != dualXCoordinates.end(); ++ it) {
                out << std::setw(14) << *it << std::setw(14) << *(dualYCoordinates.begin()) << "\n"
                    << std::setw(14) << *it << std::setw(14) << *(dualYCoordinates.rbegin()) << "\n"
                    << std::endl;
            }
            for (auto it = dualYCoordinates.begin(); it != dualYCoordinates.end(); ++ it) {
                out << std::setw(14) << *(dualXCoordinates.begin()) << std::setw(14) << *it << "\n"
                    << std::setw(14) << *(dualXCoordinates.rbegin()) << std::setw(14) << *it << "\n"
                    << std::endl;
            }
            out.close();

            out.open("particles.gpl");
869 870
            gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

871 872 873
            const unsigned int N = 100000;
            unsigned int n = 0;
            for (unsigned int i = 0; i < N; ++ i) {
874 875 876 877 878 879
                Vector_t X(0.0);
                X[0] = llc[0] + (urc[0] - llc[0]) * gsl_rng_uniform(rng);
                X[1] = llc[1] + (urc[1] - llc[1]) * gsl_rng_uniform(rng);

                for (Base* func: baseBlocks) {
                    if (func->isInside(X)) {
880
                        ++ n;
881 882 883 884 885 886 887
                        out << std::setw(14) << X[0]
                            << std::setw(14) << X[1]
                            << std::endl;
                        break;
                    }
                }
            }
888
            std::cout << (double)n / N * 100 << " % of particles passed" << std::endl;
889 890
            gsl_rng_free(rng);

kraus's avatar
kraus committed
891

892 893 894
            for (Base* func: baseBlocks) {
                delete func;
            }
kraus's avatar
kraus committed
895
        }
896 897

        out.close();
898
    }
kraus's avatar
kraus committed
899 900 901

    delete fun;

902 903
    return 0;
}