cli_pipeline.cpp

This program demonstrates how to use most features of the pipeline, as well as allowing use of the pipeline itself through execution.

/****************************************************************************************
 * SCYLLARUS : C++ HYPERSPECTRAL PROCESSING LIBRARY
 * cli_pipeline - Hyperspectral Image Pipeline usage example
 *
 * Example program showing how to use the pipeline object.
 *
 * This computer code is subject to copyright:
 * (c) National ICT Australia Limited (NICTA) 2013-2014 All Rights Reserved.
 * 
 * Jeremy Oorloff, National ICT Australia (NICTA)
 *
 ***************************************************************************************/

// External Includes
#include <armadillo>
#include <boost/program_options.hpp>
#include <boost/filesystem.hpp>
#include <opencv2/opencv.hpp>
#include <memory>
#include <chrono>

// Internal Includes
#include "scyllarus.h"

// Boost Program Options namespace alias
namespace po = boost::program_options;

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

    // Variables to store input options in
    bool        show_images;
    bool        recover_illuminant;
    bool        recover_dichromatic;
    bool        recover_materials;
    bool        process;
    bool        save = false;
    bool        keep_time;
    bool        save_refl;
    bool        nurbs;
    bool        filter;
    bool        quick_mat;
    bool        calc_error;
    int         debug;
    float       ill_alpha;
    int         ill_patch_size;
    int         dic_neigh_size;
    int         dic_gray_thres;
    int         mat_max_clust;
    int         mat_abund_num;
    float       resample;
    float       mat_temp_max;
    float       mat_temp_min;
    float       mat_cool_rate;
    float       mat_split_th;
    std::string filename;
    std::string savename;

    // Declare input options and link them to variables
    po::options_description desc("Allowed options");
    desc.add_options()
        ("help,h",           "Convey usage information")
        ("file,f",           po::value<std::string>(&filename)->default_value(""),   "Path to input file")
        ("save,o",           po::value<std::string>(&savename)->default_value(""),   "Path to output file HSZ or HDR (.hdr, .fla, .raw etc)")
        ("illuminant,i",     po::bool_switch(&recover_illuminant),                   "Perform illuminant recovery")
        ("dichromatic,d",    po::bool_switch(&recover_dichromatic),                  "Perform dichromatic decompose")
        ("material,m",       po::bool_switch(&recover_materials),                    "Perform material recovery")
        ("process,p",        po::bool_switch(&process),                              "Perform whole processing pipeline ( -i -d -m -a)")
        ("show_images,s",    po::bool_switch(&show_images),                          "Show images")
        ("print_timing,t",   po::bool_switch(&keep_time),                            "Print timing information (Walltime)")
        ("no_nurbs,n",       po::bool_switch(&nurbs),                                "Use NURBS when encoding HSZ")
        ("filter,w",         po::bool_switch(&filter),                               "Filter image before processing")
        ("quick_mat,q",      po::bool_switch(&quick_mat),                            "Use Fast material clustering method")
        ("reflectance",      po::bool_switch(&save_refl),                            "Save Reflectance Cube")
        ("calc_error,e",      po::bool_switch(&calc_error),                          "Calculate and print the RMS Error")
        ("resample,r",       po::value<float>(&resample)->default_value(1.0f),       "Resample the image by a factor ( 1 = None, 2 = /2, 4 = /4 etc.)")
        ("verbosity,v",      po::value<int>(&debug)->default_value(2),               "Quantity of output desired (0 = none, 5 = max)")
        ("ill_alpha",        po::value<float>(&ill_alpha)->default_value(50.0f),     "Illuminant Recovery Alpha value ( > 0, default = 50)")
        ("ill_patch_size",   po::value<int>(&ill_patch_size)->default_value(20),     "Illuminant patch size ( > 1, default = 20)")
        ("dic_size",         po::value<int>(&dic_neigh_size)->default_value(5),      "Dichromatic Neighbourhood Size ( > 0, default 5)")
        ("dic_gray_thres",   po::value<int>(&dic_gray_thres)->default_value(2),      "Dichromatic Gray Threshold ( > 0, default 2)")
        ("mat_max_clust",    po::value<int>(&mat_max_clust)->default_value(20),      "Material Recovery Max Clusters ( > 0, default 20)")
        ("mat_abund_num",    po::value<int>(&mat_abund_num)->default_value(5),       "Material Abundance per pixel ( > 1, default 5)")
        ("mat_temp_max",     po::value<float>(&mat_temp_max)->default_value(0.02f),  "Material Recovery Max Temperature ( > Min Temp, default 0.02)")
        ("mat_temp_min",     po::value<float>(&mat_temp_min)->default_value(0.00025f),"Material Recovery Min Temperature ( < Max Temp, default 0.00025)")
        ("mat_cool_rate",    po::value<float>(&mat_cool_rate)->default_value(0.8f),  "Material Recovery Cooling Rate ( <0, default 0.8)")
        ("mat_split_th",     po::value<float>(&mat_split_th)->default_value(-1),     "Material Recovery Cluster Split Threshold ( < 0, default cos(5pi / 180))")
        ;
    // Parse options
    po::variables_map vm;
    po::store(po::parse_command_line(argc, argv, desc), vm);
    po::notify(vm);

    // Print Version information
    //scyl::print_version();

    // Print help if requested
    if (vm.count("help")) {
        std::cout << desc << std::endl;
        return 0;
    }

    // Declare Input Object shared pointer
    std::shared_ptr<scyl::input>  input;
    std::shared_ptr<scyl::output> output;

    // Determine the file extensions of the specified files for input and output
    std::string file_extension  = boost::filesystem::extension(filename);
    std::string save_extension  = boost::filesystem::extension(savename);

    // Depending on the specified file type, create an input or print an error.
    if (file_extension.compare(".hsz") == 0 ||
        file_extension.compare(".HSZ") == 0)
    {
        // Make a HSZ input object. The act of creating this object loads the HSZ file.
        input = std::make_shared<scyl::hsz_input>(filename);
        std::cout << "Loaded " << filename << std::endl;

    } else if (!file_extension.empty())
    {
        // Make a HDR input object. Note if you want to load a .hdr style file you can specify either .hdr, (or any datafile ext associated with a .hdr)
        // The act of creating this object loads the HDR file.
        input = std::make_shared<scyl::hdr_input>(filename);
        std::cout << "Loaded " << filename << std::endl;

    } else
    {
        std::cout << "The file you have supplied for input is not supported. Try .hsz or .hdr files." << std::endl;
        std::cout << desc << std::endl;
        return -1;
    }

    // If an output filename was specified, make an output object.
    if (savename != "")
    {
        // Determine which sort of object to make
        if (save_extension.compare(".hsz") == 0 ||
            save_extension.compare(".HSZ") == 0)
        {
            // Make a HSZ output object.
            output = std::make_shared<scyl::hsz_output>();
            save = true;
        } else if (!save_extension.empty())
        {
            // Any other file extension, make a HDR output object (if the ext is not .hdr, it is assumed to be the data file extension)
            output = std::make_shared<scyl::hdr_output>(save_refl);
            save = true;
        } else
        {
            // No extension was specified so we have an error.
            std::cout << "The file you have supplied for out is not supported. Try .hsz or .hdr files." << std::endl;
            std::cout << desc << std::endl;
            return -1;
        }
    }

    // Make a Pipeline object, pass it the input(|output) object to use
    scyl::pipeline pipeline;
    // Check which constructor to use
    if (save)
    {
        // Pipeline with both input and output
        pipeline = scyl::pipeline(input, output);
    } else
    {
        // Pipeline with just input
        pipeline = scyl::pipeline(input);
    }

    // Configure settings for the pipeline
    pipeline.set_debug_level(debug);
  //pipeline.set_illuminant_method(hs::ILLUMINANT_METHOD in);
    pipeline.set_illuminant_recovery_alpha(ill_alpha);
    pipeline.set_illuminant_patch_size(ill_patch_size);
  //pipeline.set_dichromatic_method(hs::DICHROMATIC_METHOD in);
    pipeline.set_dichromatic_neighbourhood_size(dic_neigh_size);
    pipeline.set_dichromatic_gray_threshold(dic_gray_thres);
    if (quick_mat)
    {
        pipeline.set_material_method(scyl::MATERIAL_DAQ);
    } else {
        pipeline.set_material_method(scyl::MATERIAL_DA);
    }
    pipeline.set_material_max_clusters(mat_max_clust);
    pipeline.set_material_temperature_max(mat_temp_max);
    pipeline.set_material_temperature_min(mat_temp_min);
    pipeline.set_material_cooling_rate(mat_cool_rate);
    pipeline.set_material_split_threshold(mat_split_th);
    pipeline.set_abundance_num_materials(mat_abund_num);
    pipeline.set_nurbs_save_as(nurbs);


    // Filter (if desired)
    if (filter)
    {
        std::cout << "Filtering Image..." << std::endl;
        pipeline.filter(scyl::FILTER_WIENER, 5);
    }

    // Resample (if desired), float safe comparison.
    if (!(resample < 1.0f + arma::Datum<float>::eps and resample > 1.0f - arma::Datum<float>::eps) and resample > 0.0f)
    {
        pipeline.resize(pipeline.height()/resample, pipeline.width()/resample);
    }

    // Warning for large images
    if (pipeline.height() > 1500 || pipeline.width() > 1500)
    {
        std::cout << "Image is very large, may take long time and a lot of memory to process" << std::endl << "Use '-r' to resample image smaller" << std::endl;
    }

    // Perform processing. if -p was selected run the whole pipeline (by setting all the other flags to true)
    // Note we could call the pipeline.process() method for this but it is less desirable in this instance.
    if (process)
    {
        recover_illuminant  = true;
        recover_dichromatic = true;
        recover_materials   = true;
    }

    // Timing Code
    std::chrono::time_point<std::chrono::system_clock> begin, end;
    std::chrono::duration<float> time_ill = std::chrono::duration<float>(0);
    std::chrono::duration<float> time_dic = std::chrono::duration<float>(0);
    std::chrono::duration<float> time_mat = std::chrono::duration<float>(0);
    begin = std::chrono::system_clock::now();

    try
    {
        if (recover_illuminant)
        {
            // Perform illuminant recovery
            pipeline.process_recover_illuminant();
            end = std::chrono::system_clock::now();
            time_ill = end - begin;
        }
        if (recover_dichromatic)
        {
            // Perform dichrmatic parameter recovery
            pipeline.process_recover_dichromatic_parameters();
            end = std::chrono::system_clock::now();
            time_dic = end - begin;
        }
        if (recover_materials)
        {
            // Perform material recovery
            pipeline.process_recover_materials();
            end = std::chrono::system_clock::now();
            time_mat = end - begin;
        }

        if (save)
        {
            // Save the file (savename is the specified location of output)
            pipeline.save(savename);
            std::cout << "Saved " << savename << std::endl;
        }

    } catch (std::exception & e)
    {
        // There was an error!
        std::cout << e.what() << std::endl;
        return -1;
    }

    // Show images and data if requested
    if (show_images)
    {
        // Demonstrate the file was loaded.
        cv::Mat tmp[3];
        // It can sometimes be better to store the Image locally if you are going to making lots of accesses, as it may be generated each time
        // it is accessed (If you are reading from a HSZ file, for instance)
        arma::fcube I = pipeline.I();
        tmp[0] = scyl::arma_to_ocv(I.slice(1).t())/I.slice(1).max();
        tmp[1] = scyl::arma_to_ocv(I.slice(pipeline.bands() / 2).t())/I.slice(pipeline.bands() / 2).max();
        tmp[2] = scyl::arma_to_ocv(I.slice(pipeline.bands() - 1).t())/I.slice(pipeline.bands() - 2).max();
        cv::Mat showI;
        cv::merge(tmp, 3, showI);
        cv::namedWindow("I", CV_WINDOW_FREERATIO);
        cv::imshow("I", showI);

        try
        {
           std::cout << "Recovered Illuminant" << std::endl << pipeline.illuminant() << std::endl;
        } catch (std::exception & e) {}

        try
        {
            // Show the material map
            cv::Mat tmp2 = scyl::arma_to_ocv(arma::conv_to<arma::fmat>::from(pipeline.material_map().t())) / pipeline.material_map().max();
            tmp2.convertTo(tmp2, CV_8UC1, 255.0);
            cv::Mat showM;
            cv::applyColorMap(tmp2, showM, cv::COLORMAP_RAINBOW);
            cv::namedWindow("Material Map", CV_WINDOW_FREERATIO);
            cv::imshow("Material Map", showM);
        } catch (std::exception & e) {
            std::cout << e.what() << std::endl;
        }

        try
        {
            // Show the specular image.
            cv::Mat showg = scyl::arma_to_ocv(pipeline.g().t());
            cv::namedWindow("Shading", CV_WINDOW_FREERATIO);
            cv::imshow("Shading", showg/pipeline.g().max());
        } catch (std::exception & e) {}

        cv::waitKey(0);
    }

    //Save logs, print timing information
    scyl::logger::access()->save_logs("cli_pipe.log");
    if(keep_time)
    {
        std::cout << ">-------------------------------------------------------<" << std::endl
                  << "Image:            " << filename << " (" << pipeline.height() << "x" << pipeline.width() << "x" << pipeline.bands() << ")" << std::endl
                  << "Illuminant Time:  " << time_ill.count() << std::endl
                  << "Dichromatic Time: " << time_dic.count() - time_ill.count() << std::endl
                  << "Material Time:    " << time_mat.count() - time_dic.count() << std::endl
                  << "Total Time:       " << time_mat.count() << std::endl
                  << ">-------------------------------------------------------<" << std::endl;
    }
    // Calculate and print error
    if (calc_error)
    {
        try
        {
            arma::fcube reconstructed_s = scyl::reconstruct_s_elements(pipeline.material_elements(),pipeline.material_element_abundances(), pipeline.material_element_abundance_indexes());
            scyl::nan_to_zero(reconstructed_s);
            arma::fcube reconstructed_I = scyl::reconstruct_I(pipeline.illuminant(), pipeline.k(), pipeline.g(), reconstructed_s, pipeline.k_factor());
            scyl::nan_to_zero(reconstructed_I);

            reconstructed_I -= pipeline.I();

            float rms = arma::accu(arma::sqrt(arma::square(reconstructed_I)))/reconstructed_I.n_elem;

            std::cout << "RMS Error: " << rms << " (max = " << pipeline.I().max() << "), " << rms/pipeline.I().max()*100.0 << "%" << std::endl;

        } catch (std::exception & e) {}
    }

    return 0;

} // End program.