rickymexx / yocto-gl-vpt Goto Github PK

std::pair<float, vec3f> eval_spectral_tracking(vsdf& vsdf, float max_distance, rng_state& rng, const ray3f& ray) {
      auto object = vsdf.object;
      auto density_vol = object->density_vol;
      auto emission_vol = object->emission_vol;
      auto max_density = density_vol->max_voxel * object->density_mult;
      auto majorant_density = max_density;
      auto imax_density = 1.0f / max_density;
      auto path_length = 0.0f;
      auto tr = vec3f(1);
      auto f = vec3f(1);
      auto sigma_s = vsdf.scatter * max_density;
      auto Ps = sigma_s * imax_density;

      while (true) {
          path_length -= log(1 - rand1f(rng)) * imax_density;
          if (path_length >= max_distance)
              break;
          auto current_pos = ray.o + path_length * ray.d;
          auto sigma_t = vec3f(eval_vpt_density(vsdf, current_pos));
          auto sigma_n = majorant_density - sigma_t;
          auto sigma_a = sigma_t - sigma_s;
          auto Pt = sigma_t * imax_density;          
          auto Pn = sigma_n * imax_density;
          auto Pa = 1.0 - Ps - Pn;

          auto contrib = zero3f;

          tr *= exp(-sigma_t*path_length);  
 
          if (max(Pa) > rand1f(rng)) {    //ABSORPTION
              contrib = sigma_a / (majorant_density * Pa);
              contrib /= max(contrib);
              vsdf.event = EVENT_ABSORB;
              vsdf.density = sigma_t;
              return { path_length,  tr * contrib };
          }
          else if (max(Pt) > rand1f(rng)) {              // SCATTER
              contrib = sigma_s / (majorant_density * Ps);
              contrib /= max(contrib);
              vsdf.event = EVENT_SCATTER;
              vsdf.density = sigma_t;
              return { path_length, tr * contrib }; 
          }
          contrib = sigma_n;
          contrib /= max(contrib);
          f *= contrib;
      }
      return { max_distance, f };
  }

if (e == EVENT_ABSORB || e == EVENT_SCATTER) { 
    vsdf.event = e;
    vsdf.density = sigma_t;
    f *= tr;
    break;
}

std::pair<float, vec3f> eval_unidirectional_spectral_mis(vsdf& vsdf, float max_distance,
							   rng_state& rng, const ray3f& ray) {
    auto object           = vsdf.object;
    auto density_vol      = object->density_vol;
    auto emission_vol     = object->emission_vol;
    auto max_density      = density_vol->max_voxel * object->density_mult;
    auto imax_density     = 1.0f / max_density;
    auto path_length      = 0.0f;
    auto f                = vec3f(1);
    auto p                = vec3f(1);
    auto cc               = clamp((int)(rand1f(rng) * 3), 0, 2);
    auto current_pos      = ray.o; 
    auto majorant_density = max_density;
    auto tr               = vec3f(1);
    
    while (true) {
      path_length     -= majorant_density == 0.0f ? -flt_max :	log(1 - rand1f(rng)) / majorant_density;
      if (path_length >= max_distance)
	    break;
      current_pos = ray.o + path_length * ray.d;
      auto sigma_t     = vec3f(eval_vpt_density(vsdf, current_pos));
      auto sigma_s     = sigma_t * vsdf.scatter;
      auto sigma_a     = sigma_t - sigma_s;
      auto sigma_n     = vec3f(max_density) - sigma_t;
      vec3f sigma[3] = { sigma_n, sigma_s, sigma_a };
 
      tr *= exp(-sigma_t * path_length); 
      // Sample event 
      auto e = sample_event(sigma_a[cc] * imax_density,  sigma_s[cc] * imax_density, sigma_n[cc] * imax_density, rand1f(rng));

      if (e == EVENT_NULL)
          continue;
 
      f *= sigma[e];
      p  = f;
      
      // Populate vsdf with medium interaction information and return
      vsdf.event = e;
      vsdf.density = sigma_t;
      f *= tr;
      f = f / mean(p);       
      break;
    }    
    return {path_length, f };
  }

#include <yocto/yocto_cli.h>
#include <yocto/yocto_math.h>
#include <yocto/yocto_trace.h>
#include <yocto/yocto_image.h>

//openvdb includes
#include <openvdb/openvdb.h> 

#include <iostream>
  
using namespace yocto;
 
static inline bool save_yvol(
    const char* filename, int w, int h, int d, int nc, const float* voxels) {
    FILE* fp;
    fopen_s(&fp, filename, "wb");
    if (!fp) return false;
    auto fs_guard = std::unique_ptr<FILE, void (*)(FILE*)>{
        fp, [](FILE* f) { fclose(f); } };

    if (fprintf(fp, "YVOL\n") < 0) return false;
    if (fprintf(fp, "%d %d %d %d\n", w, h, d, nc) < 0) return false;
    auto nvalues = (size_t)w * (size_t)h * (size_t)d * (size_t)nc;
    if (fwrite(voxels, sizeof(float), nvalues, fp) != nvalues) return false;

    return true;
}

int main(int argc, char** argv)
{  

    if (argc != 2)
    {
        std::cout << "Usage: TestOpenVDB.exe [OpenVDB filename]" << std::endl; 
    }
    else
    {
        openvdb::initialize();
        openvdb::io::File file(argv[1]);
        file.open();
         
        openvdb::GridBase::Ptr baseGrid;

        for (openvdb::io::File::NameIterator nameIter = file.beginName(); nameIter != file.endName(); ++nameIter)
        {
            if (nameIter.gridName() == "density")
            {
                baseGrid = file.readGrid(nameIter.gridName());
            }
            else
            {
                std::cout << "Ignoring Grid: " << nameIter.gridName() << std::endl;
            }
        }

        openvdb::FloatGrid::Ptr grid = openvdb::gridPtrCast<openvdb::FloatGrid>(baseGrid);
        
        int min_x = 10000;
        int min_y = 10000;
        int min_z = 10000;

        int max_x = -10000;
        int max_y = -10000;
        int max_z = -10000;
         

        // Print all active ("on") voxels by means of an iterator.
        for (openvdb::FloatGrid::ValueOnCIter iter = grid->cbeginValueOn(); iter; ++iter)  
        {
            openvdb::Coord xyz = iter.getCoord(); 
            min_x = min(min_x, xyz.x());
            min_y = min(min_y, xyz.y());
            min_z = min(min_z, xyz.z());

            max_x = max(max_x, xyz.x());
            max_y = max(max_y, xyz.y());
            max_z = max(max_z, xyz.z());
 
        } 
        
        file.close();
        std::cout << "Input: " << argv[1] << std::endl;
         
        std::cout << "BB bounds: (" << min_x << ", " << min_y << ", " << min_z << ")" << std::endl;  
        std::cout << "BB bounds: (" << max_x << ", " << max_y << ", " << max_z << ")" << std::endl; 

        min_x = fabs(min_x);
        min_y = fabs(min_y);
        min_z = fabs(min_z);
        max_x = max_x + min_x + 1;
        max_y = max_y + min_y + 1;
        max_z = max_z + min_z + 1;
        std::cout << "BB bounds after shift: (" << max_x << ", " << max_y << ", " << max_z << ")" << std::endl;

        std::vector<float> odata(max_x * max_y * max_z, 0.0f);
        for (openvdb::FloatGrid::ValueOnCIter iter = grid->cbeginValueOn(); iter; ++iter)
        {
            openvdb::Coord xyz = iter.getCoord();
            odata[(xyz.z() + min_z) * max_y * max_x + (xyz.y() + min_y) * max_x + (xyz.x() + min_x)] = *iter;
        }
        
        std::string output = argv[1];
        output = output.substr(0, output.find_last_of(".")) + ".bin";
        save_yvol(output.c_str(), max_x, max_y, max_z, 1, odata.data());
         
        std::cout << "Volume file saved ["<< output << "]";
    }
    return 0;
}