PlantSimulationLab / Helios

The Helios simulation system is a versatile modeling framework that handles tasks such as managing geometry and associated data structures through a C++ API. Plug-ins build off of the Helios core engine, and access the geometry and data via the Helios context. The sytem comes with a visualization plug-in that can produce stunning renderings of model geometry and data with relatively little effort.
http://baileylab.ucdavis.edu/software/helios
GNU General Public License v2.0
61 stars 27 forks source link

ply file is not correctly loaded and calculation is not correctly executed #40

Open mashiro210 opened 4 months ago

mashiro210 commented 4 months ago

I ran Helios code (wrote based on tutorial 10 with slight modification) in for loop to calculate light interception rate for 120 ply files. But ply file is not correctly loaded and the calculation result was quite different between result without using for loop. How can I deal with such a problem?

bnbailey-psl commented 4 months ago

I am happy to take a look. Perhaps you can upload your main code or a snippet of the part for loading the PLY files? Based on your description I am unsure what you mean by loading the 120 files with and without the loop.

mashiro210 commented 4 months ago

Thank you for your reply. This is my code for run Helios with for loop. File path and output file name are contained in config.csv

// Author: Mashiro
// Last update: 2024/6/7
// Description: run helios for lighting simulation with ply file; for loop

// load modules
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <vector>
#include <chrono>
#include <ctime>
#include <iomanip>
#include "Visualizer.h"
#include "RadiationModel.h"

// use name space
using namespace helios; // using the helios namespace so you can omit 'helios::' before names

// define data structure for having each row of csv file

struct PlotInfo {
    std::string flightName;
    std::string filePath;
    float zMinValue;
    std::string outputFileName;

};

// define function for reading csv file

std::vector<PlotInfo> readCSV(const std::string& filename) {
    std::vector<PlotInfo> allPlotInfo;
    std::ifstream file(filename);

    if (!file.is_open()) {
        std::cerr << "Could not open the file!" << std::endl;
        return allPlotInfo;
    }

    std::string line;
    // skip header
    std::getline(file, line);

    while (std::getline(file, line)) {
        std::stringstream ss(line);
        std::string item;
        PlotInfo plotInfo;

        std::getline(ss, plotInfo.flightName, ',');
        std::getline(ss, plotInfo.filePath, ',');
        std::getline(ss, plotInfo.outputFileName, ',');
        std::getline(ss, item, ',');
        plotInfo.zMinValue = std::stof(item);
        std::getline(ss, item, ',');

        allPlotInfo.push_back(plotInfo);
    }

    file.close();
    return allPlotInfo;
}

// define function for get current date

std::string getCurrentDate8Digit() {
    // get current time
    auto now = std::chrono::system_clock::now();
    // convert to time_t
    std::time_t now_c = std::chrono::system_clock::to_time_t(now);

    // get current local time
    std::tm local_time = *std::localtime(&now_c);

    // conver to 8digit strings
    std::stringstream ss;
    ss << std::put_time(&local_time, "%Y%m%d");

    return ss.str();
}

// execution part

int main() {

    // get current date
    std::string today = getCurrentDate8Digit();

    // *** 1. set parameters *** //

    // parameters for lighting simulation

    // const float tanTheta = 0.15 / 0.69;
    const float azimuthAngle = std::atan(0.15 / 0.69);
    const float elevationAngle = deg2rad(90);
    const float pfd = 2000; // photon flux density (micro mol / m^2)
    const float irradiance = pfd / 4.6; // convert pfd to irradiance

    // path setting
    const std::string eastSensingAreaPlyPath = "/mnt/b2e33138-44a3-48b8-acb8-d5bbf2a4da69/Tanashi/2023_Soybean/Data/UAV/pointCloud/soybean_tanashi_quantam_light_sensor_sensing_area_mesh_10cm_east.ply";
    const std::string westSensingAreaPlyPath = "/mnt/b2e33138-44a3-48b8-acb8-d5bbf2a4da69/Tanashi/2023_Soybean/Data/UAV/pointCloud/soybean_tanashi_quantam_light_sensor_sensing_area_mesh_10cm_west.ply";
    const std::string outputPathBase = "/mnt/b2e33138-44a3-48b8-acb8-d5bbf2a4da69/Tanashi/2023_Soybean/Result" + std::string("/") + today;
    const std::string resultCsvPath = "/mnt/b2e33138-44a3-48b8-acb8-d5bbf2a4da69/Tanashi/2023_Soybean/Result" + std::string("/") + today + std::string("/") + today + "_light_interception_simulation_result.csv";

    // generate result file and set header
    std::ofstream outputFile(resultCsvPath);
    // outputFile << "outputFileName,eastLightInterceptionRate,westLightInterceptionRate\n";
    outputFile << "outputFileName,eastPAR,westPAR\n";

    Context context;

    // parameters for visualization of simulation result
    const float radiusForCameraPosition = 6;
    const float cmin = 0; // for colorbar
    const float cmax = 450; // for colorbar
    const uint fontSize = 20; // font size of color bar

    // *** 2. load csv file for loop *** //
    const std::string csvPath = "/mnt/b2e33138-44a3-48b8-acb8-d5bbf2a4da69/Tanashi/2023_Soybean/Result/20240710/Helios_config.csv";
    const std::vector<PlotInfo> allPlotInfo = readCSV(csvPath);

    // *** 3. execute by for loop *** //

    for (const auto& plotInfo : allPlotInfo) {

        // *** a. load ply files with a base position of (0, 0, 0) and no scaling *** //
        std::string targetPlyPath = plotInfo.filePath;
        helios::vec3 basePos = make_vec3(0, 0, plotInfo.zMinValue);
        std::vector<uint> UUIDs_plant = context.loadPLY(targetPlyPath.c_str(), basePos, 0, nullrotation, RGB::black, "ZUP");
        std::vector<uint> UUIDs_eastSensingArea = context.loadPLY(eastSensingAreaPlyPath.c_str(), basePos, 0, nullrotation, RGB::black, "ZUP");
        std::vector<uint> UUIDs_westSensingArea = context.loadPLY(westSensingAreaPlyPath.c_str(), basePos, 0, nullrotation, RGB::black, "ZUP");

        // *** b. Radiation model set-up *** //

        RadiationModel radiation(&context); // declare and initialize radiation model
        // add sun radiation source with angle (elevation & azimuth) info
        uint sourceID = radiation.addSunSphereRadiationSource(make_SphericalCoord(elevationAngle, azimuthAngle));

        radiation.addRadiationBand("PAR");
        radiation.disableEmission("PAR"); // turn off emission, no emission of primitives in solar bands
        radiation.setSourceFlux(sourceID, "PAR", irradiance); // set solar flux perpendicular to sun direction

        context.setPrimitiveData(UUIDs_eastSensingArea, "twosided_flag", uint(1)); // assume ground only absorb radiation from the top
        context.setPrimitiveData(UUIDs_westSensingArea, "twosided_flag", uint(1)); // assume ground only absorb radiation from the top
        radiation.updateGeometry(); // update geometry and radiation setting

        // *** c. run the model and calculate PAR interception *** //
        radiation.runBand("PAR");

        // Calculate PAR interception
        // east
        float PAR_eastSensingArea; // declare variable
        context.calculatePrimitiveDataAreaWeightedMean(UUIDs_eastSensingArea, "radiation_flux_PAR", PAR_eastSensingArea);
        // float eastInterceptionRate = 1 - (PAR_eastSensingArea / irradiance);
        // west
        float PAR_westSensingArea; // declare variable
        context.calculatePrimitiveDataAreaWeightedMean(UUIDs_westSensingArea, "radiation_flux_PAR", PAR_westSensingArea);
        // float westInterceptionRate = 1 - (PAR_westSensingArea / irradiance);

        // write result to outputfile; result.csv
        // outputFile << plotInfo.outputFileName << "," << eastInterceptionRate << "," << westInterceptionRate << "\n";
        outputFile << plotInfo.outputFileName << "," << PAR_eastSensingArea << "," << PAR_westSensingArea << "\n";

        // *** d. visualize the results *** //

        std::string outputPath = outputPathBase + std::string("/") + plotInfo.flightName + std::string("/") + plotInfo.outputFileName;

        Visualizer visualizer(1200); // lanch ploty window

        // define camera position
        visualizer.setCameraPosition(SphericalCoord(radiusForCameraPosition, elevationAngle, (azimuthAngle + deg2rad(180))), basePos);
        // draw result
        visualizer.buildContextGeometry(&context);
        visualizer.colorContextPrimitivesByData("radiation_flux_PAR"); // color primitives based on a pseudocolor mapping of primitive data "radiation_flux_PAR" (this is the output primitive data from the radiation model)

        // change color bar setting
        visualizer.setColorbarTitle("PAR flux [W/m^2]" ); // colorbar a title
        visualizer.setColorbarPosition(make_vec3(0.25, 0.25, 0.5)); // color bar position
        visualizer.setColorbarRange(cmin, cmax);
        visualizer.setColorbarFontSize(fontSize);

        visualizer.plotUpdate();
        visualizer.printWindow(outputPath.c_str());
        visualizer.closeWindow();

    }

    outputFile.close();

    return 0;
}
mashiro210 commented 4 months ago

result without for loop mean hard coding case for loading ply file

bnbailey-psl commented 4 months ago

I think the issue is that you are loading all of the PLY models and placing them on top of each other, which is probably not what you want. When you call Context::loadPLY() with a vec3 position argument, it will automatically translate the model so that its origin is at the position you specify. My guess is that you just want to load the model 'as-is' without applying any scaling, translation, or rotation? For that, you can just call the overloaded version with only a file path argument: std::vector<uint> UUIDs_plant = context.loadPLY(targetPlyPath.c_str());

You just need to be sure that your PLY files were exported with the z-axis as the 'up' direction in this case, since you don't have that argument to specify it anymore.

mashiro210 commented 4 months ago

Thank you for your advice. I will try it.