utility_imp.h

Go to the documentation of this file.

/****************************************************************************************
 * SCYLLARUS : C++ HYPERSPECTRAL PROCESSING LIBRARY
 * utility_imp.h - Implementations of Template functions
 *
 * This computer code is subject to copyright:
 * (c) National ICT Australia Limited (NICTA) 2013-2014 All Rights Reserved.
 *
 * Jeremy Oorloff, National ICT Australia (NICTA)
 *
 ***************************************************************************************/

#ifndef UTILITY_IMP_H_
#define UTILITY_IMP_H_

#include <armadillo>

namespace scyl
{

    template <class T> void nan_to_zero(T & in)
    {
        for(unsigned long i = 0; i < in.n_elem; i++)
        {
            if (!arma::is_finite(in(i)))
            {
                in(i) = 0;
            }
        }

    }


    template <class T> void zero_to_one(T & in)
    {
        float eps=arma::Datum<float>::eps;
        for(unsigned long i = 0; i < in.n_elem; i++)
        {
            if (in(i) <= 0.0f + eps and in(i) >= 0.0f - eps)
            {
                in(i) = 1.0f;
            }
        }
    }


    template <class T> int arma_cmp_num(const T & a,
                                        const T & b,
                                        float eps=arma::Datum<float>::eps,
                                        bool verbose=false)
    {
        unsigned int total_count = 0;
        unsigned int count = 0;
        unsigned int first_five = 0;
        if (a.n_elem != b.n_elem)
        {
            return a.n_elem; // Non matching objects, all elements are 'incorrect'
        }
        for (unsigned int i = 0; i < a.n_elem; i++)
        {
            if (!(a(i) <= b(i) + eps and a(i) >= b(i) - eps))
            {
                count++;
                if (verbose && first_five <= 5)
                {
                    std::cout.precision(15);
                    std::cout << "diff pixel " << i << " (" << a(i) << " -> " << b(i) << ")" << std::endl;

                    if (first_five == 5)
                    {
                        std::cout << "etc..." << std::endl;
                    }
                    first_five += 1;
                }
            }
            total_count++;
        }

        if (verbose)
        {
            std::cout << count << "/" << total_count << " elements differ ("
                      << (static_cast<float>(count)/total_count)*100.0 << "%)" << std::endl;
            T rms = arma::sqrt(arma::square(a) - arma::square(b));
            std::cout << "Max difference was: " << rms.max() << "(epsilon " << eps << ")" << std::endl;
        }
        return count;
    }

    template <class T> bool arma_cmp(const T & a,
                                     const T & b,
                                     const bool verbose=false,
                                     float eps=arma::Datum<float>::eps)
    {
        if (a.n_elem != b.n_elem)
        {
            return false;
        }
        for (unsigned int i = 0; i < a.n_elem; i++)
        {
            if (!(arma::is_finite(a(i))) && !(arma::is_finite(b(i))))
            {
                // Both nan or inf, so are the same
                continue;
            }
            if (!(a(i) <= b(i) + eps and a(i) >= b(i) - eps))
            {
                arma_cmp_num(a, b, eps, verbose);
                return false;
            }
        }
        return true;
    }

    template <class T> float arma_cmp_sig_figs(const T & a,
                                               const T & b,
                                               int num_sig_figs=6)
    {
        unsigned int total_count = 0;
        unsigned int count = 0;
        unsigned int first_five = 0;
        bool failed = false;
        if (a.n_elem != b.n_elem)
        {
            return 100.0;
        }
        for (unsigned int i = 0; i < a.n_elem; i++)
        {
            double scale = 0;
            if (a(i) != 0)
            {
                scale = pow(10.0, -1*ceil(log10(fabs(static_cast<double>(a(i))))));
            }
            double tolerance = pow(10.0, -1*num_sig_figs);
            double a_scaled = a(i)*scale;
            double b_scaled = b(i)*scale;

            if (!(a_scaled <= b_scaled + tolerance and a_scaled >= b_scaled - tolerance))
            {
                // Print information about which pixels have failed.
                count++;
                if (first_five <= 5)
                {
                    std::cout.precision(15);
                    std::cout << "Diff pixel " << i << " (" << a(i) << " -> " << b(i) << ")" << std::endl;

                    if (first_five == 5)
                    {
                        std::cout << "etc..." << std::endl;
                    }
                    first_five += 1;
                }
                failed = true;
            }
            total_count++;
        }
        if (failed)
        {
            float percent_count = (static_cast<float>(count)/total_count)*100.0;
            std::cout << count << "/" << total_count << " elements differ ("
                      << percent_count << "%)" << std::endl;
            T rms = arma::sqrt(arma::square(a) - arma::square(b));
            std::cout << "Max difference was: " << rms.max() << "(sig figs " << num_sig_figs << ")" << std::endl;
            return percent_count;
        }
        return 0.0;
    }
}

#endif /* UTILITY_IMP_H_ */