NurbsPerformer/src/device/Nurbs/bvh.cu

//
// Created by 14727 on 2022/12/11.
//

#include "../../../include/device/Nurbs/bvh.cuh"
/**
 * 当前层的第一个节点的坐标
 * @param layer 层号。根节点为第一层
 */
__device__ int getStartIdxOfLayerN(int layer) {
    if (layer == 1) return 0;
    int num = 1;
    // 找规律得出
    for (int i = 2; i < layer; i++) {
        num = num * 4 + 1;
    }
    return num;
}

/**
 * 根据u、v参数域中的小矩形的位置，判断它在BVH叶节点层（也就是最后一层）中的位置
 * @param ix
 * @param iy
 * @return
 */
__device__ int getChildNodeIdx(int ix, int iy) {
    int idx = 0;
    while(ix > 0) {
        int log2XCeil = int(__log2f(ix));
        idx += int(pow(4, log2XCeil));
        ix -= int(pow(2, log2XCeil));
    }
    while(iy > 0) {
        int log2XCeil = int(__log2f(iy));
        idx += 2 * int(pow(4, log2XCeil));
        iy -= int(pow(2, log2XCeil));
    }
    return idx;
}


__global__ void
buildSurfaceBvh(const float *k, int level, const float *evaluatedPoints, float lastKnot_u, float lastKnot_v,
                int sampleCnt_u, int sampleCnt_v, BVHNode *BVH) {
    // 假设还是二维grid和二维的block
    int ix = blockIdx.x * blockDim.x + threadIdx.x;
    int iy = blockIdx.y * blockDim.y + threadIdx.y;
    if (ix >= sampleCnt_u - 1 || iy >= sampleCnt_v - 1) {
        return;
    }
    float singleStepVal_u = lastKnot_u / (sampleCnt_u - 1);
    float singleStepVal_v = lastKnot_v / (sampleCnt_v - 1);

    float u = singleStepVal_u * ix, v = singleStepVal_v * iy;

    int bias = getChildNodeIdx(ix, iy);  // 每层的节点在总数组中的编号相对该层第一个节点的偏置
    for (int step = 1; step <= sampleCnt_u - 1 || step <= sampleCnt_v - 1; step = step * 2, --level) {

        float step_u = min(step, sampleCnt_u - 1);
        float step_v = min(step, sampleCnt_v - 1);

        int baseIdx = getStartIdxOfLayerN(level);  // 当前层的第一个节点的坐标
        int idx = baseIdx + bias;

        int tmpIdx = (ix * sampleCnt_v + iy) * 3;
//        AABB bound(FVec3(evaluatedPoints[tmpIdx] + *k, evaluatedPoints[tmpIdx + 1] + *k,
//                         evaluatedPoints[tmpIdx + 2] + *k));  // 最初只包含u，v点
        AABB bound(FVec3(evaluatedPoints[tmpIdx], evaluatedPoints[tmpIdx + 1],
                         evaluatedPoints[tmpIdx + 2]));  // 最初只包含u，v点
        if(idx == 75) printf("75: %g, %g, %g\n", evaluatedPoints[tmpIdx], evaluatedPoints[tmpIdx + 1], evaluatedPoints[tmpIdx + 2]);
        if (step == 1) {
            // 叶子节点
            BVH[idx].u0 = u;
            BVH[idx].u1 = u + step_u * singleStepVal_u;
            BVH[idx].v0 = v;
            BVH[idx].v1 = v + step_v * singleStepVal_v;
            BVH[idx].level = level;
//            int tmpIdx1 = ((ix + 1) * sampleCnt_v + iy) * 3;
//            bound = bound.Union(FVec3(evaluatedPoints[tmpIdx1] + *k, evaluatedPoints[tmpIdx1 + 1] + *k,
//                                      evaluatedPoints[tmpIdx1 + 2] + *k));
//            int tmpIdx2 = (ix * sampleCnt_v + iy + 1) * 3;
//            bound = bound.Union(FVec3(evaluatedPoints[tmpIdx2] + *k, evaluatedPoints[tmpIdx2 + 1] + *k,
//                                      evaluatedPoints[tmpIdx2 + 2] + *k));
//            int tmpIdx3 = ((ix + 1) * sampleCnt_v + iy + 1) * 3;
//            bound = bound.Union(FVec3(evaluatedPoints[tmpIdx3] + *k, evaluatedPoints[tmpIdx3 + 1] + *k,
//                                      evaluatedPoints[tmpIdx3 + 2] + *k));
            int tmpIdx1 = ((ix + 1) * sampleCnt_v + iy) * 3;
            bound = bound.Union(FVec3(evaluatedPoints[tmpIdx1], evaluatedPoints[tmpIdx1 + 1],
                                      evaluatedPoints[tmpIdx1 + 2]));
            int tmpIdx2 = (ix * sampleCnt_v + iy + 1) * 3;
            bound = bound.Union(FVec3(evaluatedPoints[tmpIdx2], evaluatedPoints[tmpIdx2 + 1],
                                      evaluatedPoints[tmpIdx2 + 2]));
            int tmpIdx3 = ((ix + 1) * sampleCnt_v + iy + 1) * 3;
            bound = bound.Union(FVec3(evaluatedPoints[tmpIdx3], evaluatedPoints[tmpIdx3 + 1],
                                      evaluatedPoints[tmpIdx3 + 2]));

            BVH[idx].bounds = bound;
            BVH[idx].childNum = -1;


        } else {
            if (ix % step == 0 || iy % step == 0) {
                // 非叶子节点
                BVH[idx].u0 = u;
                BVH[idx].u1 = u + step_u * singleStepVal_u;
                BVH[idx].v0 = v;
                BVH[idx].v1 = v + step_v * singleStepVal_v;
                BVH[idx].level = level;

                int firstChild = 4 * idx + 1;
                bound = bound.Union(BVH[firstChild].bounds);
                bound = bound.Union(BVH[firstChild + 1].bounds);
                bound = bound.Union(BVH[firstChild + 2].bounds);
                bound = bound.Union(BVH[firstChild + 3].bounds);
                BVH[idx].firstChild = firstChild;
                BVH[idx].bounds = bound;
            }
        }

        bias /= 4;

        __syncthreads();
    }
}

__host__ void BVH::printQuadTree() {

    for (int l = 1, levelNodesCnt = 1, baseIdx = 0; l <= maxLevel; l++, levelNodesCnt *= 4, baseIdx = baseIdx * 4 + 1) {
        for (int biasIdx = 0; biasIdx < levelNodesCnt; biasIdx++) {
            int idx = baseIdx + biasIdx;
            auto pMax = nodes[idx].bounds.pMax;
            auto pMin = nodes[idx].bounds.pMin;
            printf("<<%g, %g, %g>,<%g, %g, %g>>  ", pMin.x, pMin.y, pMin.z, pMax.x, pMax.y, pMax.z);
        }
        printf("\n =============$$$$$$$============= \n");
    }
}