src/warpcoord.h

/* This file is part of WarpTree
 *
 * Copyright (C) 2007 Jos van den Oever <jos@vandenoever.info>
 *               2001 Christophe Bouthier <christ.bouthier@free.fr>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public License
 * along with this library; see the file COPYING.LIB.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 * Boston, MA 02110-1301, USA.
 */

/*
   The functions
    void translate(const WarpCoord& t)
    void transform(const WarpCoord& c1, const WarpCoord& c2)
   contain code adapted from the program HyperTree by Christophe Bouthier and
   are used with permission.
*/

#ifndef WARPTREE_H
#define WARPTREE_H

#include <iostream>
#include <cmath>
#include <cassert>
#include <cstdlib>

/* This class provides the methods to handle the points in the graph.
 * It allows points to be translated in the unit square.
 */

#define ZONE_LENGTH 4 // size of the zone

// if WARPCOORD is not defined, the standard Poincare disk projection is used
#define WARPCOORD

class WarpCoord;
std::ostream& operator<<(std::ostream& o, const WarpCoord& e);
class WarpCoord {
private:
    double mx;
    double my;
    double md;

    /* the value for d is cached in this class. This function updates the
       value */
    inline void updateD() {
        md = sqrt(d2());
    }
    /* multiply this point with point z */
    void
    operator*=(const WarpCoord& z) {
        double tx = mx;
        mx = tx*z.mx - my * z.my;
        my = tx*z.my + my * z.mx;
        updateD();
    }
    /* divide this point by point z */
    void
    operator/=(const WarpCoord& z) {
        double d2 = z.d2();
        double tx = mx;
        mx = (tx*z.mx + my*z.my) / d2;
        my = (my*z.mx - tx*z.my) / d2;
        updateD();
    }
public:
    WarpCoord() :mx(0), my(0), md(0) {}
    WarpCoord(double x, double y) :mx(x), my(y) {
        updateD();
    }
    inline double x() const { return mx; }
    inline double y() const { return my; }
    inline double d() const { return md; }
    inline double d2() const { return mx*mx + my*my; }
    /* calculate the arg of this point */
    double arg() const {
        double a = atan(my / mx);
        if (mx < 0) {
            a += M_PI;
        } else if (my < 0) {
            a += 2 * M_PI;
        }
        return a;
    }
    /* translate this point by adding the vector t */
    void translate(const WarpCoord& t) {
        double denX = mx * t.mx + my * t.my + 1;
        double denY = my * t.mx - mx * t.my;
        double dd = denX * denX + denY * denY;

        double numX = mx + t.mx;
        double numY = my + t.my;

        mx = (numX * denX + numY * denY) / dd;
        my = (numY * denX - numX * denY) / dd;
        updateD();
    }
    /* zoom this point. this means the point move closer to or away from the
       center of the square */
    void zoom(double f) {
        assert(f > 0);
        if (md < 0.001) return;
        double r = md;
        double c = pow((1+md)/(1-md),f);
        c = (c-1)/(c+1)/r;
        mx *= c;
        my *= c;
        md *= c;
    }
    /* move this point from the disk positiong to the square position */
    void stretch() {
        if (md < 0.0001) return;
#ifdef WARPCOORD
        double fx = fabs(mx);
        double fy = fabs(my);
        if (fy > 0 && fy > fx) {
            double f = md/fy;
            my = (my>0) ?md :-md;
            mx *= f;
            md *= f;
        } else if (fx > 0) {
            double f = md/fx;
            mx = (mx>0) ?md :-md;
            my *= f;
            md *= f;
        }
#endif
    }
    /* move this point from the square position to the disk position */
    void unstretch() {
        assert(md >= 0);
        if (md < 0.0001) return;
#ifdef WARPCOORD
        double f = 1;
        double fx = fabs(mx);
        double fy = fabs(my);
        if (fy > 0 && fy > fx) {
            f = md/fy;
            md = (my>0) ?my :-my;
        } else if (fx > 0) {
            f = md/fx;
            md = (mx>0) ?mx :-mx;
        }
        mx /= f;
        my /= f;
#endif
    }
    /* is this point visible, in other words, is it not too close to the edge
       of the unit square */
    inline bool visible() const {
        return mx > -0.95 && mx < 0.95 && my > -0.95 && my < 0.95;
    }
    /* is this point valid, in other words, does it lie inside the unit
       square */
    inline bool valid() const {
        return mx > -1 && mx < 1 && my > -1 && my < 1;
    }
    /* return the negative value of this point */
    inline WarpCoord operator-() const {
        return WarpCoord(-mx, -my);
    }
    /* substract the value of point e from this point. Note that this is not
       equal to translation. */
    inline void operator-=(const WarpCoord& e) {
        mx -= e.mx;
        my -= e.my;
        updateD();
    }
    /* apply the translation that is needed to get from point c1 to c2 to this
       point */
    void transform(const WarpCoord& c1, const WarpCoord& c2) {
        WarpCoord p(c1.mx + c2.mx, c1.my + c2.my);

        WarpCoord d(c2.mx, -c2.my);
        d *= c1;
        d.mx += 1;
        p /= d;

        WarpCoord o(c1.mx, -c1.my);
        o *= c2;
        o.mx += 1;
        o /= d;

        WarpCoord z(*this);
        *this *= o;
        mx += p.mx;
        my += p.my;

        d.set(p.mx, -p.my);
        d *= z;
        d *= o;
        d.mx += 1;

        *this /= d;
        updateD();
    }
    /* change the value of this point */
    inline void set(double x, double y) {
        mx = x;
        my = y;
        updateD();
    }
};

#endif 

---