How to fill OpenVDB voxels that are inside a given plane?

Question

I have a quad defined by 4 (x,y,z) points (like a plane that has edges). I have an OpenVDB grid. I want to fill all the voxels with value 1 that are inside my quad (including edges). Is such thing possible with out setting each voxel of the quad (limited plane) manually?

Answer 1

There is an unoptimized method, for arbitrarily shaped planar quadrilaterals, you only need to input four vertices of the 3D space plane, and the output is filled planar voxels.

1. Get the vertex with the longest sum of adjacent sides as A, the adjacent vertices of A are B and D, and obtain the voxel coordinates and voxel numbers between A-B and A-D based on ray-cast, and in addition A vertex C-B and C-D progressively sample the same number of voxels.
2. Make one-to-one correspondence between the voxels adjacent to the two vertices A and D above, and fill the plane area based on ray-cast.

void VDBVolume::fillPlaneVoxel(const PlanarEquation& planar) {
  auto accessor = volume_->getUnsafeAccessor();
  const auto transform = volume_->transform();
  const openvdb::Vec3f vdbnormal(planar.normal.x(), planar.normal.y(), planar.normal.z());
  int longside_vtx[2];
  calVtxLongSide(planar.vtx, longside_vtx);
  // 1. ray cast long side
  std::vector<Eigen::Vector3f> longsidepoints;
  int neiborvtx[2];
  {
    const Eigen::Vector3f origin = planar.vtx[longside_vtx[0]];
    const openvdb::Vec3R eye(origin.x(), origin.y(), origin.z());
    
    GetRectNeiborVtx(longside_vtx[0], neiborvtx);
    for(int i = 0; i < 2; ++i) {
      Eigen::Vector3f direction = planar.vtx[neiborvtx[i]] - origin;
      openvdb::Vec3R dir(direction.x(), direction.y(), direction.z());
      dir.normalize();
      const float length = static_cast<float>(direction.norm());
      if(length > 50.f) {
        std::cout << "GetRectNeiborVtx length to large, something wrong: "<< i << "\n" << origin.transpose() << "\n"
                  << planar.vtx[neiborvtx[i]].transpose() << std::endl; 
        continue;
      }
      const float t0 =  -voxel_size_/2;
      const float t1 = length + voxel_size_/2;

      const auto ray = openvdb::math::Ray<float>(eye, dir, t0, t1).worldToIndex(*volume_);
      openvdb::math::DDA<decltype(ray)> dda(ray);
      do {
        const auto voxel = dda.voxel();
        const auto voxel_center_world = GetVoxelCenter(voxel, transform);
        longsidepoints.emplace_back(voxel_center_world);

      } while (dda.step());
    }
  }

  // 2. 在小边均匀采样相同数目点
  const int longsidepointnum = longsidepoints.size();
  std::vector<Eigen::Vector3f> shortsidepoints;
  shortsidepoints.resize(longsidepointnum);
  {
    // 输入：顶点、两邻接点，vtxs， 输出采样带年坐标，根据距离进行采样
    GenerateShortSidePoints(longside_vtx[1], neiborvtx, planar.vtx, shortsidepoints);
  }
  // 3. ray cast from longsidepoints to shortsidepoints
  // std::cout << "longsidepointnum: " << longsidepointnum << std::endl;
  for(int pid = 0; pid < longsidepointnum; ++pid) {
    const Eigen::Vector3f origin = longsidepoints[pid];
    const openvdb::Vec3R eye(origin.x(), origin.y(), origin.z());
    const Eigen::Vector3f direction = shortsidepoints[pid] - origin;
    openvdb::Vec3R dir(direction.x(), direction.y(), direction.z());
    dir.normalize();

    const float length = direction.norm();
    if(length > 50.f) {
      std::cout << "length to large, something wrong: "<< pid << "\n" << origin.transpose() << "\n"
                << shortsidepoints[pid].transpose() << std::endl; 
      continue;
    }
    const float t0 =  -voxel_size_/2;
    const float t1 = length + voxel_size_/2;
    const auto ray = openvdb::math::Ray<float>(eye, dir, t0, t1).worldToIndex(*volume_);
    openvdb::math::DDA<decltype(ray)> dda(ray);
    do {
      const auto voxel = dda.voxel();
      accessor.setValue(voxel, vdbnormal);
    } while (dda.step());
  }
}

Answer 2

If the four points build a rectangle, it could be possible using the

void fill(const CoordBBox& bbox, const ValueType& value, bool active = true);

function that exists in the Grid-class. It is not possible to transform the CoordBBox for rotations, instead you would have to do that by changing the transformation of the grid. With pseudo-code it could look like

CoordBBox plane; // created from your points
Transform old = grid.transform();
grid.setTransform(...); // Some transformation that places the grid correctly with respect to the plane
grid.fill(plane, 1);
grid.setTransform(old);

If this is not the case, you would have to set the values yourself.