Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  NTF Translator
       4             :  * Purpose:  NTF Arc to polyline stroking code.  This code is really generic,
       5             :  *           and might be moved into an OGR module at some point in the
       6             :  *           future.
       7             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       8             :  *
       9             :  ******************************************************************************
      10             :  * Copyright (c) 2001, Frank Warmerdam
      11             :  *
      12             :  * SPDX-License-Identifier: MIT
      13             :  ****************************************************************************/
      14             : 
      15             : #include <stdarg.h>
      16             : #include "ntf.h"
      17             : #include "cpl_conv.h"
      18             : #include "cpl_string.h"
      19             : 
      20             : #include <algorithm>
      21             : #include <utility>
      22             : 
      23             : /************************************************************************/
      24             : /*                     NTFArcCenterFromEdgePoints()                     */
      25             : /*                                                                      */
      26             : /*      Compute the center of an arc/circle from three edge points.     */
      27             : /************************************************************************/
      28             : 
      29           0 : int NTFArcCenterFromEdgePoints(double x_c0, double y_c0, double x_c1,
      30             :                                double y_c1, double x_c2, double y_c2,
      31             :                                double *x_center, double *y_center)
      32             : 
      33             : {
      34             : 
      35             :     /* -------------------------------------------------------------------- */
      36             :     /*      Handle a degenerate case that occurs in OSNI products by        */
      37             :     /*      making some assumptions.  If the first and third points are     */
      38             :     /*      the same assume they are intended to define a full circle,      */
      39             :     /*      and that the second point is on the opposite side of the        */
      40             :     /*      circle.                                                         */
      41             :     /* -------------------------------------------------------------------- */
      42           0 :     if (x_c0 == x_c2 && y_c0 == y_c2)
      43             :     {
      44           0 :         *x_center = (x_c0 + x_c1) * 0.5;
      45           0 :         *y_center = (y_c0 + y_c1) * 0.5;
      46             : 
      47           0 :         return TRUE;
      48             :     }
      49             : 
      50             :     /* -------------------------------------------------------------------- */
      51             :     /*      Compute the inverse of the slopes connecting the first and      */
      52             :     /*      second points.  Also compute the center point of the two        */
      53             :     /*      lines ... the point our crossing line will go through.          */
      54             :     /* -------------------------------------------------------------------- */
      55           0 :     const double m1 =
      56           0 :         (y_c1 - y_c0) != 0.0 ? (x_c0 - x_c1) / (y_c1 - y_c0) : 1e+10;
      57             : 
      58           0 :     const double x1 = (x_c0 + x_c1) * 0.5;
      59           0 :     const double y1 = (y_c0 + y_c1) * 0.5;
      60             : 
      61             :     /* -------------------------------------------------------------------- */
      62             :     /*      Compute the same for the second point compared to the third     */
      63             :     /*      point.                                                          */
      64             :     /* -------------------------------------------------------------------- */
      65           0 :     const double m2 =
      66           0 :         (y_c2 - y_c1) != 0.0 ? (x_c1 - x_c2) / (y_c2 - y_c1) : 1e+10;
      67             : 
      68           0 :     const double x2 = (x_c1 + x_c2) * 0.5;
      69           0 :     const double y2 = (y_c1 + y_c2) * 0.5;
      70             : 
      71             :     /* -------------------------------------------------------------------- */
      72             :     /*      Turn these into the Ax+By+C = 0 form of the lines.              */
      73             :     /* -------------------------------------------------------------------- */
      74           0 :     const double a1 = m1;
      75           0 :     const double a2 = m2;
      76             : 
      77           0 :     const double b1 = -1.0;
      78           0 :     const double b2 = -1.0;
      79             : 
      80           0 :     const double c1 = (y1 - m1 * x1);
      81           0 :     const double c2 = (y2 - m2 * x2);
      82             : 
      83             :     /* -------------------------------------------------------------------- */
      84             :     /*      Compute the intersection of the two lines through the center    */
      85             :     /*      of the circle, using Kramers rule.                              */
      86             :     /* -------------------------------------------------------------------- */
      87           0 :     if (a1 * b2 - a2 * b1 == 0.0)
      88           0 :         return FALSE;
      89             : 
      90           0 :     const double det_inv = 1 / (a1 * b2 - a2 * b1);
      91             : 
      92           0 :     *x_center = (b1 * c2 - b2 * c1) * det_inv;
      93           0 :     *y_center = (a2 * c1 - a1 * c2) * det_inv;
      94             : 
      95           0 :     return TRUE;
      96             : }
      97             : 
      98             : /************************************************************************/
      99             : /*                  NTFStrokeArcToOGRGeometry_Points()                  */
     100             : /************************************************************************/
     101             : 
     102           0 : OGRGeometry *NTFStrokeArcToOGRGeometry_Points(double dfStartX, double dfStartY,
     103             :                                               double dfAlongX, double dfAlongY,
     104             :                                               double dfEndX, double dfEndY,
     105             :                                               int nVertexCount)
     106             : 
     107             : {
     108           0 :     double dfStartAngle = 0.0;
     109           0 :     double dfEndAngle = 0.0;
     110           0 :     double dfCenterX = 0.0;
     111           0 :     double dfCenterY = 0.0;
     112           0 :     double dfRadius = 0.0;
     113             : 
     114           0 :     if (!NTFArcCenterFromEdgePoints(dfStartX, dfStartY, dfAlongX, dfAlongY,
     115             :                                     dfEndX, dfEndY, &dfCenterX, &dfCenterY))
     116           0 :         return nullptr;
     117             : 
     118           0 :     if (dfStartX == dfEndX && dfStartY == dfEndY)
     119             :     {
     120           0 :         dfStartAngle = 0.0;
     121           0 :         dfEndAngle = 360.0;
     122             :     }
     123             :     else
     124             :     {
     125           0 :         double dfDeltaX = dfStartX - dfCenterX;
     126           0 :         double dfDeltaY = dfStartY - dfCenterY;
     127           0 :         dfStartAngle = atan2(dfDeltaY, dfDeltaX) * 180.0 / M_PI;
     128             : 
     129           0 :         dfDeltaX = dfAlongX - dfCenterX;
     130           0 :         dfDeltaY = dfAlongY - dfCenterY;
     131           0 :         double dfAlongAngle = atan2(dfDeltaY, dfDeltaX) * 180.0 / M_PI;
     132             : 
     133           0 :         dfDeltaX = dfEndX - dfCenterX;
     134           0 :         dfDeltaY = dfEndY - dfCenterY;
     135           0 :         dfEndAngle = atan2(dfDeltaY, dfDeltaX) * 180.0 / M_PI;
     136             : 
     137           0 :         while (dfAlongAngle < dfStartAngle)
     138           0 :             dfAlongAngle += 360.0;
     139             : 
     140           0 :         while (dfEndAngle < dfAlongAngle)
     141           0 :             dfEndAngle += 360.0;
     142             : 
     143           0 :         if (dfEndAngle - dfStartAngle > 360.0)
     144             :         {
     145           0 :             std::swap(dfStartAngle, dfEndAngle);
     146             : 
     147           0 :             while (dfEndAngle < dfStartAngle)
     148           0 :                 dfStartAngle -= 360.0;
     149             :         }
     150             :     }
     151             : 
     152           0 :     dfRadius = sqrt((dfCenterX - dfStartX) * (dfCenterX - dfStartX) +
     153           0 :                     (dfCenterY - dfStartY) * (dfCenterY - dfStartY));
     154             : 
     155           0 :     return NTFStrokeArcToOGRGeometry_Angles(
     156           0 :         dfCenterX, dfCenterY, dfRadius, dfStartAngle, dfEndAngle, nVertexCount);
     157             : }
     158             : 
     159             : /************************************************************************/
     160             : /*                  NTFStrokeArcToOGRGeometry_Angles()                  */
     161             : /************************************************************************/
     162             : 
     163           0 : OGRGeometry *NTFStrokeArcToOGRGeometry_Angles(double dfCenterX,
     164             :                                               double dfCenterY, double dfRadius,
     165             :                                               double dfStartAngle,
     166             :                                               double dfEndAngle,
     167             :                                               int nVertexCount)
     168             : 
     169             : {
     170           0 :     OGRLineString *poLine = new OGRLineString;
     171             : 
     172           0 :     nVertexCount = std::max(2, nVertexCount);
     173           0 :     const double dfSlice = (dfEndAngle - dfStartAngle) / (nVertexCount - 1);
     174             : 
     175           0 :     poLine->setNumPoints(nVertexCount);
     176             : 
     177           0 :     for (int iPoint = 0; iPoint < nVertexCount; iPoint++)
     178             :     {
     179           0 :         const double dfAngle = (dfStartAngle + iPoint * dfSlice) * M_PI / 180.0;
     180             : 
     181           0 :         const double dfArcX = dfCenterX + cos(dfAngle) * dfRadius;
     182           0 :         const double dfArcY = dfCenterY + sin(dfAngle) * dfRadius;
     183             : 
     184           0 :         poLine->setPoint(iPoint, dfArcX, dfArcY);
     185             :     }
     186             : 
     187           0 :     return poLine;
     188             : }