Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Microstation DGN Access Library
       4             :  * Purpose:  Code to stroke Arcs/Ellipses into polylines.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2001, Avenza Systems Inc, http://www.avenza.com/
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "dgnlibp.h"
      14             : #include <cmath>
      15             : 
      16             : constexpr double DEG_TO_RAD = M_PI / 180.0;
      17             : 
      18             : /************************************************************************/
      19             : /*                         ComputePointOnArc()                          */
      20             : /************************************************************************/
      21             : 
      22         438 : static void ComputePointOnArc2D(double dfPrimary, double dfSecondary,
      23             :                                 double dfAxisRotation, double dfAngle,
      24             :                                 double *pdfX, double *pdfY)
      25             : 
      26             : {
      27             :     // dfAxisRotation and dfAngle are supposed to be in Radians
      28         438 :     const double dfCosRotation = cos(dfAxisRotation);
      29         438 :     const double dfSinRotation = sin(dfAxisRotation);
      30         438 :     const double dfEllipseX = dfPrimary * cos(dfAngle);
      31         438 :     const double dfEllipseY = dfSecondary * sin(dfAngle);
      32             : 
      33         438 :     *pdfX = dfEllipseX * dfCosRotation - dfEllipseY * dfSinRotation;
      34         438 :     *pdfY = dfEllipseX * dfSinRotation + dfEllipseY * dfCosRotation;
      35         438 : }
      36             : 
      37             : /************************************************************************/
      38             : /*                            DGNStrokeArc()                            */
      39             : /************************************************************************/
      40             : 
      41             : /**
      42             :  * Generate a polyline approximation of an arc.
      43             :  *
      44             :  * Produce a series of equidistant (actually equi-angle) points along
      45             :  * an arc.  Currently this only works for 2D arcs (and ellipses).
      46             :  *
      47             :  * @param hFile the DGN file to which the arc belongs (currently not used).
      48             :  * @param psArc the arc to be approximated.
      49             :  * @param nPoints the number of points to use to approximate the arc.
      50             :  * @param pasPoints the array of points into which to put the results.
      51             :  * There must be room for at least nPoints points.
      52             :  *
      53             :  * @return TRUE on success or FALSE on failure.
      54             :  */
      55             : 
      56           6 : int DGNStrokeArc(CPL_UNUSED DGNHandle hFile, DGNElemArc *psArc, int nPoints,
      57             :                  DGNPoint *pasPoints)
      58             : {
      59           6 :     if (nPoints < 2)
      60           0 :         return FALSE;
      61             : 
      62           6 :     if (psArc->primary_axis == 0.0 || psArc->secondary_axis == 0.0)
      63             :     {
      64           0 :         CPLError(CE_Warning, CPLE_AppDefined,
      65             :                  "Zero primary or secondary axis in DGNStrokeArc().");
      66           0 :         return FALSE;
      67             :     }
      68             : 
      69           6 :     const double dfAngleStep = psArc->sweepang / (nPoints - 1);
      70         444 :     for (int i = 0; i < nPoints; i++)
      71             :     {
      72         438 :         const double dfAngle = (psArc->startang + dfAngleStep * i) * DEG_TO_RAD;
      73             : 
      74         438 :         ComputePointOnArc2D(psArc->primary_axis, psArc->secondary_axis,
      75         438 :                             psArc->rotation * DEG_TO_RAD, dfAngle,
      76         438 :                             &(pasPoints[i].x), &(pasPoints[i].y));
      77         438 :         pasPoints[i].x += psArc->origin.x;
      78         438 :         pasPoints[i].y += psArc->origin.y;
      79         438 :         pasPoints[i].z = psArc->origin.z;
      80             :     }
      81             : 
      82           6 :     return TRUE;
      83             : }
      84             : 
      85             : /************************************************************************/
      86             : /*                           DGNStrokeCurve()                           */
      87             : /************************************************************************/
      88             : 
      89             : /**
      90             :  * Generate a polyline approximation of an curve.
      91             :  *
      92             :  * Produce a series of equidistant points along a microstation curve element.
      93             :  * Currently this only works for 2D.
      94             :  *
      95             :  * @param hFile the DGN file to which the arc belongs (currently not used).
      96             :  * @param psCurve the curve to be approximated.
      97             :  * @param nPoints the number of points to use to approximate the curve.
      98             :  * @param pasPoints the array of points into which to put the results.
      99             :  * There must be room for at least nPoints points.
     100             :  *
     101             :  * @return TRUE on success or FALSE on failure.
     102             :  */
     103             : 
     104           0 : int DGNStrokeCurve(CPL_UNUSED DGNHandle hFile, DGNElemMultiPoint *psCurve,
     105             :                    int nPoints, DGNPoint *pasPoints)
     106             : {
     107           0 :     const int nDGNPoints = psCurve->num_vertices;
     108             : 
     109           0 :     if (nDGNPoints < 6)
     110           0 :         return FALSE;
     111             : 
     112           0 :     if (nPoints < nDGNPoints - 4)
     113           0 :         return FALSE;
     114             : 
     115             :     /* -------------------------------------------------------------------- */
     116             :     /*      Compute the Compute the slopes/distances of the segments.       */
     117             :     /* -------------------------------------------------------------------- */
     118             :     double *padfMx =
     119           0 :         static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
     120             :     double *padfMy =
     121           0 :         static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
     122             :     double *padfD =
     123           0 :         static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
     124             :     double *padfTx =
     125           0 :         static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
     126             :     double *padfTy =
     127           0 :         static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
     128             : 
     129           0 :     double dfTotalD = 0.0;
     130             : 
     131           0 :     DGNPoint *pasDGNPoints = psCurve->vertices;
     132             : 
     133           0 :     for (int k = 0; k < nDGNPoints - 1; k++)
     134             :     {
     135             :         /* coverity[overrun-local] */
     136           0 :         padfD[k] = sqrt((pasDGNPoints[k + 1].x - pasDGNPoints[k].x) *
     137           0 :                             (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) +
     138           0 :                         (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) *
     139           0 :                             (pasDGNPoints[k + 1].y - pasDGNPoints[k].y));
     140           0 :         if (padfD[k] == 0.0)
     141             :         {
     142           0 :             padfD[k] = 0.0001;
     143           0 :             padfMx[k] = 0.0;
     144           0 :             padfMy[k] = 0.0;
     145             :         }
     146             :         else
     147             :         {
     148           0 :             padfMx[k] = (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k];
     149           0 :             padfMy[k] = (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k];
     150             :         }
     151             : 
     152           0 :         if (k > 1 && k < nDGNPoints - 3)
     153           0 :             dfTotalD += padfD[k];
     154             :     }
     155             : 
     156             :     /* -------------------------------------------------------------------- */
     157             :     /*      Compute the Tx, and Ty coefficients for each segment.           */
     158             :     /* -------------------------------------------------------------------- */
     159           0 :     for (int k = 2; k < nDGNPoints - 2; k++)
     160             :     {
     161           0 :         if (fabs(padfMx[k + 1] - padfMx[k]) == 0.0 &&
     162           0 :             fabs(padfMx[k - 1] - padfMx[k - 2]) == 0.0)
     163             :         {
     164           0 :             padfTx[k] = (padfMx[k] + padfMx[k - 1]) / 2;
     165             :         }
     166             :         else
     167             :         {
     168           0 :             padfTx[k] = (padfMx[k - 1] * fabs(padfMx[k + 1] - padfMx[k]) +
     169           0 :                          padfMx[k] * fabs(padfMx[k - 1] - padfMx[k - 2])) /
     170           0 :                         (std::abs(padfMx[k + 1] - padfMx[k]) +
     171           0 :                          std::abs(padfMx[k - 1] - padfMx[k - 2]));
     172             :         }
     173             : 
     174           0 :         if (fabs(padfMy[k + 1] - padfMy[k]) == 0.0 &&
     175           0 :             fabs(padfMy[k - 1] - padfMy[k - 2]) == 0.0)
     176             :         {
     177           0 :             padfTy[k] = (padfMy[k] + padfMy[k - 1]) / 2;
     178             :         }
     179             :         else
     180             :         {
     181           0 :             padfTy[k] = (padfMy[k - 1] * fabs(padfMy[k + 1] - padfMy[k]) +
     182           0 :                          padfMy[k] * fabs(padfMy[k - 1] - padfMy[k - 2])) /
     183           0 :                         (std::abs(padfMy[k + 1] - padfMy[k]) +
     184           0 :                          std::abs(padfMy[k - 1] - padfMy[k - 2]));
     185             :         }
     186             :     }
     187             : 
     188             :     /* -------------------------------------------------------------------- */
     189             :     /*      Determine a step size in D.  We scale things so that we have    */
     190             :     /*      roughly equidistant steps in D, but assume we also want to      */
     191             :     /*      include every node along the way.                               */
     192             :     /* -------------------------------------------------------------------- */
     193           0 :     double dfStepSize = dfTotalD / (nPoints - (nDGNPoints - 4) - 1);
     194             : 
     195             :     /* ==================================================================== */
     196             :     /*      Process each of the segments.                                   */
     197             :     /* ==================================================================== */
     198           0 :     double dfD = dfStepSize;
     199           0 :     int iOutPoint = 0;
     200             : 
     201           0 :     for (int k = 2; k < nDGNPoints - 3; k++)
     202             :     {
     203             :         /* --------------------------------------------------------------------
     204             :          */
     205             :         /*      Compute the "x" coefficients for this segment. */
     206             :         /* --------------------------------------------------------------------
     207             :          */
     208           0 :         const double dfCx = padfTx[k];
     209           0 :         const double dfBx =
     210           0 :             (3.0 * (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k] -
     211           0 :              2.0 * padfTx[k] - padfTx[k + 1]) /
     212           0 :             padfD[k];
     213           0 :         const double dfAx =
     214           0 :             (padfTx[k] + padfTx[k + 1] -
     215           0 :              2 * (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k]) /
     216           0 :             (padfD[k] * padfD[k]);
     217             : 
     218             :         /* --------------------------------------------------------------------
     219             :          */
     220             :         /*      Compute the Y coefficients for this segment. */
     221             :         /* --------------------------------------------------------------------
     222             :          */
     223           0 :         const double dfCy = padfTy[k];
     224           0 :         const double dfBy =
     225           0 :             (3.0 * (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k] -
     226           0 :              2.0 * padfTy[k] - padfTy[k + 1]) /
     227           0 :             padfD[k];
     228           0 :         const double dfAy =
     229           0 :             (padfTy[k] + padfTy[k + 1] -
     230           0 :              2 * (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k]) /
     231           0 :             (padfD[k] * padfD[k]);
     232             : 
     233             :         /* --------------------------------------------------------------------
     234             :          */
     235             :         /*      Add the start point for this segment. */
     236             :         /* --------------------------------------------------------------------
     237             :          */
     238           0 :         pasPoints[iOutPoint].x = pasDGNPoints[k].x;
     239           0 :         pasPoints[iOutPoint].y = pasDGNPoints[k].y;
     240           0 :         pasPoints[iOutPoint].z = 0.0;
     241           0 :         iOutPoint++;
     242             : 
     243             :         /* --------------------------------------------------------------------
     244             :          */
     245             :         /*      Step along, adding intermediate points. */
     246             :         /* --------------------------------------------------------------------
     247             :          */
     248           0 :         while (dfD < padfD[k] && iOutPoint < nPoints - (nDGNPoints - k - 1))
     249             :         {
     250           0 :             pasPoints[iOutPoint].x = dfAx * dfD * dfD * dfD + dfBx * dfD * dfD +
     251           0 :                                      dfCx * dfD + pasDGNPoints[k].x;
     252           0 :             pasPoints[iOutPoint].y = dfAy * dfD * dfD * dfD + dfBy * dfD * dfD +
     253           0 :                                      dfCy * dfD + pasDGNPoints[k].y;
     254           0 :             pasPoints[iOutPoint].z = 0.0;
     255           0 :             iOutPoint++;
     256             : 
     257           0 :             dfD += dfStepSize;
     258             :         }
     259             : 
     260           0 :         dfD -= padfD[k];
     261             :     }
     262             : 
     263             :     /* -------------------------------------------------------------------- */
     264             :     /*      Add the start point for this segment.                           */
     265             :     /* -------------------------------------------------------------------- */
     266           0 :     while (iOutPoint < nPoints)
     267             :     {
     268           0 :         pasPoints[iOutPoint].x = pasDGNPoints[nDGNPoints - 3].x;
     269           0 :         pasPoints[iOutPoint].y = pasDGNPoints[nDGNPoints - 3].y;
     270           0 :         pasPoints[iOutPoint].z = 0.0;
     271           0 :         iOutPoint++;
     272             :     }
     273             : 
     274             :     /* -------------------------------------------------------------------- */
     275             :     /*      Cleanup.                                                        */
     276             :     /* -------------------------------------------------------------------- */
     277           0 :     CPLFree(padfMx);
     278           0 :     CPLFree(padfMy);
     279           0 :     CPLFree(padfD);
     280           0 :     CPLFree(padfTx);
     281           0 :     CPLFree(padfTy);
     282             : 
     283           0 :     return TRUE;
     284             : }
     285             : 
     286             : /************************************************************************/
     287             : /*                                main()                                */
     288             : /*                                                                      */
     289             : /*      test mainline                                                   */
     290             : /************************************************************************/
     291             : #ifdef notdef
     292             : int main(int argc, char **argv)
     293             : 
     294             : {
     295             :     if (argc != 5)
     296             :     {
     297             :         printf(  // ok
     298             :             "Usage: stroke primary_axis secondary_axis axis_rotation angle\n");
     299             :         exit(1);
     300             :     }
     301             : 
     302             :     const double dfPrimary = CPLAtof(argv[1]);
     303             :     const double dfSecondary = CPLAtof(argv[2]);
     304             :     const double dfAxisRotation = CPLAtof(argv[3]) / 180 * M_PI;
     305             :     const double dfAngle = CPLAtof(argv[4]) / 180 * M_PI;
     306             : 
     307             :     double dfX = 0.0;
     308             :     double dfY = 0.0;
     309             :     ComputePointOnArc2D(dfPrimary, dfSecondary, dfAxisRotation, dfAngle, &dfX,
     310             :                         &dfY);
     311             : 
     312             :     printf("X=%.2f, Y=%.2f\n", dfX, dfY);  // ok
     313             : 
     314             :     return 0;
     315             : }
     316             : 
     317             : #endif