[work 8] Delaunay diagram

Movie

Source code

about

・ドロネー図の作図過程を表示する
・ドロネー図についてはこちらのページを参考にさせていただきました。（注：このページの三角形の外接円の中心座標をそのまま使ったら上手く動作しなかったので、コードでは自分で修正したものを使用してます）

file

・上部にあるファイル名が表示されているボタンを押すと、表示されるファイルが切り替わります
・別ウィンドウ表示したい時や行番号などが無いRawMode表示したい時は、コード内右上のボタンを押してください（ボタンはマウスオーバーすると表示されます）

#include "ofMain.h"
#include "ofApp.h"

//======================================
int main( ){

    // 4K:4096x2160
    // 2K:2048x1080
    // FullHD:1920x1080
    // HD:1440x1080
    // HD720p:1280x720
    // DVD:720x480
	ofSetupOpenGL(1280,720, OF_WINDOW);

	// this kicks off the running of my app
	// can be OF_WINDOW or OF_FULLSCREEN
	// pass in width and height too:
	ofRunApp( new ofApp());

}

#pragma once

#include "ofMain.h"
#include <random>

#include "Delaunay.hpp"

#include "ImageRecorder.hpp"    // ImageRecorder


class ofApp : public ofBaseApp{
	public:
    ofApp();
    ~ofApp();
    
		void setup();
		void update();
		void draw();
		
		void keyPressed(int key);
		void keyReleased(int key);
		void mouseMoved(int x, int y);
		void mouseDragged(int x, int y, int button);
		void mousePressed(int x, int y, int button);
		void mouseReleased(int x, int y, int button);
		void mouseEntered(int x, int y);
		void mouseExited(int x, int y);
		void windowResized(int w, int h);
		void dragEvent(ofDragInfo dragInfo);
		void gotMessage(ofMessage msg);
    
private:
    Delaunay2d delaunay;
};

#include "ofApp.h"


ofApp::ofApp(){
    
}

ofApp::~ofApp(){
    
}

//--------------------------------------------------------------
void ofApp::setup(){
    double fps = 1;
    
    ofSetFrameRate(fps);
    ofBackground(255,255,255);
    ofSetBackgroundAuto(true);
    
    delaunay.setup();
}

//--------------------------------------------------------------
void ofApp::update(){
    delaunay.update();
}


//--------------------------------------------------------------
void ofApp::draw(){
    delaunay.display();
}

#ifndef Delaunay_hpp
#define Delaunay_hpp

#include <stdio.h>
#include "ofMain.h"


struct DVector {
    int x;
    int y;
    
    DVector operator+(const DVector& v) const
    {
        DVector tmp;
        tmp.x = x + v.x;
        tmp.y = y + v.y;
        return tmp;
    }
    
    DVector operator-(const DVector& v) const
    {
        DVector tmp;
        tmp.x = x - v.x;
        tmp.y = y - v.y;
        return tmp;
    }
    
    bool operator==(const DVector& v) const
    {
        return (x == v.x && y == v.y);
    }
    
    // for set/map
    bool operator<(const DVector& v) const
    {
        return x != v.x ? x < v.x : y < v.y;
    }
};


struct DCircle {
    DVector center;
    unsigned int radius;
};

class DTriangle {
public:
    DVector p1;
    DVector p2;
    DVector p3;
    DCircle circle;    // Circumscribed circle
    
    bool operator==(const DTriangle& t) const
    {
        return ((p1 == t.p1 && p2 == t.p2 && p3 == t.p3) ||
        (p1 == t.p2 && p2 == t.p3 && p3 == t.p1) ||
        (p1 == t.p3 && p2 == t.p1 && p3 == t.p2) ||
        
        (p1 == t.p3 && p2 == t.p2 && p3 == t.p1) ||
        (p1 == t.p2 && p2 == t.p1 && p3 == t.p3) ||
        (p1 == t.p1 && p2 == t.p3 && p3 == t.p2));
    }
    
    // for set/map
    bool operator<(const DTriangle& t) const
    {
        // Compare triangle areas
        DVector a1, a2;
        a1 = p2 - p1;
        a2 = p3 - p1;
        
        float s1 = abs(a1.x * a2.y + a1.y * a2.x) / 2;
        
        DVector a3, a4;
        a3 = t.p2 - t.p1;
        a4 = t.p3 - t.p1;
        
        float s2 = abs(a3.x * a4.y + a3.y * a4.x) / 2;
        
        return s1 < s2;
    }
};

typedef enum {
    D_MAKE_POINTS = 0,
    D_MAKE_HUGE_TRIANGLE,
    D_DIV_TRIANGLE,
    D_CLEAR_HUGE_TRIANGLE,
    D_STATE_END
} delaunayState;


class Delaunay2d {
public:
    Delaunay2d();
    ~Delaunay2d();
    void setup();
    void update();
    void display();
    
private:
    void getCircumscribedCircle(DVector p1, DVector p2, DVector p3, DCircle &c);
    void setDTriangle(DTriangle inT, std::map<DTriangle, bool> &outT);
    void delTriangleWithPoint(DVector p);
    void makeHugeTriangle(std::set<DVector> p, DTriangle &t);
    void divTriangles(DVector p, std::map<DTriangle, bool> inT, std::vector<DTriangle> &outT);
    
    const unsigned int max_points = 25;
    std::set<DVector> points;               // input points
    std::set<DVector>::iterator pitr;
    std::map<DTriangle, bool> triangles;    // output triangles
    delaunayState state;
    DTriangle hugeT;
    DVector currPoint;
};
#endif /* Delaunay_hpp */

#include "Delaunay.hpp"


Delaunay2d::Delaunay2d()
{
    
}

Delaunay2d::~Delaunay2d()
{
    
}

void Delaunay2d::setup()
{
    state = D_MAKE_POINTS;
}

void Delaunay2d::update()
{
    if (state == D_MAKE_POINTS) {
        for (int i = 0; i < max_points; i++) {
            DVector v;
            v.x = (int)ofRandom(400, 800);
            v.y = (int)ofRandom(100, 400);
            
            points.insert(v);
        }
        
        pitr = points.begin();
        
        state = D_MAKE_HUGE_TRIANGLE;
    } else if (state == D_MAKE_HUGE_TRIANGLE) {
        
        // To encircle the points with a triangle
        makeHugeTriangle(points, hugeT);
        
        // Set huge triangle to collection
        setDTriangle(hugeT, triangles);
        
        state = D_DIV_TRIANGLE;
    } else if (state == D_DIV_TRIANGLE) {
        DVector pt = *pitr;
        
        std::vector<DTriangle> tmpTriangles; // temporary
        
        // Divide the triangle into triangles
        divTriangles(pt, triangles, tmpTriangles);
        
        // Set triangles to collection
        if (tmpTriangles.empty() == false) {
            for (auto titr = tmpTriangles.begin(); titr != tmpTriangles.end(); titr++) {
                DTriangle t = *titr;
                setDTriangle(t, triangles);
            }
            
            tmpTriangles.clear();
        }
        
        currPoint = pt;
        
        if (pitr != points.end()) {
            pitr++;
        } else {
            state = D_CLEAR_HUGE_TRIANGLE;
        }
    } else if (state == D_CLEAR_HUGE_TRIANGLE) {
        delTriangleWithPoint(hugeT.p1);
        delTriangleWithPoint(hugeT.p2);
        delTriangleWithPoint(hugeT.p3);
        state = D_STATE_END;
    } else {
        // D_STATE_END
    }
}

void Delaunay2d::display()
{
    if (triangles.empty() == false) {
        for (auto titr = triangles.begin(); titr != triangles.end(); titr++) {
            if (titr->second != false) {
                ofSetColor(0, 0, 0);
                ofNoFill();
                ofDrawTriangle(titr->first.p1.x, titr->first.p1.y, titr->first.p2.x, titr->first.p2.y, titr->first.p3.x, titr->first.p3.y);
            }
        }
    }
    
    for (auto vitr = points.begin(); vitr != points.end(); vitr++) {
        ofSetColor(0, 255, 255);
        ofFill();
        ofDrawCircle(vitr->x, vitr->y, 5);
        
        auto p = points.find(currPoint);
        if (p != points.end()) {
            ofSetColor(255, 0, 0);
            ofDrawCircle(p->x, p->y, 6);
        }
    }
}


void Delaunay2d::getCircumscribedCircle(DVector p1, DVector p2, DVector p3, DCircle &c)
{
    // length(p1 - center) == length(p2 - center) == length(p3 - center)
    double k = 2.0 * ((p1.y - p2.y) * (p3.x - p2.x) - (p2.y - p3.y) * (p2.x - p1.x));
    double a1 = p1.x * p1.x + p1.y * p1.y;
    double a2 = p2.x * p2.x + p2.y * p2.y;
    double a3 = p3.x * p3.x + p3.y * p3.y;
    
    double x = (a1 * (p3.y - p2.y) + a2 * (p1.y - p3.y) + a3 * (p2.y - p1.y)) / (-1 * k);
    double y = (a1 * (p3.x - p2.x) + a2 * (p1.x - p3.x) + a3 * (p2.x - p1.x)) / k;

    double dx = x - p1.x;
    double dy = y - p1.y;
    double radius = sqrt(dx * dx + dy * dy);
    
    c.center.x = (int)(x + 0.5);    // roundup
    c.center.y = (int)(y + 0.5);    // roundup
    c.radius = (unsigned int)(radius + 0.5); // roundup
}

void Delaunay2d::setDTriangle(DTriangle inT, std::map<DTriangle, bool> &outT)
{
    auto itr = outT.find(inT);
    if (itr != outT.end()) {
        // This triangle is invalid
        outT.erase(itr);
        outT.insert(std::make_pair(inT, false));
    } else {
        // This triangle is valid
        outT.insert(std::make_pair(inT, true));
    }
}


void Delaunay2d::delTriangleWithPoint(DVector p)
{
    std::vector<DTriangle> tmpT;
    
    // Delete triangle include this point (p)
    for (auto titr = triangles.begin(); titr != triangles.end(); titr++) {
        if (p == titr->first.p1 || p == titr->first.p2 || p == titr->first.p3) {
            tmpT.push_back(titr->first);
        }
    }
    
    if (tmpT.empty() == false) {
        for (auto titr = tmpT.begin(); titr != tmpT.end(); titr++) {
            triangles.erase(triangles.find(*titr));
            triangles.insert(std::make_pair(*titr, false));
        }
    }
}


void Delaunay2d::makeHugeTriangle(std::set<DVector> p, DTriangle &t)
{
    int x_min = INT_MAX;
    int x_max = INT_MIN;
    int y_min = INT_MAX;
    int y_max = INT_MIN;
    
    // Rectangle with all points
    for (auto pitr = p.begin(); pitr != p.end(); pitr++) {
        if (pitr->x > x_max) {
            x_max = pitr->x;
        }
        if (pitr->x < x_min) {
            x_min = pitr->x;
        }
        
        if (pitr->y > y_max) {
            y_max = pitr->y;
        }
        if (pitr->y < y_min) {
            y_min = pitr->y;
        }
    }
    
    // A circle circumscribing the rectangle
    {
        int cx = x_min + ((x_max - x_min) / 2);
        int cy = y_min + ((y_max - y_min) / 2);
        int radius = (int)sqrt((cx - x_min) * (cx - x_min) + (cy - y_min) * (cy - y_min));
        
        DVector p1, p2, p3;
        p1.x = (int)(cx - sqrt(3.0) * radius);
        p1.y = (int)(cy - radius);
        p2.x = (int)(cx + sqrt(3.0) * radius);
        p2.y = p1.y;
        p3.x = cx;
        p3.y = (int)(cy + 2 * radius);
        
        DCircle c;
        getCircumscribedCircle(p1, p2, p3, c);
        
        t.p1 = p1;
        t.p2 = p2;
        t.p3 = p3;
        t.circle = c;
    }
}


void Delaunay2d::divTriangles(DVector p, std::map<DTriangle, bool> inT, std::vector<DTriangle> &outT)
{
    for (auto titr = inT.begin(); titr != inT.end(); titr++) {
        if (titr->second != false) {
            DTriangle t = titr->first;
            
            int dx = t.circle.center.x - p.x;
            int dy = t.circle.center.y - p.y;
            
            if ((unsigned long)((dx * dx) + (dy * dy)) < (unsigned long)(t.circle.radius * t.circle.radius)) {
                
                // Divide the triangle (t)
                
                DTriangle t1;
                DCircle c1;
                
                getCircumscribedCircle(p, t.p1, t.p2, c1);
                
                t1.p1 = p;
                t1.p2 = t.p1;
                t1.p3 = t.p2;
                t1.circle = c1;
                
                outT.push_back(t1);
                
                DTriangle t2;
                DCircle c2;
                
                getCircumscribedCircle(p, t.p2, t.p3, c2);
                
                t2.p1 = p;
                t2.p2 = t.p2;
                t2.p3 = t.p3;
                t2.circle = c2;
                
                outT.push_back(t2);
                
                
                DTriangle t3;
                DCircle c3;
                
                getCircumscribedCircle(p, t.p3, t.p1, c3);
                
                t3.p1 = p;
                t3.p2 = t.p3;
                t3.p3 = t.p1;
                t3.circle = c3;
                
                outT.push_back(t3);
                
                // This triangle (t) will be invalid at setDTriangle().
                outT.push_back(t);
            }
        }
    }
}