distanceMapGenerationTool.cpp

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// -                        CloudViewer: www.cloudViewer.org                  -
// ----------------------------------------------------------------------------
// Copyright (c) 2018-2024 www.cloudViewer.org
// SPDX-License-Identifier: MIT
// ----------------------------------------------------------------------------
 
#include "distanceMapGenerationTool.h"
 
// qCC
#include <ecvMainAppInterface.h>
 
// CV_DB_LIB
#include <ecvMesh.h>
#include <ecvPointCloud.h>
#include <ecvPolyline.h>
#include <ecvProgressDialog.h>
#include <ecvScalarField.h>
 
// cloudViewer
#include <Delaunay2dMesh.h>
 
// Qt
#include <QFile>
#include <QMainWindow>
// Qt5/Qt6 Compatibility
#include <QtCompat.h>
 
// Meta-data key for profile (polyline) origin
const char PROFILE_ORIGIN_KEY[] = "ProfileOrigin";
// Meta-data key for profile (polyline) axis
const char REVOLUTION_AXIS_KEY[] = "RevolutionAxis";
// Meta-data key for the profile (polyline) height shift
const char PROFILE_HEIGHT_SHIFT_KEY[] = "ProfileHeightShift";
 
// shortcuts
static const double M_PI_DIV_2 = M_PI / 2;
static const double M_PI_DIV_4 = M_PI / 4;
 
// helper
static inline double ComputeLatitude_rad(PointCoordinateType x,
                                         PointCoordinateType y,
                                         PointCoordinateType z) {
    double r = static_cast<double>(x * x + y * y);
 
    if (r < 1.0e-8) {
        return z < 0.0 ? -M_PI_DIV_2 : M_PI_DIV_2;
    }
 
    return atan(z / sqrt(static_cast<double>(r)));
}
 
// helper
static bool GetPolylineMetaVector(const ccPolyline* polyline,
                                  const QString& key,
                                  CCVector3& P) {
    assert(polyline);
    if (polyline) {
        // we try to get the right meta-data 'vector'
        QVariant x = polyline->getMetaData(key + QString(".x"));
        QVariant y = polyline->getMetaData(key + QString(".y"));
        QVariant z = polyline->getMetaData(key + QString(".z"));
        if (x.isValid() && y.isValid() && z.isValid()) {
            bool ok[3] = {true, true, true};
            P.x = static_cast<PointCoordinateType>(x.toDouble(ok));
            P.y = static_cast<PointCoordinateType>(y.toDouble(ok + 1));
            P.z = static_cast<PointCoordinateType>(z.toDouble(ok + 2));
            return (ok[0] && ok[1] && ok[2]);
        }
    }
 
    return false;
}
 
// helper
static void SetPoylineMetaVector(ccPolyline* polyline,
                                 const QString& key,
                                 const CCVector3& P) {
    assert(polyline);
    if (polyline) {
        // add revolution axis as meta-data
        polyline->setMetaData(key + QString(".x"), QVariant(P.x));
        polyline->setMetaData(key + QString(".y"), QVariant(P.y));
        polyline->setMetaData(key + QString(".z"), QVariant(P.z));
    }
}
 
// helper
ccGLMatrix
DistanceMapGenerationTool::ProfileMetaData::computeCloudToSurfaceOriginTrans()
        const {
    ccGLMatrix cloudToSurfaceOrigin;
    cloudToSurfaceOrigin.setTranslation(-origin);
 
    // we must take the axis of the Surface of Revolution into account (if any
    // and if it is not colinear with X, Y or Z)
    if (hasAxis &&
        axis.u[revolDim] + std::numeric_limits<PointCoordinateType>::epsilon() <
                1.0) {
        ccGLMatrix rotation;
        CCVector3 Z(0, 0, 0);
        Z.u[revolDim] = PC_ONE;
        rotation = ccGLMatrix::FromToRotation(axis, Z);
        cloudToSurfaceOrigin = rotation * cloudToSurfaceOrigin;
    }
 
    return cloudToSurfaceOrigin;
}
 
// helper
ccGLMatrix
DistanceMapGenerationTool::ProfileMetaData::computeCloudToProfileOriginTrans()
        const {
    ccGLMatrix cloudToPolylineOrigin = computeCloudToSurfaceOriginTrans();
 
    // add the height shift along the revolution axis (if any)
    cloudToPolylineOrigin.getTranslation()[revolDim] -= heightShift;
 
    return cloudToPolylineOrigin;
}
 
void DistanceMapGenerationTool::SetPoylineRevolDim(ccPolyline* polyline,
                                                   int revolDim) {
    assert(polyline);
    if (polyline) {
        // add revolution dimension as meta-data
        QVariant dim(revolDim);
        polyline->setMetaData(REVOLUTION_AXIS_KEY, dim);
    }
}
 
int DistanceMapGenerationTool::GetPoylineRevolDim(const ccPolyline* polyline) {
    assert(polyline);
    if (polyline) {
        // we try to get the revolution dimension from the polyline meta-data
        QVariant axis = polyline->getMetaData(REVOLUTION_AXIS_KEY);
        if (axis.isValid()) {
            bool ok = true;
            int dim = axis.toInt(&ok);
            if (ok && dim >= 0 && dim <= 2) return dim;
        }
    }
 
    return -1;
}
 
void DistanceMapGenerationTool::SetPoylineOrigin(ccPolyline* polyline,
                                                 const CCVector3& origin) {
    SetPoylineMetaVector(polyline, PROFILE_ORIGIN_KEY, origin);
}
 
bool DistanceMapGenerationTool::GetPoylineOrigin(const ccPolyline* polyline,
                                                 CCVector3& origin) {
    // we try to get the profile origin from the polyline meta-data
    return GetPolylineMetaVector(polyline, PROFILE_ORIGIN_KEY, origin);
}
 
void DistanceMapGenerationTool::SetPoylineAxis(ccPolyline* polyline,
                                               const CCVector3& axis) {
    SetPoylineMetaVector(polyline, REVOLUTION_AXIS_KEY, axis);
}
 
bool DistanceMapGenerationTool::GetPoylineAxis(const ccPolyline* polyline,
                                               CCVector3& axis) {
    // we try to get the profile origin from the polyline meta-data
    return GetPolylineMetaVector(polyline, REVOLUTION_AXIS_KEY, axis);
}
 
void DistanceMapGenerationTool::SetPolylineHeightShift(
        ccPolyline* polyline, PointCoordinateType heightShift) {
    assert(polyline);
    if (polyline) {
        // add 'height shift' as meta-data
        polyline->setMetaData(PROFILE_HEIGHT_SHIFT_KEY, QVariant(heightShift));
    }
}
 
bool DistanceMapGenerationTool::GetPolylineHeightShift(
        const ccPolyline* polyline, PointCoordinateType& heightShift) {
    assert(polyline);
    if (polyline) {
        // retrieve the right meta-data
        QVariant shift = polyline->getMetaData(PROFILE_HEIGHT_SHIFT_KEY);
        if (shift.isValid()) {
            bool ok;
            heightShift = static_cast<PointCoordinateType>(shift.toDouble(&ok));
            return ok;
        }
    }
 
    return false;
}
 
bool DistanceMapGenerationTool::GetPoylineMetaData(const ccPolyline* polyline,
                                                   ProfileMetaData& data) {
    if (!polyline) {
        assert(false);
        return false;
    }
 
    data.revolDim = GetPoylineRevolDim(polyline);
    if (data.revolDim < 0 || data.revolDim > 2) {
        return false;
    }
 
    if (!GetPoylineOrigin(polyline, data.origin)) {
        return false;
    }
 
    if (!GetPolylineHeightShift(polyline, data.heightShift)) {
        data.heightShift = 0;
    }
 
    data.hasAxis = GetPoylineAxis(polyline, data.axis);
 
    return true;
}
 
bool DistanceMapGenerationTool::ComputeRadialDist(
        ccPointCloud* cloud,
        ccPolyline* profile,
        bool storeRadiiAsSF /*=false*/,
        ecvMainAppInterface* app /*=0*/) {
    // check input cloud and profile/polyline
    if (!cloud || !profile) {
        if (app)
            app->dispToConsole(
                    QString("Internal error: invalid input parameters"),
                    ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return false;
    }
    assert(cloud && profile);
 
    // number of vertices for the profile
    cloudViewer::GenericIndexedCloudPersist* vertices =
            profile->getAssociatedCloud();
    unsigned vertexCount = vertices->size();
    if (vertexCount < 2) {
        if (app)
            app->dispToConsole(
                    QString("Invalid polyline (not enough vertices)"),
                    ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return false;
    }
 
    // profile meta-data
    ProfileMetaData profileDesc;
    if (!GetPoylineMetaData(profile, profileDesc)) {
        if (app)
            app->dispToConsole(
                    QString("Invalid polyline (bad or missing meta-data)"),
                    ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return false;
    }
 
    // reserve a new scalar field (or take the old one if it already exists)
    int sfIdx = cloud->getScalarFieldIndexByName(RADIAL_DIST_SF_NAME);
    if (sfIdx < 0) sfIdx = cloud->addScalarField(RADIAL_DIST_SF_NAME);
    if (sfIdx < 0) {
        if (app)
            app->dispToConsole(
                    QString("Failed to allocate a new scalar field for "
                            "computing distances! Try to free some memory ..."),
                    ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return false;
    }
    ccScalarField* sf =
            static_cast<ccScalarField*>(cloud->getScalarField(sfIdx));
    unsigned pointCount = cloud->size();
    sf->resize(pointCount);  // should always be ok
    assert(sf);
 
    ccScalarField* radiiSf = 0;
    if (storeRadiiAsSF) {
        int sfIdxRadii = cloud->getScalarFieldIndexByName(RADII_SF_NAME);
        if (sfIdxRadii < 0) sfIdxRadii = cloud->addScalarField(RADII_SF_NAME);
        if (sfIdxRadii < 0) {
            if (app)
                app->dispToConsole(
                        QString("Failed to allocate a new scalar field for "
                                "storing radii! You should try to free some "
                                "memory ..."),
                        ecvMainAppInterface::WRN_CONSOLE_MESSAGE);
            // return false;
        } else {
            radiiSf = static_cast<ccScalarField*>(
                    cloud->getScalarField(sfIdxRadii));
            radiiSf->resize(pointCount);  // should always be ok
        }
    }
 
    bool success = true;
 
    // now compute the distance between the cloud and the (implicit) surface of
    // revolution
    {
        ccGLMatrix cloudToProfile =
                profileDesc.computeCloudToProfileOriginTrans();
 
        // we deduce the horizontal dimensions from the revolution axis
        const unsigned char dim1 = static_cast<unsigned char>(
                profileDesc.revolDim < 2 ? profileDesc.revolDim + 1 : 0);
        const unsigned char dim2 = (dim1 < 2 ? dim1 + 1 : 0);
 
        ecvProgressDialog dlg(true, app ? app->getMainWindow() : 0);
        dlg.setMethodTitle(QObject::tr("Cloud to profile radial distance"));
        dlg.setInfo(QObject::tr("Polyline: %1 vertices\nCloud: %2 points")
                            .arg(vertexCount)
                            .arg(pointCount));
        dlg.start();
        cloudViewer::NormalizedProgress nProgress(
                static_cast<cloudViewer::GenericProgressCallback*>(&dlg),
                pointCount);
 
        for (unsigned i = 0; i < pointCount; ++i) {
            const CCVector3* P = cloud->getPoint(i);
 
            // relative point position
            CCVector3 Prel = cloudToProfile * (*P);
 
            // deduce point height and radius (i.e. in profile 2D coordinate
            // system)
            double height = Prel.u[profileDesc.revolDim];
            // TODO FIXME: we assume the surface of revolution is smooth!
            double radius = sqrt(Prel.u[dim1] * Prel.u[dim1] +
                                 Prel.u[dim2] * Prel.u[dim2]);
 
            if (radiiSf) {
                ScalarType radiusVal = static_cast<ScalarType>(radius);
                radiiSf->setValue(i, radiusVal);
            }
 
            // search nearest "segment" in polyline
            ScalarType minDist = NAN_VALUE;
            for (unsigned j = 1; j < vertexCount; ++j) {
                const CCVector3* A = vertices->getPoint(j - 1);
                const CCVector3* B = vertices->getPoint(j);
 
                double alpha = (height - A->y) / (B->y - A->y);
                if (alpha >= 0.0 && alpha <= 1.0) {
                    // we deduce the right radius by linear interpolation
                    double radius_th = A->x + alpha * (B->x - A->x);
                    double dist = radius - radius_th;
 
                    // we look at the closest segment (if the polyline is
                    // concave!)
                    if (!cloudViewer::ScalarField::ValidValue(minDist) ||
                        dist * dist < minDist * minDist) {
                        minDist = static_cast<ScalarType>(dist);
                    }
                }
            }
 
            sf->setValue(i, minDist);
 
            if (!nProgress.oneStep()) {
                // cancelled by user
                for (unsigned j = i; j < pointCount; ++j)
                    sf->setValue(j, NAN_VALUE);
 
                success = false;
                break;
            }
 
            // TEST
            //*const_cast<CCVector3*>(P) = Prel;
        }
 
        // TEST
        // cloud->invalidateBoundingBox();
    }
 
    sf->computeMinAndMax();
    cloud->setCurrentDisplayedScalarField(sfIdx);
    cloud->showSF(true);
 
    return success;
}
 
bool DistanceMapGenerationTool::ComputeMinAndMaxLatitude_rad(
        ccPointCloud* cloud,
        double& minLat_rad,
        double& maxLat_rad,
        const ccGLMatrix& cloudToSurfaceOrigin,  // e.g. translation to the
                                                 // revolution origin
        unsigned char revolutionAxisDim) {
    minLat_rad = maxLat_rad = 0.0;
 
    assert(cloud);
    // invalid parameters?
    if (!cloud || revolutionAxisDim > 2) return false;
 
    unsigned count = cloud->size();
    if (count == 0) return true;
 
    // revolution axis
    const unsigned char Z = revolutionAxisDim;
    // we deduce the 2 other ('horizontal') dimensions from the revolution axis
    const unsigned char X = (Z < 2 ? Z + 1 : 0);
    const unsigned char Y = (X < 2 ? X + 1 : 0);
 
    for (unsigned n = 0; n < count; ++n) {
        const CCVector3* P = cloud->getPoint(n);
        CCVector3 relativePos = cloudToSurfaceOrigin * (*P);
 
        // latitude between 0 and pi/2
        double lat_rad = ComputeLatitude_rad(relativePos.u[X], relativePos.u[Y],
                                             relativePos.u[Z]);
 
        if (n) {
            if (lat_rad < minLat_rad)
                minLat_rad = lat_rad;
            else if (lat_rad > maxLat_rad)
                maxLat_rad = lat_rad;
        } else {
            minLat_rad = maxLat_rad = lat_rad;
        }
    }
 
    return true;
}
 
double DistanceMapGenerationTool::ConicalProjectN(double phi1, double phi2) {
    if (phi1 >= phi2) return 1.0;
 
    assert(fabs(phi1) < M_PI_DIV_2);
    assert(fabs(phi2) < M_PI_DIV_2);
 
    double tan_pl1 = tan(M_PI_DIV_4 - phi1 / 2);
    double tan_pl2 = tan(M_PI_DIV_4 - phi2 / 2);
 
    return (log(cos(phi1)) - log(cos(phi2))) / (log(tan_pl1) - log(tan_pl2));
}
 
double DistanceMapGenerationTool::ConicalProject(double phi,
                                                 double phi1,
                                                 double n) {
    assert(phi1 >= -M_PI_DIV_2 && phi1 < M_PI_DIV_2);
 
    double tan_pl1 = tan(M_PI_DIV_4 - phi1 / 2);
    double tan_pl = tan(M_PI_DIV_4 - phi / 2);
 
    return cos(phi1) * pow(tan_pl / tan_pl1, n) / n;
}
 
QSharedPointer<DistanceMapGenerationTool::Map>
DistanceMapGenerationTool::CreateMap(ccPointCloud* cloud,
                                     ccScalarField* sf,
                                     const ccGLMatrix& cloudToSurface,
                                     unsigned char revolutionAxisDim,
                                     double xStep_rad,
                                     double yStep,
                                     double yMin,
                                     double yMax,
                                     bool conical,
                                     bool counterclockwise,
                                     FillStrategyType fillStrategy,
                                     EmptyCellFillOption emptyCellfillOption,
                                     ecvMainAppInterface* app /*=0*/) {
    assert(cloud && sf);
    if (!cloud || !sf) {
        if (app)
            app->dispToConsole(QString("[DistanceMapGenerationTool] Internal "
                                       "error: invalid input structures!"),
                               ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return QSharedPointer<Map>(0);
    }
 
    // invalid parameters?
    if (xStep_rad <= 0.0 || yStep <= 0.0 || yMax <= yMin ||
        revolutionAxisDim > 2) {
        if (app)
            app->dispToConsole(QString("[DistanceMapGenerationTool] Internal "
                                       "error: invalid grid parameters!"),
                               ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return QSharedPointer<Map>(0);
    }
 
    unsigned count = cloud->size();
    if (count == 0) {
        if (app)
            app->dispToConsole(QString("[DistanceMapGenerationTool] Cloud is "
                                       "empty! Nothing to do!"),
                               ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return QSharedPointer<Map>(0);
    }
 
    // revolution axis
    const unsigned char Z = revolutionAxisDim;
    // we deduce the 2 other ('horizontal') dimensions from the revolution axis
    const unsigned char X = (Z < 2 ? Z + 1 : 0);
    const unsigned char Y = (X < 2 ? X + 1 : 0);
 
    // grid dimensions
    unsigned xSteps = 0;
    {
        if (xStep_rad > 0)
            xSteps = static_cast<unsigned>(ceil((2 * M_PI) / xStep_rad));
        if (xSteps == 0) {
            if (app)
                app->dispToConsole(QString("[DistanceMapGenerationTool] "
                                           "Invalid longitude step/boundaries! "
                                           "Can't generate a proper map!"),
                                   ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
            return QSharedPointer<Map>(0);
        }
    }
 
    unsigned ySteps = 0;
    {
        if (yStep > 0)
            ySteps = static_cast<unsigned>(ceil((yMax - yMin) / yStep));
        if (ySteps == 0) {
            if (app)
                app->dispToConsole(QString("[DistanceMapGenerationTool] "
                                           "Invalid latitude step/boundaries! "
                                           "Can't generate a proper map!"),
                                   ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
            return QSharedPointer<Map>(0);
        }
    }
 
    unsigned cellCount = xSteps * ySteps;
    if (app)
        app->dispToConsole(QString("[DistanceMapGenerationTool] Projected map "
                                   "size: %1 x %2 (%3 cells)")
                                   .arg(xSteps)
                                   .arg(ySteps)
                                   .arg(cellCount),
                           ecvMainAppInterface::STD_CONSOLE_MESSAGE);
 
    // reserve memory for the output matrix
    QSharedPointer<Map> grid(new Map);
    try {
        grid->resize(cellCount);
    } catch (const std::bad_alloc&) {
        if (app)
            app->dispToConsole(
                    QString("[DistanceMapGenerationTool] Not enough memory!"),
                    ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return QSharedPointer<Map>(0);
    }
 
    // update grid info ("for the records")
    grid->xSteps = xSteps;
    grid->xMin = 0.0;
    grid->xMax = 2 * M_PI;
    grid->xStep = xStep_rad;
    grid->ySteps = ySteps;
    grid->yMin = yMin;
    grid->yMax = yMax;
    grid->yStep = yStep;
    grid->conical = conical;
 
    // motion direction
    grid->counterclockwise = counterclockwise;
    double ccw = (counterclockwise ? -1.0 : 1.0);
 
    for (unsigned n = 0; n < count; ++n) {
        // we skip invalid values
        const ScalarType& val = sf->getValue(n);
        if (!cloudViewer::ScalarField::ValidValue(val)) continue;
 
        const CCVector3* P = cloud->getPoint(n);
        CCVector3 relativePos = cloudToSurface * (*P);
 
        // convert to cylindrical or conical (spherical) coordinates
        double x =
                ccw * atan2(relativePos.u[X], relativePos.u[Y]);  // longitude
        if (x < 0.0) {
            x += 2 * M_PI;
        }
 
        double y = 0.0;
        if (conical) {
            y = ComputeLatitude_rad(
                    relativePos.u[X], relativePos.u[Y],
                    relativePos.u[Z]);  // latitude between 0 and pi/2
        } else {
            y = relativePos.u[Z];  // height
        }
 
        int i = static_cast<int>((x - grid->xMin) / grid->xStep);
        int j = static_cast<int>((y - grid->yMin) / grid->yStep);
 
        // if we fall exactly on the max corner of the grid box
        if (i == static_cast<int>(grid->xSteps)) --i;
        if (j == static_cast<int>(grid->ySteps)) --j;
 
        // we skip points outside the box!
        if (i < 0 || i >= static_cast<int>(grid->xSteps) || j < 0 ||
            j >= static_cast<int>(grid->ySteps)) {
            continue;
        }
        assert(i >= 0 && j >= 0);
 
        MapCell& cell = (*grid)[j * static_cast<int>(grid->xSteps) + i];
        if (cell.count)  // if there's already values projected in this cell
        {
            switch (fillStrategy) {
                case FILL_STRAT_MIN_DIST:
                    // Set the minimum SF value
                    if (val < cell.value) cell.value = val;
                    break;
                case FILL_STRAT_AVG_DIST:
                    // Sum the values
                    cell.value += static_cast<double>(val);
                    break;
                case FILL_STRAT_MAX_DIST:
                    // Set the maximum SF value
                    if (val > cell.value) cell.value = val;
                    break;
                default:
                    assert(false);
                    break;
            }
        } else {
            // for the first point, we simply have to store its associated value
            // (whatever the case)
            cell.value = val;
        }
        ++cell.count;
 
        // if (progressCb)
        //  progressCb->update(30.0 * (float)n / (float)cloud->size());
    }
 
    // we need to finish the average values computation
    if (fillStrategy == FILL_STRAT_AVG_DIST) {
        MapCell* cell = &grid->at(0);
        for (unsigned i = 0; i < cellCount; ++i, ++cell)
            if (cell->count > 1)
                cell->value /= static_cast<double>(cell->count);
    }
 
    // fill empty cells with zero?
    if (emptyCellfillOption == FILL_WITH_ZERO) {
        MapCell* cell = &grid->at(0);
        for (unsigned i = 0; i < cellCount; ++i, ++cell) {
            if (cell->count == 0) {
                cell->value = 0.0;
                cell->count = 1;
            }
        }
    } else if (emptyCellfillOption == FILL_INTERPOLATE) {
        // convert the non-empty cells to a 2D point cloud
        unsigned fillCount = 0;
        {
            MapCell* cell = &grid->at(0);
            for (unsigned i = 0; i < cellCount; ++i, ++cell)
                if (cell->count != 0) ++fillCount;
        }
 
        // do we really need to interpolate empty grid cells?
        if (fillCount) {
            std::vector<CCVector2> the2DPoints;
            try {
                the2DPoints.reserve(fillCount);
            } catch (...) {
                // out of memory
                if (app)
                    app->dispToConsole(
                            QString("[DistanceMapGenerationTool] Not enough "
                                    "memory to interpolate!"),
                            ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
            }
 
            if (the2DPoints.capacity() == fillCount) {
                // fill 2D vector with non-empty cell indexes
                {
                    const MapCell* cell = &grid->at(0);
                    for (unsigned j = 0; j < grid->ySteps; ++j)
                        for (unsigned i = 0; i < grid->xSteps; ++i, ++cell)
                            if (cell->count)
                                the2DPoints.push_back(CCVector2(
                                        static_cast<PointCoordinateType>(i),
                                        static_cast<PointCoordinateType>(j)));
                }
 
                // mesh the '2D' points
                cloudViewer::Delaunay2dMesh* dm =
                        new cloudViewer::Delaunay2dMesh();
                std::string errorStr;
                if (!dm->buildMesh(the2DPoints,
                                   cloudViewer::Delaunay2dMesh::USE_ALL_POINTS,
                                   errorStr)) {
                    if (app)
                        app->dispToConsole(
                                QString("[DistanceMapGenerationTool] "
                                        "Interpolation failed: Triangle lib. "
                                        "said '%1'")
                                        .arg(QString::fromStdString(errorStr)),
                                ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
                } else {
                    unsigned triNum = dm->size();
                    MapCell* cells = &grid->at(0);
                    // now we are going to 'project' all triangles on the grid
                    dm->placeIteratorAtBeginning();
                    for (unsigned k = 0; k < triNum; ++k) {
                        const cloudViewer::VerticesIndexes* tsi =
                                dm->getNextTriangleVertIndexes();
                        // get the triangle bounding box (in grid coordinates)
                        int P[3][2];
                        int xMin = 0, yMin = 0, xMax = 0, yMax = 0;
                        {
                            for (unsigned j = 0; j < 3; ++j) {
                                const CCVector2& P2D = the2DPoints[tsi->i[j]];
                                P[j][0] = static_cast<int>(P2D.x);
                                P[j][1] = static_cast<int>(P2D.y);
                            }
                            xMin = std::min(std::min(P[0][0], P[1][0]),
                                            P[2][0]);
                            yMin = std::min(std::min(P[0][1], P[1][1]),
                                            P[2][1]);
                            xMax = std::max(std::max(P[0][0], P[1][0]),
                                            P[2][0]);
                            yMax = std::max(std::max(P[0][1], P[1][1]),
                                            P[2][1]);
                        }
                        // now scan the cells
                        {
                            // pre-computation for barycentric coordinates
                            const double& valA =
                                    cells[P[0][0] + P[0][1] * grid->xSteps]
                                            .value;
                            const double& valB =
                                    cells[P[1][0] + P[1][1] * grid->xSteps]
                                            .value;
                            const double& valC =
                                    cells[P[2][0] + P[2][1] * grid->xSteps]
                                            .value;
                            int det =
                                    (P[1][1] - P[2][1]) * (P[0][0] - P[2][0]) +
                                    (P[2][0] - P[1][0]) * (P[0][1] - P[2][1]);
 
                            for (int j = yMin; j <= yMax; ++j) {
                                MapCell* cell =
                                        cells +
                                        static_cast<unsigned>(j) * grid->xSteps;
 
                                for (int i = xMin; i <= xMax; ++i) {
                                    // if the cell is empty
                                    if (!cell[i].count) {
                                        // we test if it's included or not in
                                        // the current triangle Point Inclusion
                                        // in Polygon Test (inspired from W.
                                        // Randolph Franklin - WRF)
                                        bool inside = false;
                                        for (int ti = 0; ti < 3; ++ti) {
                                            const int* P1 = P[ti];
                                            const int* P2 = P[(ti + 1) % 3];
                                            if ((P2[1] <= j && j < P1[1]) ||
                                                (P1[1] <= j && j < P2[1])) {
                                                int t = (i - P2[0]) * (P1[1] -
                                                                       P2[1]) -
                                                        (P1[0] - P2[0]) *
                                                                (j - P2[1]);
                                                if (P1[1] < P2[1]) t = -t;
                                                if (t < 0) inside = !inside;
                                            }
                                        }
                                        // can we interpolate?
                                        if (inside) {
                                            double l1 =
                                                    static_cast<double>(
                                                            (P[1][1] -
                                                             P[2][1]) *
                                                                    (i -
                                                                     P[2][0]) +
                                                            (P[2][0] -
                                                             P[1][0]) *
                                                                    (j -
                                                                     P[2][1])) /
                                                    det;
                                            double l2 =
                                                    static_cast<double>(
                                                            (P[2][1] -
                                                             P[0][1]) *
                                                                    (i -
                                                                     P[2][0]) +
                                                            (P[0][0] -
                                                             P[2][0]) *
                                                                    (j -
                                                                     P[2][1])) /
                                                    det;
                                            double l3 = 1.0 - l1 - l2;
 
                                            cell[i].count = 1;
                                            cell[i].value = l1 * valA +
                                                            l2 * valB +
                                                            l3 * valC;
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
 
                delete dm;
                dm = 0;
            }
        }
    }
 
    // update min and max values
    {
        const MapCell* cell = &grid->at(0);
        grid->minVal = grid->maxVal = cell->value;
        ++cell;
        for (unsigned i = 1; i < cellCount; ++i, ++cell) {
            if (cell->value < grid->minVal)
                grid->minVal = cell->value;
            else if (cell->value > grid->maxVal)
                grid->maxVal = cell->value;
        }
    }
 
    // end of process
    return grid;
}
 
CCVector3 DistanceMapGenerationTool::ProjectPointOnCone(double lon_rad,
                                                        double lat_rad,
                                                        double latMin_rad,
                                                        double nProj,
                                                        bool counterclockwise) {
    double theta = nProj * (lon_rad - M_PI);
    double r = ConicalProject(lat_rad, latMin_rad, nProj);
 
    CCVector3 P(static_cast<PointCoordinateType>((counterclockwise ? -r : r) *
                                                 sin(theta)),
                static_cast<PointCoordinateType>(-r * cos(theta)), 0);
 
    return P;
}
 
ccMesh* DistanceMapGenerationTool::ConvertConicalMapToMesh(
        const QSharedPointer<Map>& map,
        bool counterclockwise,
        QImage mapTexture /*=QImage()*/) {
    if (!map) return 0;
 
    unsigned meshVertCount = map->xSteps * map->ySteps;
    unsigned meshFaceCount = (map->xSteps - 1) * (map->ySteps - 1) * 2;
    ccPointCloud* cloud = new ccPointCloud();
    ccMesh* mesh = new ccMesh(cloud);
    mesh->addChild(cloud);
    if (!cloud->reserve(meshVertCount) || !mesh->reserve(meshFaceCount)) {
        // not enough memory
        delete mesh;
        return 0;
    }
 
    // compute projection constant
    double nProj =
            ConicalProjectN(map->yMin, map->yMax) * map->conicalSpanRatio;
    assert(nProj >= -1.0 && nProj <= 1.0);
 
    // create vertices
    {
        double cwSign = (counterclockwise ? -1.0 : 1.0);
        for (unsigned j = 0; j < map->xSteps; ++j) {
            // longitude
            double lon_rad =
                    static_cast<double>(j) / map->xSteps * (2.0 * M_PI);
 
            double theta = nProj *
                           (lon_rad -
                            M_PI);  //-Pi shift so that the map is well centered
            double sin_theta = sin(theta);
            double cos_theta = cos(theta);
 
            for (unsigned i = 0; i < map->ySteps; ++i) {
                double lat_rad =
                        map->yMin + static_cast<double>(i) * map->yStep;
                double r = ConicalProject(lat_rad, map->yMin, nProj);
 
                CCVector3 Pxyz(static_cast<PointCoordinateType>(cwSign * r *
                                                                sin_theta),
                               static_cast<PointCoordinateType>(-r * cos_theta),
                               0);
                cloud->addPoint(Pxyz);
            }
        }
    }
 
    // create facets
    {
        for (unsigned j = 0; j + 1 < map->xSteps; ++j) {
            for (unsigned i = 0; i + 1 < map->ySteps; ++i) {
                unsigned vertA = j * map->ySteps + i;
                unsigned vertB = vertA + map->ySteps;
                unsigned vertC = vertB + 1;
                unsigned vertD = vertA + 1;
 
                mesh->addTriangle(vertB, vertC, vertD);
                mesh->addTriangle(vertB, vertD, vertA);
            }
        }
    }
 
    // do we have a texture as well?
    if (true )  // we force tex. coordinates and indexes
                                        // creation!
    {
        // texture coordinates
        TextureCoordsContainer* texCoords = new TextureCoordsContainer();
        if (!texCoords->reserveSafe(meshVertCount)) {
            // not enough memory to finish the job!
            delete texCoords;
            return mesh;
        }
 
        // create default texture coordinates
        for (unsigned j = 0; j < map->xSteps; ++j) {
            TexCoords2D T(static_cast<float>(j) / (map->xSteps - 1), 0.0f);
            for (unsigned i = 0; i < map->ySteps; ++i) {
                T.ty = static_cast<float>(i) / (map->ySteps - 1);
                texCoords->addElement(T);
            }
        }
 
        if (!mesh->reservePerTriangleTexCoordIndexes()) {
            // not enough memory to finish the job!
            delete texCoords;
            return mesh;
        }
 
        // set texture indexes
        {
            for (unsigned j = 0; j + 1 < map->xSteps; ++j) {
                for (unsigned i = 0; i + 1 < map->ySteps; ++i) {
                    unsigned vertA = j * map->ySteps + i;
                    unsigned vertB = vertA + map->ySteps;
                    unsigned vertC = vertB + 1;
                    unsigned vertD = vertA + 1;
 
                    mesh->addTriangleTexCoordIndexes(vertB, vertC, vertD);
                    mesh->addTriangleTexCoordIndexes(vertB, vertD, vertA);
                }
            }
        }
 
        // set material indexes
        if (!mesh->reservePerTriangleMtlIndexes()) {
            // not enough memory to finish the job!
            delete texCoords;
            mesh->removePerTriangleTexCoordIndexes();
            return mesh;
        }
        for (unsigned i = 0; i < meshFaceCount; ++i) {
            mesh->addTriangleMtlIndex(0);
        }
 
        // set material
#if 0
        {
            ccMaterial::Shared material(new ccMaterial("texture"));
            material->setTexture(mapTexture, QString(), false);
 
            ccMaterialSet* materialSet = new ccMaterialSet();
            materialSet->addMaterial(material);
 
            mesh->setMaterialSet(materialSet);
        }
#endif
 
        mesh->setTexCoordinatesTable(texCoords);
        mesh->showMaterials(true);
        mesh->setVisible(true);
    }
 
    return mesh;
}
 
bool DistanceMapGenerationTool::ComputeSurfacesAndVolumes(
        const QSharedPointer<Map>& map,
        ccPolyline* profile,
        Measures& surface,
        Measures& volume) {
    if (!map || !profile)
        // invalid input!
        return false;
 
    cloudViewer::GenericIndexedCloudPersist* vertices =
            profile->getAssociatedCloud();
    unsigned vertexCount = vertices ? vertices->size() : 0;
    if (vertexCount < 2) {
        // invalid profile!
        return false;
    }
 
    const ccPointCloud* pcVertices =
            dynamic_cast<ccPointCloud*>(profile->getAssociatedCloud());
    if (!pcVertices) {
        return false;
    }
 
    // surface measures
    surface = Measures();
    // volume measures
    volume = Measures();
 
    // theoretical surface and volumes
    {
        double surfaceProd = 0.0;
        double volumeProd = 0.0;
        const double yMax = map->yMin + map->yStep * map->ySteps;
        for (unsigned i = 1; i < pcVertices->size(); ++i) {
            const CCVector3* P0 = pcVertices->getPoint(i - 1);
            const CCVector3* P1 = pcVertices->getPoint(i);
 
            // polyline: X = radius, Y = height
            double r0 = P0->x;
            double y0 = P0->y;
            double r1 = P1->x;
            double y1 = P1->y;
 
            // without loss of generality ;)
            if (y0 > y1) {
                std::swap(y0, y1);
                std::swap(r0, r1);
            }
 
            // segment is totally outside the map?
            if (y1 < map->yMin || y0 > yMax) {
                // we skip it
                continue;
            }
 
            if (y0 < map->yMin) {
                // interpolate r0 @ map->yMin
                double alpha = (map->yMin - y0) / (y1 - y0);
                assert(alpha >= 0.0 && alpha <= 1.0);
                r0 = r0 + alpha * (r1 - r0);
                y0 = map->yMin;
            } else if (y1 > yMax) {
                // interpolate r1 @ map->yMax
                double alpha = (yMax - y0) / (y1 - y0);
                assert(alpha >= 0.0 && alpha <= 1.0);
                r1 = r0 + alpha * (r1 - r0);
                y1 = yMax;
            }
 
            // product for truncated cone surface (see
            // http://en.wikipedia.org/wiki/Frustum)
            double segmentLength =
                    sqrt((r1 - r0) * (r1 - r0) + (y1 - y0) * (y1 - y0));
            surfaceProd += (r0 + r1) * segmentLength;
 
            // product for truncated cone volume (see
            // http://en.wikipedia.org/wiki/Frustum)
            volumeProd += (y1 - y0) * (r0 * r0 + r1 * r1 + r0 * r1);
        }
 
        surface.theoretical = M_PI * surfaceProd;
        volume.theoretical = M_PI / 3.0 * volumeProd;
    }
 
    int revolDim = GetPoylineRevolDim(profile);
    if (revolDim < 0) return false;
 
    // constant factors
    const double surfPart =
            map->xStep /
            2.0;  // perimeter of a portion of circle of angle alpha = alpha * r
                  // (* height to get the external surface)
    const double volPart = map->yStep * map->xStep /
                           6.0;  // area of a portion of circle of angle alpha =
                                 // alpha/2 * r^2 (* height to get the volume)
 
    const MapCell* cell = &map->at(0);
    // for each row
    for (unsigned j = 0; j < map->ySteps; ++j) {
        // corresponding heights
        double height1 = map->yMin + j * map->yStep;
        double height2 = height1 + map->yStep;
        double r_th1 = -1.0;
        double r_th2 = -1.0;
 
        // search nearest "segment" in polyline
        double height_middle = (height1 + height2) / 2.0;
        for (unsigned k = 1; k < vertexCount; ++k) {
            const CCVector3* A = vertices->getPoint(k - 1);
            const CCVector3* B = vertices->getPoint(k);
 
            double alpha = (height_middle - A->y) / (B->y - A->y);
            if (alpha >= 0.0 && alpha <= 1.0) {
                r_th1 = A->x + (height1 - A->y) / (B->y - A->y) * (B->x - A->x);
                r_th2 = A->x + (height2 - A->y) / (B->y - A->y) * (B->x - A->x);
                break;  // FIXME: we hope that there's only one segment facing
                        // this particular height?!
            }
        }
 
        if (r_th1 >= 0.0 /* && r_th2 >= 0.0*/) {
            // for each column
            for (unsigned i = 0; i < map->xSteps; ++i, ++cell) {
                // deviation from theory
                double d = (cell->count != 0 ? cell->value : 0.0);
 
                //"real" radius
                double r1 = r_th1 + d;  // see ComputeRadialDist --> << double
                                        // dist = radius - radius_th; >>
                double r2 = r_th2 +
                            d;  // FIXME: works only if the "scalar field" used
                                // for map creation was radial distances!!!
 
                // surface of the element (truncated cone external face)
                {
                    double s = sqrt((r2 - r1) * (r2 - r1) +
                                    map->yStep * map->yStep);
                    double externalSurface = /*surfPart * */ (r1 + r2) * s;
                    surface.total += externalSurface;
                    // dispatch in 'positive' and 'negative' surface
                    if (d >= 0.0)
                        surface.positive += externalSurface;
                    else
                        surface.negative += externalSurface;
                }
 
                // volume of the element
                {
                    volume.total +=
                            /*volPart * */ (r1 * r1 + r2 * r2 + r1 * r2);
                    // volume of the gain (or loss) of matter
                    double diffVolume = /*volPart * */ fabs(
                            3.0 * d *
                            (r_th1 + r_th2 +
                             d));  // = (r*r) * part - (r_th*r_th) * part =
                                   // [(r_th+d)*(r_th+d)-r_th*r_th] * part
                    if (d >= 0.0)
                        volume.positive += diffVolume;
                    else
                        volume.negative += diffVolume;
                }
            }
        } else {
            cell += map->xSteps;
        }
    }
 
    // don't forget to mult. by constants
    surface.total *= surfPart;
    surface.positive *= surfPart;
    surface.negative *= surfPart;
 
    volume.total *= volPart;
    volume.positive *= volPart;
    volume.negative *= volPart;
 
    return true;
}
 
bool DistanceMapGenerationTool::ConvertCloudToCylindrical(
        ccPointCloud* cloud,
        const ccGLMatrix&
                cloudToSurface,  // e.g. translation to the revolution origin
        unsigned char revolutionAxisDim,
        bool counterclockwise /*=false*/) {
    assert(cloud);
    if (!cloud || cloud->size() == 0) return false;
 
    // revolution axis
    const unsigned char Z = revolutionAxisDim;
    // we deduce the 2 other ('horizontal') dimensions from the revolution axis
    const unsigned char X = (Z < 2 ? Z + 1 : 0);
    const unsigned char Y = (X < 2 ? X + 1 : 0);
 
    // motion direction
    PointCoordinateType ccw = (counterclockwise ? -PC_ONE : PC_ONE);
 
    // get projection height
    for (unsigned n = 0; n < cloud->size(); ++n) {
        CCVector3* P = const_cast<CCVector3*>(cloud->getPoint(n));
        CCVector3 relativePos = cloudToSurface * (*P);
 
        // convert to cylindrical coordinates
        double lon_rad =
                ccw * atan2(relativePos.u[X], relativePos.u[Y]);  // longitude
        if (lon_rad < 0.0) {
            lon_rad += 2 * M_PI;
        }
 
        PointCoordinateType height = relativePos.u[Z];
 
        P->x = static_cast<PointCoordinateType>(lon_rad);
        P->y = height;
        P->z = 0;
    }
 
    cloud->refreshBB();
    if (cloud->getOctree()) {
        cloud->deleteOctree();
    }
    // TODO FIXME: and kd-trees? etc. We need a better way to handle those
    // cases...
 
    return true;
}
 
bool DistanceMapGenerationTool::ConvertCloudToConical(
        ccPointCloud* cloud,
        const ccGLMatrix&
                cloudToSurface,  // e.g. translation to the revolution origin
        unsigned char revolutionAxisDim,
        double latMin_rad,
        double latMax_rad,
        double conicalSpanRatio /*=1.0*/,
        bool counterclockwise /*=false*/) {
    assert(cloud);
    if (!cloud || cloud->size() == 0) return false;
 
    // revolution axis
    const unsigned char Z = revolutionAxisDim;
    // we deduce the 2 other ('horizontal') dimensions from the revolution axis
    const unsigned char X = (Z < 2 ? Z + 1 : 0);
    const unsigned char Y = (X < 2 ? X + 1 : 0);
 
    // motion direction
    PointCoordinateType ccw = (counterclockwise ? -PC_ONE : PC_ONE);
    // projection factor
    double nProj = ConicalProjectN(latMin_rad, latMax_rad) * conicalSpanRatio;
 
    // get projection height
    for (unsigned n = 0; n < cloud->size(); ++n) {
        CCVector3* P = const_cast<CCVector3*>(cloud->getPoint(n));
        CCVector3 relativePos = cloudToSurface * (*P);
 
        // convert to cylindrical coordinates
        PointCoordinateType ang_rad =
                ccw * atan2(relativePos.u[X], relativePos.u[Y]);
        if (ang_rad < 0.0)
            ang_rad += static_cast<PointCoordinateType>(2 * M_PI);
 
        double lat_rad =
                ComputeLatitude_rad(relativePos.u[X], relativePos.u[Y],
                                    relativePos.u[Z]);  // between 0 and pi/2
 
        *P = ProjectPointOnCone(ang_rad, lat_rad, latMin_rad, nProj,
                                counterclockwise);
    }
 
    cloud->refreshBB();
    if (cloud->getOctree()) {
        cloud->deleteOctree();
    }
    // TODO FIXME: and kd-trees? etc. We need a better way to handle those
    // cases...
 
    return true;
}
 
bool DistanceMapGenerationTool::SaveMapAsCSVMatrix(
        const QSharedPointer<Map>& map,
        QString filename,
        QString xUnit,
        QString yUnit,
        double xConversionFactor /*=1.0*/,
        double yConversionFactor /*=1.0*/,
        ecvMainAppInterface* app /*=0*/) {
    if (!map) {
        if (app)
            app->dispToConsole(QString("[SaveMapAsCSVMatrix] Internal error: "
                                       "invalid input map!"),
                               ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return false;
    }
 
    // write file
    QFile file(filename);
    if (!file.open(QFile::WriteOnly | QFile::Text)) {
        if (app)
            app->dispToConsole(QString("[SaveMapAsCSVMatrix] Failed to open "
                                       "file for writing! Check access rights"),
                               ecvMainAppInterface::ERR_CONSOLE_MESSAGE);
        return false;
    }
    QTextStream stream(&file);
 
    // write CSV header
    {
        // min and max height (for all lines)
        stream << QString("Height min (%1);").arg(yUnit);
        stream << QString("Height max (%1);").arg(yUnit);
 
        // for each column
        for (unsigned i = 0; i < map->xSteps; ++i) {
            // min and max angle for the current column
            double minX = xConversionFactor * (map->xMin + i * map->xStep);
            double maxX =
                    xConversionFactor * (map->xMin + (i + 1) * map->xStep);
            stream << QString("%1-%2 (%3);").arg(minX).arg(maxX).arg(xUnit);
        }
 
        // eol
        stream << QtCompat::endl;
    }
 
    // for each line
    for (unsigned j = 0; j < map->ySteps; ++j) {
        // min and max height (for the current line)
        double minY = yConversionFactor *
                      (map->yMin + (map->ySteps - 1 - j) * map->yStep);
        double maxY = yConversionFactor *
                      (map->yMin + (map->ySteps - j) * map->yStep);
        stream << QString::number(minY) << QString(";");
        stream << QString::number(maxY) << QString(";");
 
        // for each column
        for (unsigned i = 0; i < map->xSteps; ++i) {
            // write the grid value
            stream << QString::number(map->at(i + j * map->xSteps).value)
                   << QString(";");
        }
        // eol
        stream << QString("\n");
    }
 
    file.close();
 
    return true;
}
 
ccMesh* DistanceMapGenerationTool::ConvertProfileToMesh(
        ccPolyline* profile,
        const ccGLMatrix& cloudToProfile,
        bool counterclockwise,
        unsigned angularSteps /*=36*/,
        QImage mapTexture /*=QImage()*/) {
    if (!profile || angularSteps < 3) {
        return 0;
    }
 
    // profile vertices
    cloudViewer::GenericIndexedCloudPersist* profileVertices =
            profile->getAssociatedCloud();
    unsigned profVertCount = profileVertices->size();
    if (profVertCount < 2) {
        return 0;
    }
 
    // profile meta-data
    ProfileMetaData profileDesc;
    if (!GetPoylineMetaData(profile, profileDesc)) {
        assert(false);
        return 0;
    }
 
    unsigned char Z = static_cast<unsigned char>(profileDesc.revolDim);
    // we deduce the 2 other ('horizontal') dimensions
    const unsigned char X = (Z < 2 ? Z + 1 : 0);
    const unsigned char Y = (X < 2 ? X + 1 : 0);
 
    unsigned meshVertCount = profVertCount * angularSteps;
    unsigned meshFaceCount = (profVertCount - 1) * angularSteps * 2;
    ccPointCloud* cloud = new ccPointCloud("vertices");
    ccMesh* mesh = new ccMesh(cloud);
    if (!cloud->reserve(meshVertCount) || !mesh->reserve(meshFaceCount)) {
        // not enough memory
        delete cloud;
        delete mesh;
        return 0;
    }
 
    ccGLMatrix profileToCloud = cloudToProfile.inverse();
 
    // create vertices
    {
        double cwSign = (counterclockwise ? -1.0 : 1.0);
        for (unsigned j = 0; j < angularSteps; ++j) {
            double angle_rad =
                    static_cast<double>(j) / angularSteps * (2 * M_PI);
 
            CCVector3d N(sin(angle_rad) * cwSign, cos(angle_rad), 0);
 
            for (unsigned i = 0; i < profVertCount; ++i) {
                const CCVector3* P = profileVertices->getPoint(i);
                double radius = static_cast<double>(P->x);
 
                CCVector3 Pxyz;
                Pxyz.u[X] = static_cast<PointCoordinateType>(radius * N.x);
                Pxyz.u[Y] = static_cast<PointCoordinateType>(radius * N.y);
                Pxyz.u[Z] = P->y;
 
                profileToCloud.apply(Pxyz);
 
                cloud->addPoint(Pxyz);
            }
        }
        mesh->addChild(cloud);
    }
 
    PointCoordinateType h0 = profileVertices->getPoint(0)->y;
    PointCoordinateType dH =
            profileVertices->getPoint(profVertCount - 1)->y - h0;
    bool invertedHeight = (dH < 0);
 
    // create facets
    {
        for (unsigned j = 0; j < angularSteps; ++j) {
            unsigned nextJ = ((j + 1) % angularSteps);
            for (unsigned i = 0; i + 1 < profVertCount; ++i) {
                unsigned vertA = j * profVertCount + i;
                unsigned vertB = nextJ * profVertCount + i;
                unsigned vertC = vertB + 1;
                unsigned vertD = vertA + 1;
 
                if (invertedHeight) {
                    mesh->addTriangle(vertB, vertC, vertD);
                    mesh->addTriangle(vertB, vertD, vertA);
                } else {
                    mesh->addTriangle(vertB, vertD, vertC);
                    mesh->addTriangle(vertB, vertA, vertD);
                }
            }
        }
    }
 
    // do we have a texture as well?
    if (!mapTexture.isNull()) {
        // texture coordinates
        TextureCoordsContainer* texCoords = new TextureCoordsContainer();
        mesh->addChild(texCoords);
        if (!texCoords->reserveSafe(
                    meshVertCount +
                    profVertCount))  // we add a column for correct wrapping!
        {
            // not enough memory to finish the job!
            return mesh;
        }
 
        // create default texture coordinates
        for (unsigned j = 0; j <= angularSteps; ++j) {
            TexCoords2D T(static_cast<float>(j) / angularSteps, 0.0f);
            for (unsigned i = 0; i < profVertCount; ++i) {
                T.ty = (profileVertices->getPoint(i)->y - h0) / dH;
                if (invertedHeight) T.ty = 1.0f - T.ty;
                texCoords->addElement(T);
            }
        }
 
        if (!mesh->reservePerTriangleTexCoordIndexes()) {
            // not enough memory to finish the job!
            return mesh;
        }
 
        // set texture indexes
        {
            for (unsigned j = 0; j < angularSteps; ++j) {
                unsigned nextJ = ((j + 1) /*% angularSteps*/);
                for (unsigned i = 0; i + 1 < profVertCount; ++i) {
                    unsigned vertA = j * profVertCount + i;
                    unsigned vertB = nextJ * profVertCount + i;
                    unsigned vertC = vertB + 1;
                    unsigned vertD = vertA + 1;
 
                    if (invertedHeight) {
                        mesh->addTriangleTexCoordIndexes(vertB, vertC, vertD);
                        mesh->addTriangleTexCoordIndexes(vertB, vertD, vertA);
                    } else {
                        mesh->addTriangleTexCoordIndexes(vertB, vertD, vertC);
                        mesh->addTriangleTexCoordIndexes(vertB, vertA, vertD);
                    }
                }
            }
        }
 
        // set material indexes
        if (!mesh->reservePerTriangleMtlIndexes()) {
            // not enough memory to finish the job!
            mesh->removeChild(texCoords);
            mesh->removePerTriangleTexCoordIndexes();
            return mesh;
        }
        for (unsigned i = 0; i < meshFaceCount; ++i) {
            mesh->addTriangleMtlIndex(0);
        }
 
        // set material
#if 0
        {
            ccMaterial::Shared material(new ccMaterial("texture"));
            material->setTexture(mapTexture, QString(), false);
 
            ccMaterialSet* materialSet = new ccMaterialSet();
            materialSet->addMaterial(material);
 
            mesh->setMaterialSet(materialSet);
        }
#endif
 
        mesh->setTexCoordinatesTable(texCoords);
        mesh->showMaterials(true);
        mesh->setVisible(true);
        cloud->setVisible(false);
    }
 
    return mesh;
}
 
ccPointCloud* DistanceMapGenerationTool::ConvertMapToCloud(
        const QSharedPointer<Map>& map,
        ccPolyline* profile,
        double baseRadius /*=1.0*/,
        bool keepNaNPoints /*=true*/) {
    if (!map || !profile) return 0;
 
    unsigned count = map->ySteps * map->xSteps;
 
    ccPointCloud* cloud = new ccPointCloud("map");
    ccScalarField* sf = new ccScalarField("values");
    if (!cloud->reserve(count) || !sf->reserveSafe(count)) {
        // not enough memory
        delete cloud;
        sf->release();
        return nullptr;
    }
 
    // number of vertices
    cloudViewer::GenericIndexedCloudPersist* polyVertices =
            profile->getAssociatedCloud();
    unsigned polyVertCount = polyVertices->size();
    if (polyVertCount < 2) {
        delete cloud;
        sf->release();
        return nullptr;
    }
 
    // profile meta-data
    ProfileMetaData profileDesc;
    if (!GetPoylineMetaData(profile, profileDesc)) {
        delete cloud;
        sf->release();
        return nullptr;
    }
 
    unsigned char Z = static_cast<unsigned char>(profileDesc.revolDim);
    // we deduce the 2 other ('horizontal') dimensions
    const unsigned char X = (Z < 2 ? Z + 1 : 0);
    const unsigned char Y = (X < 2 ? X + 1 : 0);
 
    const double xStep =
            baseRadius * (2 * M_PI) / static_cast<double>(map->xSteps);
 
    const MapCell* cell = &map->at(0);
    for (unsigned j = 0; j < map->ySteps; ++j) {
        CCVector3 P(0, 0, 0);
        P.u[Z] = static_cast<PointCoordinateType>(map->yMin +
                                                  (j + 0.5) * map->yStep);
 
        // for each column
        for (unsigned i = 0; i < map->xSteps; ++i, ++cell) {
            if (keepNaNPoints || cell->count != 0) {
                P.u[X] = static_cast<PointCoordinateType>(map->xMin +
                                                          (i + 0.5) * xStep);
 
                // search nearest "segment" in polyline
                for (unsigned k = 1; k < polyVertCount; ++k) {
                    const CCVector3* A = polyVertices->getPoint(k - 1);
                    const CCVector3* B = polyVertices->getPoint(k);
 
                    double alpha = (P.u[Z] - profileDesc.heightShift - A->y) /
                                   (B->y - A->y);
                    if (alpha >= 0.0 && alpha <= 1.0) {
                        // we deduce the right radius by linear interpolation
                        double radius_th = A->x + alpha * (B->x - A->x);
                        // TODO: we take the first radius (even if there are
                        // other segments at this particular height, because we
                        // can't guess which one is the 'right' one!
                        P.u[Y] = static_cast<PointCoordinateType>(radius_th);
                        break;
                    }
                }
 
                cloud->addPoint(profileDesc.origin + P);
 
                ScalarType val = cell->count
                                         ? static_cast<ScalarType>(cell->value)
                                         : NAN_VALUE;
                sf->addElement(val);
            }
        }
    }
 
    sf->computeMinAndMax();
    int sfIdx = cloud->addScalarField(sf);
    cloud->setCurrentDisplayedScalarField(sfIdx);
    cloud->showSF(true);
    cloud->resize(cloud->size());  // if we have skipped NaN values!
 
    return cloud;
}
 
QImage DistanceMapGenerationTool::ConvertMapToImage(
        const QSharedPointer<Map>& map,
        ccColorScale::Shared colorScale,
        unsigned colorScaleSteps /*=ccColorScale::MAX_STEPS*/) {
    if (!map || !colorScale) return QImage();
 
    // create image
    QImage image(QSize(map->xSteps, map->ySteps), QImage::Format_ARGB32);
    if (image.isNull()) {
        // not enough memory!
        return QImage();
    }
 
    // convert map cells to pixels
    {
        bool csIsRelative = colorScale->isRelative();
 
        const MapCell* cell = &map->at(0);
        for (unsigned j = 0; j < map->ySteps; ++j) {
            // for each column
            for (unsigned i = 0; i < map->xSteps; ++i, ++cell) {
                const ecvColor::Rgb* rgb = &ecvColor::lightGrey;
 
                if (cell->count != 0) {
                    double relativePos =
                            csIsRelative ? (cell->value - map->minVal) /
                                                   (map->maxVal - map->minVal)
                                         : colorScale->getRelativePosition(
                                                   cell->value);
                    if (relativePos < 0.0)
                        relativePos = 0.0;
                    else if (relativePos > 1.0)
                        relativePos = 1.0;
                    rgb = colorScale->getColorByRelativePos(
                            relativePos, colorScaleSteps, &ecvColor::lightGrey);
                }
 
                // DGM FIXME: QImage::sePixel is quite slow!
                image.setPixel(static_cast<int>(i), static_cast<int>(j),
                               qRgb(rgb->r, rgb->g, rgb->b));
            }
        }
    }
 
    return image;
}