Quadtree ile Hızlandırılmış Geometrik Arama

Bir çok problemin çözümünde geometrik çözümlerden istifade edilebilir. Geometrik çözüm, bir meseleyi geometrik perspektif vasıtasıyla çözmektir, velev ki problemimiz doğrudan geometrik olarak tanımlanmış olmasın. Biraz daha açarsak, doğrudan geometrik olarak tanımlanmış problemlere örnek vermek gerekirse, kapalı bir sistemde hareket eden hava moleküllerinin bulunduğu kap yüzeyiyle ve birbirleriyle olan çarpışma testlerinin yapılması olabilir. Doğrudan geometrik olarak tanımlanmayan problemlere ise bir cemiyetteki insanların birbiriyle olan münasebetlerinin mevcut skaler veriler üzerinden geometrik yorumlama vasıtasıyla tespit edilmesi olabilir. İster doğrudan ister yorumlama vasıtasıyla olsun geometrik hale gelen verilerde en büyük problem hızlı bir şekilde arama yapabilmektir. Arama ile kastedilen işlem, belirli bir bölgedeki veyahut belirli bir noktanın belirli bir mesafedeki komşularının tespitidir. Bu işlemi gerçekleştirmek için akla ilk gelen çözüm en ilkel hali ile zorlama yöntemidir (brute-force). Bu yöntemde tahmin edebileceğiniz gibi aranan noktaların bulunması için tüm noktalar sırayla taranır ve referans noktamıza olan uzaklığına bakılır. Bu usul her zaman doğru sonuca kati suretle varacaktır, fakat muhtemelen en uzun yoldan, zira muhtemel olsun olmasın tüm noktaları taradık ve maliyeti azami tuttuk.

 

İnsan zihni bu problemi çözerken aslında çok da fazla zorlanmaz, zira sezgilerindeki gelişmişlik arama uzayında orta kısımda bulunan bir referans noktası için yapacağımız komşu aramasını gidipte çok uzaklarda aramaz, doğrudan ilgili kısma bakar. Fakat bilgisayar sistemlerinde genelde durum böyle değildir, veriler hafızada karışık bir sırayla tutulmaktadır (genelde üretim sırasına göre). Bu yüzden hangi nokta hangi noktaya ne kadar yakın bunu verinin geliş sırasında anlamak mümkün değildir. İşte hızlandırma tam da bu noktada mümkün olmaktadır. Aslında hızlandırma yöntemlerinde yapılan genel olarak düzensiz veriyi daha düzenli bir hale getirmekten ibarettir. Meseleyi daha iyi idrak etmek için doğrudan kullanacağımız veriyapısını detaylandıralım. Quadtree veri yapısı alsında bir çeşit ağaç yapısıdır. 2 boyutlu uzayda geçerli olmak üzere bir kök düğümden başlayarak her düğüm 4 alt düğüme ayrılmaktadır. Bu bölümleme kapsadığı noktalardan bağımsız olarak 4 eşit parçaya bölme şeklinde gerçekleşir. Nokta ihtiva eden düğümlere yaprak denmektedir. Bizim tasarımımızda bir yaprakta en fazla 1 nokta mevcut olabilir. Tarama işlemi ise, tıpkı ağacı oluşturduğumuz gibi, aradığımız nokta mevcut düğümün sol (alt-üst) ya da sağ (alt-üst) kısımlarından hangisine ait ise o bölgedeki alt ağacı tarayarak devam ediyoruz, böylelikle yakınlık ihtimali olmayan diğer bölgelerdeki noktaları devre dışı bırakmış ve bu yükten kurtulmuş oluyoruz, aynı insanın sezgisel olarak yaptığı gibi. Örnek ağaçlara bir göz atalım:

Veriyapısı ve algoritma hakkında daha detaylı bilgi için: http://en.wikipedia.org/wiki/Quadtree

Amacım algoritmayı tam anlamıyla anlatmak değil, sadece bir gerçekleme örneği sağlamak. Bu sebeple daha fazla detaylandırmak istemiyorum. Gerçekleme için C++ dili ve görselleştirme için OpenGL kütüphanesini kullandım. Performans arttırımı için dar boğaz olan dinamik bellek ayrımı (dynamic memory allocation) için bir iyileştirme gerçekleştirdim, fakat bu iyileştirmenin çalışabilmesi için maksimum nokta sınırının belirlenmesi gerekiyor, bu bilinmiyorsa daha farklı bir usul geliştirilebilir veya bu iyileştirme kaldırılabilir. Eğer kararlılık (stability) sizin için daha önemli ise bu iyileştirmeyi kaldırmanızı veya gözden geçirmenizi tavsiye ederim, zira bazı durumlarda (nokta yoğunlaşması gibi) hatalarla karşılaştım ayrıca aynı anda birden fazla ağaç üretmesine mani oluyor (no multiple instance).

Not: Nokta tanımında açık kaynak kodlu ve ücretsiz temin edilebilen GTL (https://code.google.com/p/geometry/) kütüphanesinden istifade ettim. 

/*
 * Point.h
 *
 *  Created on: Oct 24, 2013
 *      Author: ramazan
 */

#ifndef POINT_H_
#define POINT_H_

#include <gtl/vec2.hpp>
#include <gtl/vec3.hpp>

#include "Utility.h"
#include <algorithm>

using namespace gtl;

class Point {
public:
 Point() {
  color.setValue(0, 0, 1);
  location.setValue(randomFloat() * SCREEN_WIDTH, randomFloat() * SCREEN_HEIGHT);
 }

 void walk() {
  location += Vec2f(randomSignedFloat() * 30, randomSignedFloat() * 30);
  location[0] = std::max(location[0], 0.f);
  location[0] = std::min(location[0], (float)SCREEN_WIDTH);

  location[1] = std::max(location[1], 0.f);
  location[1] = std::min(location[1], (float)SCREEN_HEIGHT);
 }

 Vec2f location;
 Vec3f color;
};


#endif /* POINT_H_ */

/*
 * Utility.h
 *
 *  Created on: Oct 24, 2013
 *      Author: ramazan
 */

#ifndef UTILITY_H_
#define UTILITY_H_

#include <stdlib.h>
#include <stdio.h>

#define SCREEN_WIDTH 800
#define SCREEN_HEIGHT 600

#define THREESHOLD 0.000001
#define PARMEQ(x, y) (fabs(x - y) <= THREESHOLD)

static inline float randomFloat() {
 return (float)rand()/(float)RAND_MAX;
}

static inline float randomSignedFloat() {
 const float r = randomFloat();
 return r - 0.5f;
}


#endif /* UTILITY_H_ */

#ifndef QUADTREE_H_
#define QUADTREE_H_

#include <assert.h>
#include <vector>

#include "Point.h"

// You can remove all optimization lines together. 

#define MAX_POINT_CAPACITY 1000000 // Optimization line!

class QuadNode {
public:
 QuadNode() {
  leaf = true;
  point = NULL;
  children[0] = children[1] = children[2] = children[3] = NULL;
 }

 ~QuadNode() {
  for (int i = 0; i < 4; i++) {
   QuadNode *child = children[i];
   if (child != NULL)
    delete child;
  }
 }

 void insert(const Point &point) {
  if (leaf && this->point == NULL) {
   this->point = &point;
  }
  else {
   makeChild(-1, -1);
   makeChild(-1, 1);
   makeChild(1, -1);
   makeChild(1, 1);

   QuadNode *node = NULL;

   if (leaf) {
    const Point *thisPoint = this->point;
    this->point = NULL;
    node = findChild(*thisPoint);
    node->insert(*thisPoint);
   }

   node = findChild(point);
   node->insert(point);

   leaf = false;
  }
 }

 void query(const Vec2f &center, const float radius, std::vector<const Point *> &points) {
  const float absDX = fabs(center[0] - this->center[0]);
  const float absDY = fabs(center[1] - this->center[1]);

  if (absDX <= (halfSize[0] + radius) && absDY <= (halfSize[1] + radius)) {
   if (point != NULL) {
    Vec2f disn = point->location - center;
    const float dist = disn.length();

    if (dist <= radius)
     points.push_back(point);
   }
   else {
    for (int i = 0; i < 4; i++) {
     QuadNode *child = children[i];
     if (child != NULL)
      child->query(center, radius, points);
    }
   }
  }
 }

 void render() {
  const Vec2f p0 = center + Vec2f(-halfSize[0], halfSize[1]);
  const Vec2f p1 = center + Vec2f(-halfSize[0], -halfSize[1]);
  const Vec2f p2 = center + Vec2f(halfSize[0], -halfSize[1]);
  const Vec2f p3 = center + Vec2f(halfSize[0], halfSize[1]);

  glColor3f(1, 0, 0);
  glBegin(GL_LINES);
  glVertex2fv(p0.getValue());
  glVertex2fv(p1.getValue());

  glVertex2fv(p1.getValue());
  glVertex2fv(p2.getValue());

  glVertex2fv(p2.getValue());
  glVertex2fv(p3.getValue());

  glVertex2fv(p3.getValue());
  glVertex2fv(p0.getValue());
  glEnd();

  for (int i = 0; i < 4; i++) {
   QuadNode *child = children[i];
   if (child != NULL)
    child->render();
  }
 }
private:
 int makeChild(const float dX, const float dY) {
  int index = -1;
  float newX = halfSize[0] / 2;
  float newY = halfSize[1] / 2;
  if (dX < 0) {
   newX = -newX;
   if (dY < 0) {
    index = 1;
    newY = -newY;
   }
   else
    index = 0;
  } else {
   if (dY < 0) {
    index = 2;
    newY = -newY;
   }
   else
    index = 3;
  }

  verifyChild(index, center + Vec2f(newX, newY));

  return index;
 }
 QuadNode *findChild(const Point &point) {
  const float dX = point.location[0] - center[0];
  const float dY = point.location[1] - center[1];

  return children[makeChild(dX, dY)];
 }

 void verifyChild(const int index, const Vec2f newCenter) {
  assert(index >= 0 && index <= 3);

  if (children[index] == NULL) {
   children[index] = new QuadNode();
   assert(children[index] != NULL);
   children[index]->center = newCenter;
   children[index]->halfSize = halfSize / 2;
  }
 }

private:
 static char memory[]; // Optimization line!
 static int memoryIndex; // Optimization line!

public:
 void *operator new(size_t size) { // Optimization line!
  assert(memoryIndex < MAX_POINT_CAPACITY); // Optimization line!
  return &memory[memoryIndex++ * size]; // Optimization line!
 } // Optimization line!

 void operator delete(void *p) { // Optimization line!
  p = NULL; // Optimization line!
 } // Optimization line!

public:

 // 0 3
 // 1 2
 QuadNode *children[4];
 const Point *point;
 Vec2f center;
 Vec2f halfSize;
 bool leaf;
};

int QuadNode::memoryIndex = 0; // Optimization line!
char QuadNode::memory[sizeof(QuadNode) * MAX_POINT_CAPACITY] = { }; // Optimization line!

class Quadtree {
public:
 Quadtree(const Vec2f &center, const float w, const float h) {
  root = new QuadNode();
  root->center = center;
  root->halfSize.setValue(w / 2, h / 2);
 }

 ~Quadtree() {
  if (root != NULL) {
   delete root;
  }
 }

 void insert(const Point &point) {
  assert(root != NULL);

  root->insert(point);
 }

 void render() {
  assert(root != NULL);
  root->render();
 }

 void query(const Vec2f &center, const float radius, std::vector<const Point *> &points) {
  assert(root != NULL);

  root->query(center, radius, points);
 }

 QuadNode *root;
};


#endif /* QUADTREE_H_ */

// main.cpp

#include <iostream>
#include <glut\\glut.h>
#include <time.h>

#include "Quadtree.h"

#define RENDER_WAIT_TIME 25
#define NUMBER_OF_POINTS 30000

using namespace std;

Point *points = NULL;
Quadtree *qt;

Vec2f qCenter(300, 400);
float radius = 100;

void display(void) {
    glClear(GL_COLOR_BUFFER_BIT);
    glLoadIdentity();

    glBegin(GL_POINTS);
    for (int i = 0; i < NUMBER_OF_POINTS; i++) {
     glColor3fv(points[i].color.getValue());
     glVertex2fv(points[i].location.getValue());
    }
    glEnd();

    glColor3f(0, 1, 0);
    glBegin(GL_LINES);
    for (int i = 0; i < 360; i++) {
     Vec2f p1 = qCenter + Vec2f(cos((i * M_PI / 180.f)) * radius, sin((i * M_PI / 180.f)) * radius);
     Vec2f p2 = qCenter + Vec2f(cos(((i + 1) * M_PI / 180.f)) * radius, sin(((i + 1) * M_PI / 180.f)) * radius);
     glVertex2fv(p1.getValue());
     glVertex2fv(p2.getValue());
    }
    glEnd();

    if (qt != NULL)
     qt->render();

    //glFlush();
    glutSwapBuffers();
}

void mouse(int btn, int state, int x, int y) {
    //if(btn==GLUT_LEFT_BUTTON && state == GLUT_DOWN) axis = 0;
}

void search() {
 if (qt == NULL)
  return;
 std::vector<const Point *> qpoints;
    qt->query(qCenter, radius, qpoints);

    std::cout<<qpoints.size()<<std::endl;
    for (int i = 0; i < NUMBER_OF_POINTS; i++) {
  points[i].color.setValue(0, 0, 1);
     for (unsigned int j = 0; j < qpoints.size(); j++) {
      if (&points[i] == qpoints[j])
       points[i].color.setValue(0, 1, 0);
     }
    }
}

void mainLoop(int value) {
    glutPostRedisplay();
    glutTimerFunc(RENDER_WAIT_TIME, mainLoop, 0);
}

void build() {
 qt = new Quadtree(Vec2f(SCREEN_WIDTH / 2, SCREEN_HEIGHT / 2), SCREEN_WIDTH, SCREEN_HEIGHT);
    for (int i = 0; i < NUMBER_OF_POINTS; i++)
     qt->insert(points[i]);
}

void keyPressed(unsigned char key, int x, int y) {
    if (key == 'x')
        exit(0);

 if (key == 'w')
  qCenter += Vec2f(0, 10);

 if (key == 's')
  qCenter += Vec2f(0, -10);

 if (key == 'a')
  qCenter += Vec2f(-10, 0);

 if (key == 'd')
  qCenter += Vec2f(10, 0);

 if (key == 'r') {
  for (int i = 0; i < NUMBER_OF_POINTS; i++)
   points[i].walk();
  build();
 }

 search();
}

void myReshape(int w, int h) {
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(0.0f, 800, 0, 600, 1.0f, -1.0f);
    glMatrixMode(GL_MODELVIEW);
}

void keySpecial(int key, int x, int y) {
    if (key == GLUT_KEY_F4)
        exit(0);
}

void init() {
    glShadeModel(GL_FLAT);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glPointSize(3);

    points = new Point[NUMBER_OF_POINTS];
 //points[0].location.setValue(100, 100);
 //points[1].location.setValue(110, 120);
 //points[2].location.setValue(105, 105);
 build();
 search();
}

int main(int argc, char** argv) {
    srand(time(NULL));
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
    glutInitWindowSize(SCREEN_WIDTH, SCREEN_HEIGHT);
    glutCreateWindow("Quadtree");
    glutReshapeFunc(myReshape);
    glutDisplayFunc(display);
    glutKeyboardFunc(keyPressed);
    glutSpecialFunc(keySpecial);
    init();
    glutTimerFunc(RENDER_WAIT_TIME, mainLoop, 0);
    glutMouseFunc(mouse);
    glutMainLoop();
    return 0;
}

Bir kaç deneme yapalım:

10 adet nokta

 100 adet nokta

 5000 adet nokta

 30000 adet nokta
 

Bu bölümde 2 boyutlu düzlem üzerindeki geometrik problemlerin çözümünde kullanılan arama işlemini hızlandıran bir veriyapısını inceledik. Aynı veri yapısı 3 boyutlu uzayda da uygulanabilir (Octree - http://en.wikipedia.org/wiki/Octree). Mantık olarak aynı yalnızca 4 yerine 8 eşit parçaya ayrılmaktadır. Performans değerlendirmesini zaman cinsinden değerleri belirterek yapmak istemiyorum, kullandığım bilgisayar üzerinde brute-force yöntemine göre tahmin edebileceğiniz gibi (beklendiği üzere) çok büyük oranda hızlanma sağladı. Bunu kendi probleminize uygulayarak görebilirsiniz. Basit düzeyde anlatmaya çalıştım, umarım işinize yarar.