average.cxx

Go to the documentation of this file.

/*
         Copyright 2008-2019 Pierre SMARS (smars@yuntech.edu.tw)
         This program is free software: you can redistribute it and/or modify
         it under the terms of the GNU General Public License as published by
         the Free Software Foundation, either version 2 of the License, or
         (at your option) any later version.
 
         This program is distributed in the hope that it will be useful,
         but WITHOUT ANY WARRANTY; without even the implied warranty of
         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
         GNU General Public License for more details.
 
         You should have received a copy of the GNU General Public License
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
         */
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
 
#include "a_point.h"
#include "a_segment.h"
#include "a_triangle.h"
#include "a_mat_sq.h"
#include "a_mat_r.h"
#include "a_mat_c.h"
#include "SimpleOpt.h"
 
enum { OPT_HELP, OPT_P, OPT_L, OPT_T, OPT_BOX, OPT_VERBOSE, OPT_COV, OPT_SCALAR, OPT_MASS };
 
CSimpleOpt::SOption g_rgOptions[] = {
        { OPT_P, _T("-p"),     SO_NONE    },
        { OPT_P, _T("--point"),     SO_NONE    },
        { OPT_L, _T("-l"),     SO_NONE    },
        { OPT_L, _T("--line"),     SO_NONE    },
        { OPT_T, _T("-t"),     SO_NONE    },
        { OPT_T, _T("--triangle"),     SO_NONE    },
        { OPT_BOX, _T("-b"),     SO_NONE    },
        { OPT_BOX, _T("--box"),     SO_NONE    },
        { OPT_MASS, _T("-m"),     SO_NONE    },
        { OPT_MASS, _T("--mass"),     SO_NONE    },
        { OPT_VERBOSE, _T("-v"),     SO_NONE    },
        { OPT_VERBOSE, _T("--verbose"),     SO_NONE    },
        { OPT_COV, _T("-c"),     SO_NONE    },
        { OPT_COV, _T("--covariance"),     SO_NONE    },
        { OPT_SCALAR, _T("-sf"), SO_REQ_SEP },
        { OPT_SCALAR, _T("--scalar-file"), SO_REQ_SEP },
        { OPT_HELP, _T("-?"),     SO_NONE    },
        { OPT_HELP, _T("-h"),     SO_NONE    },
        { OPT_HELP, _T("--help"), SO_NONE    },
        SO_END_OF_OPTIONS
};
//*****************************************
int error(int val)
{
#include "average.help"
        return val;
}
//*****************************************
double xmin[] = {1.e30,1.e30,1.e30};
double xmax[] = {-1.e30,-1.e30,-1.e30};
a_mat_sq x2(3);
//*****************************************
void extremes(const a_point& p)
{
        for (short i=0;i<3;i++)
        {
                double x = p[i];
                if (x<xmin[i]) xmin[i]=x;
                if (x>xmax[i]) xmax[i]=x;
        }
}
//*****************************************
int main( int argc, char *argv[] )
{
        char code = 'p';
        bool box = false;
        bool verbose = false;
        bool cov = false;
        double mass;
        bool with_scalars = false;
        bool with_mass = false;
        std::string scalar_file;
        unsigned int n_scalars;
        a_point average;
        CSimpleOpt args(argc, argv, g_rgOptions);
        while (args.Next())
        {
                if (args.LastError() == SO_SUCCESS)
                {
                        if (args.OptionId() == OPT_HELP)
                                return error(0);
                        else if (args.OptionId() == OPT_P)
                                code = 'p';
                        else if (args.OptionId() == OPT_L)
                                code = 'l';
                        else if (args.OptionId() == OPT_T)
                                code = 't';
                        else if (args.OptionId() == OPT_BOX)
                                box = true;     
                        else if (args.OptionId() == OPT_COV)
                                cov = true;     
                        else if (args.OptionId() == OPT_SCALAR)
                        {
                                with_scalars = true;
                                scalar_file = args.OptionArg();
                        }
                        else if (args.OptionId() == OPT_MASS)
                                with_mass = true;       
                        else if (args.OptionId() == OPT_VERBOSE)
                                verbose = true; 
                        else
                                return error(-1);
                        // handle option, using OptionId(), OptionText() and OptionArg()
                } else {
                        return error(args.LastError());
                        // handle error, one of: SO_OPT_INVALID, SO_OPT_MULTIPLE,
                        // SO_ARG_INVALID, SO_ARG_INVALID_TYPE, SO_ARG_MISSING
                }
        }
        std::vector<double> scalars;
        if (with_scalars)
        {
                std::ifstream in(scalar_file.c_str());
                if (!in)
                {
                        std::cerr << "Cannot open file '" << scalar_file << "'" << std::endl;
                        return error(-2);
                }
                in >> n_scalars;
                for (unsigned int i=0; i<n_scalars; i++)
                {
                        double val;
                        in >> val;
                        scalars.push_back(val);
                }
                in.close();
        }
        a_point pt;
        switch (code) {
                case 'p':
                        {
                                unsigned int n_pts;
                                std::cin >> n_pts;
                                if ((with_scalars)&&(n_pts!=n_scalars))
                                {
                                        std::cerr << "Number of scalars does not correspond to number of points" << std::endl;
                                        return error(-2);
                                }
                                for (unsigned int k = 0; k < n_pts; k++)
                                {       
                                        a_point p;
                                        std::cin >> p;
                                        if (with_scalars)
                                        {
                                                mass += scalars[k];
                                                pt += p*scalars[k];
                                        }
                                        else
                                        {
                                                mass += 1;
                                                pt += p;
                                        }
                                        if (box)
                                        {
                                                if (!with_scalars||(scalars[k]!=0))
                                                        extremes(p);
                                        }
                                        if (cov)
                                        {
                                                if (with_scalars)
                                                        x2 += p.inertia()*scalars[k];
                                                else
                                                        x2 += p.inertia();
                                        }
                                }
                        }
                        break;
                case 'l':
                        {
                                int n_lis;
                                std::cin >> n_lis;
                                if ((with_scalars)&&(n_lis!=n_scalars))
                                {
                                        std::cerr << "Number of scalars does not correspond to number of lines" << std::endl;
                                        return error(-2);
                                }
                                for (unsigned int k = 0; k < n_lis; k++)
                                {
                                        a_segment s;
                                        std::cin >> s;
                                        double l = s.length();
                                        if (with_scalars)
                                                l *= scalars[k];
                                        mass += l;
                                        a_point centre = s.c();
                                        pt += centre*l;
                                        if (box)
                                        {
                                                if (!with_scalars||(scalars[k]!=0))
                                                {
                                                        extremes(s.p1());
                                                        extremes(s.p2());
                                                }
                                        }
                                        if (cov)
                                        {
                                                if (with_scalars)
                                                        x2 += s.inertia()*scalars[k];
                                                else
                                                        x2 += s.inertia();
                                        }
                                }
                        }
                        break;
                case 't':
                        {
                                int npts;
                                std::cin >> npts;
                                std::vector<a_point> points(npts);
                                for (unsigned int k = 0; k < npts; k++)
                                {
                                        a_point  p;
                                        std::cin >> p;
                                        if (box)
                                                extremes(p);
                                        points[k] = p;
                                }
                                int n_tris;
                                std::cin >> n_tris;
                                if ((with_scalars)&&(n_tris!=n_scalars))
                                {
                                        std::cerr << "Number of scalars does not correspond to number of triangles" << std::endl;
                                        return error(-2);
                                }
                                for (unsigned int k = 0; k < n_tris; k++)
                                {
                                        int pr[3];
                                        for (short l = 0; l < 3; l++)
                                                std::cin >> pr[l];
                                        a_triangle tri(points[pr[0]],points[pr[1]],points[pr[2]]);
                                        double surface = tri.S();
                                        if (with_scalars)
                                                surface *= scalars[k];
                                        a_point centre = tri.c();
                                        mass += surface;
                                        pt += centre*surface;
                                        if (cov)
                                        {
                                                if (with_scalars)
                                                        x2 += tri.inertia()*scalars[k];
                                                else
                                                        x2 += tri.inertia();
                                        }
                                        points.clear();
                                }
                                break;
                        }
        }
        average = pt/mass;
        if (verbose)
                std::cout << "center_of_mass: ";
        std::cout << average << std::endl;
        if (with_mass)
        {
                if (verbose)
                        std::cout << "mass: ";
                std::cout << mass << std::endl;
        }
        if (cov)
        {
                if (verbose)
                        std::cout << "inertia: " << std::endl;
                a_mat_sq m(3);
                m = average.inertia();
                m *= mass;
                std::cout << x2-m;
        }
        if (box)
        {
                if (verbose)
                {
                        std::cout << "bounding_box: ";
                        std::cout << std::endl << "min: ";
                }
                std::cout << xmin[0] << " " << xmin[1] << " " << xmin[2] << std::endl;
                if (verbose)
                        std::cout << "max: ";
                std::cout << xmax[0] << " " << xmax[1] << " " << xmax[2] << std::endl;
        }
        return 0;
}