배열 쿼드트리

#include <iostream>
#include <vector>
#include <queue>

const int NULL_INDEX = -1;

constexpr int MAX_NODES = 1000;
constexpr int MAX_PARTICLES = 1000;

float nodeX[MAX_NODES];
float nodeY[MAX_NODES];
float nodeWidth[MAX_NODES];
float nodeHeight[MAX_NODES];
int nodeParticleIndex[MAX_NODES];
int nodeChildren[MAX_NODES][4];

float particleX[MAX_PARTICLES];
float particleY[MAX_PARTICLES];
float particleMass[MAX_PARTICLES];
int nodeCount = 0;


unsigned int expandBits(unsigned int v)
{
    v = (v | (v << 8)) & 0x00FF00FF;
    v = (v | (v << 4)) & 0x0F0F0F0F;
    v = (v | (v << 2)) & 0x33333333;
    v = (v | (v << 1)) & 0x55555555;

    return v;
}

unsigned int morton2D(float x, float y)
{
    x += 32768.0f;
    y += 32768.0f;

    auto scaledX = (unsigned int) x;
    auto scaledY = (unsigned int) y;

    scaledX = std::max(0U, std::min(scaledX, 65535U));
    scaledY = std::max(0U, std::min(scaledY, 65535U));

    unsigned int xx = expandBits(scaledX);
    unsigned int yy = expandBits(scaledY);

    return xx | (yy << 1);
}

int createNode(float x, float y, float width, float height)
{
    int index = nodeCount++;
    nodeX[index] = x;
    nodeY[index] = y;
    nodeWidth[index] = width;
    nodeHeight[index] = height;
    nodeParticleIndex[index] = NULL_INDEX;
    for (int i = 0; i < 4; i++)
    {
        nodeChildren[index][i] = NULL_INDEX;
    }
    return index;
}

void insertParticleToNode(int nodeIndex, int particleIndex)
{
    if (nodeParticleIndex[nodeIndex] == NULL_INDEX && nodeChildren[nodeIndex][0] == NULL_INDEX)
    {
        nodeParticleIndex[nodeIndex] = particleIndex;
        return;
    }

    if (nodeParticleIndex[nodeIndex] != NULL_INDEX)
    {
        int existingParticleIndex = nodeParticleIndex[nodeIndex];
        nodeParticleIndex[nodeIndex] = NULL_INDEX;

        float halfWidth = nodeWidth[nodeIndex] / 2.0f;
        float halfHeight = nodeHeight[nodeIndex] / 2.0f;

        int childIndex = 0;
        if (particleX[existingParticleIndex] > nodeX[nodeIndex] + halfWidth)
        {
            childIndex |= 1;
        }
        if (particleY[existingParticleIndex] > nodeY[nodeIndex] + halfHeight)
        {
            childIndex |= 2;
        }

        if (nodeChildren[nodeIndex][childIndex] == NULL_INDEX)
        {
            float childX = nodeX[nodeIndex] + (childIndex & 1 ? halfWidth : 0);
            float childY = nodeY[nodeIndex] + (childIndex & 2 ? halfHeight : 0);
            nodeChildren[nodeIndex][childIndex] = createNode(childX, childY, halfWidth, halfHeight);
        }
        insertParticleToNode(nodeIndex, existingParticleIndex);
    }

    float halfWidth = nodeWidth[nodeIndex] / 2.0f;
    float halfHeight = nodeHeight[nodeIndex] / 2.0f;

    int childIndex = 0;
    if (particleX[particleIndex] > nodeX[nodeIndex] + halfWidth)
    {
        childIndex |= 1;
    }
    if (particleY[particleIndex] > nodeY[nodeIndex] + halfHeight)
    {
        childIndex |= 2;
    }

    if (nodeChildren[nodeIndex][childIndex] == NULL_INDEX)
    {
        float childX = nodeX[nodeIndex] + (childIndex & 1 ? halfWidth : 0);
        float childY = nodeY[nodeIndex] + (childIndex & 2 ? halfHeight : 0);

        nodeChildren[nodeIndex][childIndex] = createNode(childX, childY, halfWidth, halfHeight);
    }

    insertParticleToNode(nodeChildren[nodeIndex][childIndex], particleIndex);
}

void buildTree()
{
    nodeCount = 0;
    int rootNodeIndex = createNode(-32768.0f, -32768.0f, 65536.0f, 65536.0f);


    for (int i = 0; i < MAX_PARTICLES; i++)
    {
        if (particleMass[i] > 0)
        {
            insertParticleToNode(rootNodeIndex, i);
        }
    }
}

void printNode(int nodeIndex, int depth = 0)
{
    for (int i = 0; i < depth; i++)
    {
        std::cout << "  ";
    }

    std::cout << "Node [" << nodeX[nodeIndex] << ", " << nodeY[nodeIndex] << ", "
              << nodeWidth[nodeIndex] << ", " << nodeHeight[nodeIndex] << "]";

    if (nodeParticleIndex[nodeIndex] != NULL_INDEX)
    {
        std::cout << " Particle: (" << particleX[nodeParticleIndex[nodeIndex]] << ", "
                  << particleY[nodeParticleIndex[nodeIndex]] << ")";
    }

    std::cout << std::endl;

    for (int i = 0; i < 4; i++)
    {
        if (nodeChildren[nodeIndex][i] != NULL_INDEX)
        {
            printNode(nodeChildren[nodeIndex][i], depth + 1);
        }
    }
}

int main()
{
    particleX[0] = -10000.0f;
    particleY[0] = -10000.0f;
    particleMass[0] = 1.0f;

    particleX[1] = 10000.0f;
    particleY[1] = 10000.0f;
    particleMass[1] = 1.0f;

    particleX[2] = -10000.0f;
    particleY[2] = 10000.0f;
    particleMass[2] = 1.0f;

    particleX[3] = 10000.0f;
    particleY[3] = -10000.0f;
    particleMass[3] = 1.0f;


    buildTree();

    printNode(0);

    return 0;
}

Barnes-Hut n체 시뮬레이션 만들 때 기존의 포인터 트리로는 캐시 접근이 비효율적이라서 배열로 캐시 히트를 최대한으로 높인 트리.

 

파티클 배열 인덱스를 morton 코드 순서로 정렬해주면 공간 상에서 생성되는 노드의 위치를 한 곳으로 모을 수 있어 캐시 히트를 더욱 향상시킬 수 있음.

 

나중에는 힐베르트 커브도 구현해보고 싶지만, 일단 GPU와 CPU에서 빠르게 동작하는 코드를 만들고 OpenGL 태양계 시뮬레이션(CPU)와 은하 시뮬레이션(GPU)를 만드는 게 우선이므로 넘어감.

 

기존에 만들어본 Barnes-Hut n체는 다른 논문들에 비하면 속도가 좀 많이 느려서 한 10만체 정도는 만들고 싶기도 하고 게임에서 토성의 고리 등을 실시간 시뮬레이션 할 때 꼭 필요한지라 만들어봄.

 

하고 싶은 프로젝트는 많은데 다 할 수 있을지 모르겠다.

 

하고싶은 거 목록:

이징 모델, 인공신경망 예측

Smooth Particle Hydrodynamics

나비에 스톡스 시뮬레이션