ilihelper.cpp

来自「支持各种栅格图像和矢量图像读取的库」· C++ 代码 · 共 111 行

/******************************************************************************
 * $Id: ilihelper.cpp 10646 2007-01-18 02:38:10Z warmerdam $
 *
 * Project:  Interlis 1/2 Translator
 * Purpose:  Helper functions for Interlis reader
 * Author:   Pirmin Kalberer, Sourcepole AG
 *
 ******************************************************************************
 * Copyright (c) 2004, Pirmin Kalberer, Sourcepole AG
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/


#include "ilihelper.h"


CPL_CVSID("$Id: ilihelper.cpp 10646 2007-01-18 02:38:10Z warmerdam $");


OGRPoint *getARCCenter(OGRPoint *ptStart, OGRPoint *ptArc, OGRPoint *ptEnd) {  
  // FIXME precision
  double bx = ptStart->getX(); double by = ptStart->getY();
  double cx = ptArc->getX(); double cy = ptArc->getY();
  double dx = ptEnd->getX(); double dy = ptEnd->getY();
  double temp, bc, cd, det, x, y;
  
  temp = cx*cx+cy*cy;
  bc = (bx*bx + by*by - temp)/2.0;
  cd = (temp - dx*dx - dy*dy)/2.0;
  det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
  
  OGRPoint *center = new OGRPoint();
  
  if (fabs(det) < 1.0e-6) { // could not determin the determinante: too small
    return center;
  }    
  det = 1/det;
  x = (bc*(cy-dy)-cd*(by-cy))*det;
  y = ((bx-cx)*cd-(cx-dx)*bc)*det;
  
  center->setX(x);
  center->setY(y);
  
  return center;
  //r = sqrt((x-bx)*(x-bx)+(y-by)*(y-by));
}

double getPhi(OGRPoint *center, OGRPoint *pt) {
  double cx = center->getX(); double cy = center->getY();
  double px = pt->getX(); double py = pt->getY();
  double r = sqrt((cx-px)*(cx-px)+(cy-py)*(cy-py));
  double phi = acos((px-cx)/r);
  return (py>cy) ? phi : -phi;
}

void interpolateArc(OGRLineString* line, OGRPoint *ptStart, OGRPoint *ptOnArc, OGRPoint *ptEnd, double arcIncr) {
  OGRPoint *center = getARCCenter(ptStart, ptOnArc, ptEnd);

  double cx = center->getX(); double cy = center->getY();
  double px = ptOnArc->getX(); double py = ptOnArc->getY();
  double r = sqrt((cx-px)*(cx-px)+(cy-py)*(cy-py));

  double phiPtStart = getPhi(center, ptStart);
  double phiPtOnArc = getPhi(center, ptOnArc);
  double phiPtEnd = getPhi(center, ptEnd);

  double deltaPhi = phiPtEnd - phiPtStart;
  if (deltaPhi < 0) deltaPhi += 2*PI;
  if (deltaPhi < PI) { //clockwise, if phiPtEnd - phiPtStart < PI
    if (phiPtEnd < phiPtStart) phiPtEnd += 2*PI;
    if (phiPtOnArc < phiPtStart) phiPtOnArc += 2*PI;
    for (double angle = phiPtStart; angle<phiPtEnd; angle += arcIncr) {
      if (angle != phiPtStart && angle != phiPtEnd)
        line->addPoint(center->getX()+r*cos(angle), center->getY()+r*sin(angle), 0);
      //add ptOnArc, if angle < phiPtOnArc < angle + arcIncr
      if (angle < phiPtOnArc && angle + arcIncr > phiPtOnArc)
        line->addPoint(ptOnArc);
    }
  } else {
    if (phiPtStart < phiPtEnd) phiPtStart += 2*PI;
    if (phiPtOnArc < phiPtEnd) phiPtOnArc += 2*PI;
    for (double angle = phiPtStart; angle>phiPtEnd; angle -= arcIncr) {
      if (angle != phiPtStart && angle != phiPtEnd)
        line->addPoint(center->getX()+r*cos(angle), center->getY()+r*sin(angle), 0);
      //add ptOnArc, if angle > phiPtOnArc > angle - arcIncr
      if (angle > phiPtOnArc && angle - arcIncr < phiPtOnArc)
        line->addPoint(ptOnArc);
    }
  }
  line->addPoint(ptEnd);
  delete center;
}