GeographicLib: Intersect.hpp Source File (original) (raw)

72 private:

74 public:

79 typedef std::pair<Math::real, Math::real> Point;

84 static const unsigned LineCaps = Geodesic::LATITUDE | Geodesic::LONGITUDE |

85 Geodesic::AZIMUTH | Geodesic::REDUCEDLENGTH | Geodesic::GEODESICSCALE |

86 Geodesic::DISTANCE_IN;

87 private:

88 static const int numit_ = 100;

90 real _a, _f,

91 _rR,

92 _d,

93 _eps,

94 _tol,

95 _delta,

96 _t1,

97 _t2,

98 _t3,

99 _t4,

100 _t5,

101 _d1,

102 _d2,

103 _d3;

104

106 { using std::fabs; return fabs(x) + fabs(y); }

107

108 class XPoint {

109 public:

111 int c;

113 : x(x), y(y), c(c)

114 {}

115 XPoint()

116 : x(Math::NaN()), y(Math::NaN()), c(0)

117 {}

118 XPoint(const Point& p)

119 : x(p.first), y(p.second), c(0)

120 {}

121 XPoint& operator+=(const XPoint& p) {

122 x += p.x; y += p.y;

123 if (p.c) c = p.c;

124 return *this;

125 }

126 XPoint operator+(const XPoint& p) const {

127 XPoint t = *this; t += p; return t;

128 }

129 Math::real Dist() const { return d1(x, y); }

130 Math::real Dist(const XPoint& p) const { return d1(x - p.x, y - p.y); }

131 Point data() const { return Point(x, y); }

132 };

133

134

136 private:

137 const real _delta;

138 public:

139 SetComp(Math::real delta) : _delta(delta) {}

140 bool eq(const XPoint& p, const XPoint& q) const {

141 return d1(p.x - q.x, p.y - q.y) <= _delta;

142 }

143 bool operator()(const XPoint& p, const XPoint& q) const {

144 return !eq(p, q) && ( p.x != q.x ? p.x < q.x : p.y < q.y );

145 }

146 };

147 SetComp _comp;

148

149 class RankPoint {

150 private:

151 const real _x, _y;

152 public:

153 RankPoint(const Point& p0) : _x(p0.first), _y(p0.second) {}

154 RankPoint(const XPoint& p0) : _x(p0.x), _y(p0.y) {}

155 bool operator()(const XPoint& p, const XPoint& q) const {

156 real dp = d1(p.x - _x, p.y - _y),

157 dq = d1(q.x - _x, q.y - _y);

158 return dp != dq ? (dp < dq) :

159 (p.x != q.x ? (p.x < q.x) : (p.y < q.y));

160 }

161 };

162

163 XPoint Spherical(const GeodesicLine& lineX, const GeodesicLine& lineY,

164 const XPoint& p) const;

165

166 XPoint Basic(const GeodesicLine& lineX, const GeodesicLine& lineY,

167 const XPoint& p0) const;

168

169 XPoint ClosestInt(const GeodesicLine& lineX, const GeodesicLine& lineY,

170 const XPoint& p0) const;

171

172 XPoint NextInt(const GeodesicLine& lineX, const GeodesicLine& lineY) const;

173

174 XPoint SegmentInt(const GeodesicLine& lineX, const GeodesicLine& lineY,

175 int& segmode) const;

176

177 std::vector

178 AllInt0(const GeodesicLine& lineX, const GeodesicLine& lineY,

179 Math::real maxdist, const XPoint& p0) const;

180 std::vector

181 AllInternal(const GeodesicLine& lineX, const GeodesicLine& lineY,

182 Math::real maxdist, const Point& p0,

183 std::vector& c, bool cp) const;

184

185

186 Math::real ConjugateDist(const GeodesicLine& line, Math::real s3, bool semi,

187 Math::real m12 = 0, Math::real M12 = 1,

188 Math::real M21 = 1) const;

189 Math::real distpolar(Math::real lat1, Math::real* lat2 = nullptr) const;

190 Math::real polarb(Math::real* lata = nullptr, Math::real* latb = nullptr)

191 const;

192 Math::real conjdist(Math::real azi, Math::real* ds = nullptr,

193 Math::real* sp = nullptr, Math::real* sm = nullptr)

194 const;

195 Math::real distoblique(Math::real* azi = nullptr, Math::real* sp = nullptr,

196 Math::real* sm = nullptr) const;

197

198

199

200 static XPoint fixcoincident(const XPoint& p0, const XPoint& p);

201 static XPoint fixcoincident(const XPoint& p0, const XPoint& p, int c);

202 static XPoint fixsegment(Math::real sx, Math::real sy, const XPoint& p);

203 static int segmentmode(Math::real sx, Math::real sy, const XPoint& p) {

204 return (p.x < 0 ? -1 : p.x <= sx ? 0 : 1) * 3

205 + (p.y < 0 ? -1 : p.y <= sy ? 0 : 1);

206 }

207 mutable long long _cnt0, _cnt1, _cnt2, _cnt3, _cnt4;

208 public:

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228 Intersect(const Geodesic& geod);

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251 Point Closest(Math::real latX, Math::real lonX, Math::real aziX,

252 Math::real latY, Math::real lonY, Math::real aziY,

253 const Point& p0 = Point(0, 0), int* c = nullptr) const;

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271 Point Closest(const GeodesicLine& lineX, const GeodesicLine& lineY,

272 const Point& p0 = Point(0, 0), int* c = nullptr) const;

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307 Point Segment(Math::real latX1, Math::real lonX1,

308 Math::real latX2, Math::real lonX2,

309 Math::real latY1, Math::real lonY1,

310 Math::real latY2, Math::real lonY2,

311 int& segmode, int* c = nullptr) const;

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337 Point Segment(const GeodesicLine& lineX, const GeodesicLine& lineY,

338 int& segmode, int* c = nullptr) const;

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359 Point Next(Math::real latX, Math::real lonX,

360 Math::real aziX, Math::real aziY, int* c = nullptr) const;

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383 Point Next(const GeodesicLine& lineX, const GeodesicLine& lineY,

384 int* c = nullptr) const;

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411 std::vector All(Math::real latX, Math::real lonX, Math::real aziX,

412 Math::real latY, Math::real lonY, Math::real aziY,

413 Math::real maxdist, std::vector& c,

414 const Point& p0 = Point(0, 0))

415 const;

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437 std::vector All(Math::real latX, Math::real lonX, Math::real aziX,

438 Math::real latY, Math::real lonY, Math::real aziY,

439 Math::real maxdist, const Point& p0 = Point(0, 0))

440 const;

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462 std::vector All(const GeodesicLine& lineX, const GeodesicLine& lineY,

463 Math::real maxdist, std::vector& c,

464 const Point& p0 = Point(0, 0))

465 const;

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486 std::vector All(const GeodesicLine& lineX, const GeodesicLine& lineY,

487 Math::real maxdist, const Point& p0 = Point(0, 0))

488 const;

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

505

506

507

508

509

510

511

512

513

514 long long NumBasic() const { return _cnt1; }

515

516

517

518

519

520

521

522

523

524

525

526

527

528 long long NumChange() const { return _cnt2; }

529

530

531

532

533

534

535

536

537 long long NumCorner() const { return _cnt3; }

538

539

540

541

542

543

544

545

546

547

548

549

551

552

553

554

555

556

557

558

559

560

561

562

564

565

566

567

568

569

570

571

572

573

574

575

577 using std::fabs;

578 return fabs(p.first - p0.first) + fabs(p.second - p0.second);

579 }

580 };