四分木

Pythonで四分木を実装してみた。ランダムな座標を持ったを四分木に挿入したとき、下のような分割になった。プロット数は10000個。

まぁ全体的に均等に分散していて、深さもそれほど変わってないことが見て分かるね。大数の法則だよね。具体的に言うと、深さは7〜14になった。(トップレベルの深さを1としてる)
サンプル数が1000のときは5〜12、100では4〜7になった。スケールしてるっぽい。


こちらはx,yそれぞれについてガウス分布にしたがわせたもの。具体的に言うと、平均0.5,分散0.1と指定した。

真ん中に集中しているのは当たり前なのだけど、いい感じにデータが取れてるね。こちらは深さは3〜17と、まあまあ深くなっているもののまだまだ平和な値。
サンプル数を1000に減らしたときは深さ3〜14,100では3〜8という値になった。はじっこの方は全然データが来ないので、サンプル数が増えてもあまり深さのミニマムは変わらない。案外と、ランダムの場合と比較しても深さはそれほど変わらない。


プロット数100000のときは、ランダムで8〜18、ガウス分布で4〜19。そろそろ精度限界が疑われる。ぱいそんだし。


以下ソース。きたない。あとdeleteをテストしてないのでたぶんバグがある。【あとで直す】

class Region ():
    def __init__ (self):
        pass
    def includes(x, y):
        pass
    def rectangles(x,y):
        pass



class RoundRegion (Region):
    def __init__ (self, x, y, r):
        self.x = x
        self.y = y
        self.r = r


    def includes(self, x, y):
        xDiff = self.x - x
        yDiff = self.y - y
        if xDiff * xDiff + yDiff * yDiff <= self.r * self.r:
            return True
        return False


    def rectangles(self,  qa):
        if self._rect(qa.minX,qa.minY) or self._rect(qa.minX, qa.maxY) or self._rect(qa.maxX, qa.minY) or self._rect(qa.maxX, qa.maxY):
            return True
        else:
            return False


    def _rect(self, x, y):
        if x >= self.x - self.r and x < self.x + self.r :
            if y >= self.y - self.r and y < self.y + self.r :
                return True
        return False



class RectRegion (Region):
    def __init__ (self, area):
        self.minX, self.maxX = area.minX, area.maxX
        self.minY, self.maxY = area.minY, area.maxY


    def includes (self, x, y):
        if x >= self.minX and x < self.maxX:
            if y >= self.minX and y < self.maxY:
                return True
        return False


    def _rectangles (self, x, y):
        return self.includes(x,y)



class QElements ():
    def __init__ (self, x, y, obj):
        self.x = x
        self.y = y
        self.obj = obj


    def equals (self, elem):
        thr = 0.0
        xDiff = self.x - elem.x
        yDiff = self.y - elem.y
        if xDiff <= thr and yDiff <= thr:
            return True
        else:
            return False



class QArea ():
    def __init__ (self, minX, maxX, minY, maxY):
        self.minX, self.maxX = minX, maxX
        self.midX            = (minX + maxX) / 2.0
        self.minY, self.maxY = minY, maxY
        self.midY            = (minY + maxY) / 2.0


    def contains (self, elem):
        if elem.x >= self.minX and elem.x < self.maxX:
            if elem.y >= self.minY and elem.y < self.maxY:
                return true
        return false



class QNode (QArea):
    def __init__ (self, elem, qa):
        self.leaf            = elem
        QArea.__init__ (self, qa.minX, qa.maxX, qa.minY, qa.maxY)
        self.region = {'NE':None,'NW':None,'SW':None,'SE':None}
#    def ___init__ (self, elem, minX, maxX, minY, maxY):
#        self.leaf            = elem
#        QArea.__init__ (self, minX, maxX, minY, maxY)
#        self.region = {'NE':None,'NW':None,'SW':None,'SE':None}


    def kids (self):
        res = []
        for key, val in self.region.iteritems:
            res.append(key)
        return res


    def getRegion (self, reg):
        if reg[0] == 'N':
            minY, maxY = self.midY, self.maxY
        elif reg[0] == 'S':
            minY, maxY = self.minY, self.midY
        else:
            return None

        if reg[1] == 'W':
            minX, maxX = self.minX, self.midX
        elif reg[1] == 'E':
            minX, maxX = self.midX, self.maxX
        else:
            return None

        return QArea(minX, maxX, minY, maxY)

    def getArea(self, x, y):
#        print self.minY, self.midY, self.maxY, y

        if y >= self.minY and y < self.midY:
            reg = 'S'
        elif y >= self.midY and y < self.maxY:
            reg = 'N'
        else:
            return None

#        print self.minX, self.midX, self.maxX, x
        if x >= self.minX and x < self.midX:
            reg += 'W'
        elif x >= self.midX and x < self.maxX:
            reg += 'E'
        else:
            return None

#        print reg
        qa =self.getRegion(reg)

        return reg, qa


    def insert (self, elem):
        if self.leaf:
#            print self.leaf.x + self.leaf.y
#            print "      leaf:", self.minX, self.midX, self.maxX, elem.x
#            print "           ", self.minY, self.midY, self.maxY, elem.y
            (reg, qa) = self.getArea(self.leaf.x, self.leaf.y)
            self.region[reg] = QNode(self.leaf, qa)
            self.leaf = None

        else:
#            print "      node:", self.minX, self.midX, self.maxX, elem.x
#            print "           ", self.minY, self.midY, self.maxY, elem.y
            pass

        (reg, qa) = self.getArea(elem.x, elem.y)
        if self.region[reg] is None:
            self.region[reg] = QNode(elem, qa)
        else:
            self.region[reg].insert(elem)


    def delete (self, elem):
        if self.leaf:
            if self.leaf.equals(elem):
                return self.leaf, True, True
            else:
                return None, False, None

        (reg, qa) = self.getArea(elem.x, elem.y)
        if self.region[reg] is None:
            return None, False, None

        (res, desc, leaf) = self.delete (elem)
        if not desc:
            return res, desc, leaf

        elif self.region[reg].leaf:
            # next to leaf node; the beginning
            kids = self.kids()
            if kids == 1:
                return res, True, None
            if kids == 2:
                for k in self.region:
                    if reg != k:
                        break
                self.region[reg] = None
                if self.region[k].leaf:
                    return res, True, self.region[k].leaf
                else:
                    return res, False, None
            else:
                r, qa = self.getArea(leaf.x, leaf.y)
                self.region[req] = QNode(leaf, qa)
                return res, False, None

        else:
            # passby node
            kids = self.kids()
            if kids == 1:
                return res, desc, leaf
            if kids == 2 and not leaf:
                for k in self.region:
                    if reg != k:
                        break
                self.region[reg] = None
                if self.region[k].leaf:
                    return res, True, self.region[k].leaf
                else:
                    return res, False, None
            else:
                r, qa = self.getArea(leaf.x, leaf.y)
                self.region[req] = QNode(leaf, qa)
                return res, False, None


    def searchItem(self, elem):
        if self.leaf:
            if self.leaf.equals(elem):
                return self.leaf
            else:
                return None
        (reg, qa) = self.getArea(elem.x, elem.y)
        if self.region[reg] is None:
            return None, False, None
        else:
            return self.search(elem)


    def searchRegion(self, reg):
        if self.leaf:
            if reg.includes (self.leaf.x, self.leaf.y):
                return [self.leaf]
            else:
                return []

        lst = []
        for k in self.region:
            if not self.region[k]:
                continue
            qa = self.getRegion(k)
            if reg.rectangles (qa):
                l = self.region[k].searchRegion(reg)
                if l:
                    lst.extend(l)
        return lst


    def dump(self):
        if self.leaf:
            print "leaf", self.leaf.x, self.leaf.y
            return

        for k in self.region:
            if self.region[k]:
                print self.region[k].minX, self.region[k].maxX, self.region[k].minY, self.region[k].maxY
                self.region[k].dump()


    def getDepth(self, min, max, depth):
        if self.leaf:
            if min > depth:
                min = depth
            if max < depth:
                max = depth
            return (min, max)

        else:
            for k in self.region:
                if not self.region[k]:
                    continue
                min_x, max_x = self.region[k].getDepth(min, max, depth+1)
                if min_x < min:
                    min = min_x
                if max_x > max:
                    max = max_x
            return (min, max)
                


class QTree (QNode):
    def __init__ (self):
        self.tree = None


    def insert (self, x, y, obj):
        elem = QElements(x, y, obj)
#        print "insert:", x, y
        if self.tree :
            self.tree.insert(elem)
        else:
            self.tree = QNode (elem, QArea(0,1,0,1))


    def delete (self, elem):
        if not self.tree:
            return None

        res, susp, leaf = self.tree.delete(elem)
        if susp:
            self.tree = QNode (leaf, QArea(0,1,0,1))
        return res


    def get (self, elem):
        if not self.tree:
            return None

        return self.tree.searchItem(elem)


    def search (self, x, y, r):
        if not self.tree:
            return None

        region = RoundRegion(x, y, r)
        return self.tree.searchRegion(region)


    def dump (self):
        if not self.tree:
            return None
        
        return self.tree.dump()


    def getDepth (self):
        if not self.tree:
            return None

        return self.tree.getDepth (999999, 0, 1)

def okay(x,y):
    if x < 0 or x >= 1:
        return False
    if y < 0 or y >= 1:
        return False
    return True


if __name__ == '__main__':
    import random

    tree = QTree()

    for x in xrange(10000):
        x, y = -1,-1
        while not okay(x,y):
            x, y = random.random(), random.random()
#            x, y = random.gauss(0.5, 0.1), random.gauss(0.5, 0.1)
        tree.insert(x, y, None)

    min,max = tree.getDepth()
    print "Depth:",min,max
    tree.dump()

gdで描くソース。てきとう。最終出力が白黒なのに線は赤だとか、そゆところを直さず公開してしまうのはひどい。

import gd

sz = 4000
im = gd.image((sz, sz))

for l in open("gauss~"):
    lst = l.split(" ")
    if lst[0] == "leaf":
#        x,y = float(lst[1]), float(lst[2])
#        x,y = int(x*sz), int(x*sz)
#        im.setPixel((x,y), im.colorAllocate((0,255,0)))
        pass
    elif lst[0] == "Depth:":
        pass
    else:
        sx, bx = int(sz * float(lst[0])), int(sz * float(lst[1]))
        sy, by = int(sz * float(lst[2])), int(sz * float(lst[3]))
        im.lines(((sx,sy),(sx,by),(bx,by),(bx,sy),(sx,sy)), im.colorAllocate((255,0,0)))

f = open("gauss.wbmp", "w")
im.writeWbmp(f)
f.close()