~ubuntu-branches/ubuntu/trusty/qgis/trusty : revision 18

src/core/pal/pointset.cpp

* libpal - Automated Placement of Labels Library

* University of Applied Sciences, Western Switzerland

* http://www.hes-so.ch

* Contact:

* maxence.laurent <at> heig-vd <dot> ch

* or

* eric.taillard <at> heig-vd <dot> ch

* This file is part of libpal.

* libpal is free software: you can redistribute it and/or modify

* it under the terms of the GNU General Public License as published by

* the Free Software Foundation, either version 3 of the License, or

* (at your option) any later version.

* libpal is distributed in the hope that it will be useful,

* but WITHOUT ANY WARRANTY; without even the implied warranty of

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

* GNU General Public License for more details.

* You should have received a copy of the GNU General Public License

* along with libpal. If not, see <http://www.gnu.org/licenses/>.

#ifndef _POINT_SET_H_

#define _POINT_SET_H_

#if defined(_VERBOSE_) || (_DEBUG_) || (_DEBUG_FULL_)

#include <iostream>

#endif

#ifdef HAVE_CONFIG_H

#include <config.h>

#endif

#include "pointset.h"

#include "util.h"

#include <pal/pal.h>

#include "geomfunction.h"

namespace pal

{

PointSet::PointSet()

{

nbPoints = cHullSize = 0;

x = NULL;

y = NULL;

cHull = NULL;

type = -1;

}

PointSet::PointSet( int nbPoints, double *x, double *y )

{

this->nbPoints = nbPoints;

this->x = new double[nbPoints];

this->y = new double[nbPoints];

int i;

for ( i = 0; i < nbPoints; i++ )

{

this->x[i] = x[i];

this->y[i] = y[i];

}

type = GEOS_POLYGON;

cHull = NULL;

}

PointSet::PointSet( double x, double y )

{

nbPoints = cHullSize = 1;

this->x = new double[1];

this->y = new double[1];

this->x[0] = x;

this->y[0] = y;

cHull = NULL;

parent = NULL;

holeOf = NULL;

type = GEOS_POINT;

}

PointSet::PointSet( PointSet &ps )

{

int i;

nbPoints = ps.nbPoints;

x = new double[nbPoints];

y = new double[nbPoints];

100

101

102

for ( i = 0; i < nbPoints; i++ )

103

{

104

x[i] = ps.x[i];

105

y[i] = ps.y[i];

106

}

107

108

if ( ps.cHull )

109

{

110

cHullSize = ps.cHullSize;

111

for ( i = 0; i < cHullSize; i++ )

112

{

113

cHull[i] = ps.cHull[i];

114

}

115

}

116

else

117

{

118

cHull = NULL;

119

cHullSize = 0;

120

}

121

122

type = ps.type;

123

124

holeOf = ps.holeOf;

125

}

126

127

PointSet::~PointSet()

128

{

129

deleteCoords();

130

131

if ( cHull )

132

delete[] cHull;

133

}

134

135

void PointSet::deleteCoords()

136

{

137

if ( x )

138

delete[] x;

139

if ( y )

140

delete[] y;

141

x = NULL;

142

y = NULL;

143

}

144

145

146

147

148

PointSet* PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )

149

{

150

151

int i, j;

152

153

PointSet *newShape = new PointSet();

154

155

newShape->type = GEOS_POLYGON;

156

157

newShape->nbPoints = nbPtSh;

158

159

newShape->x = new double[newShape->nbPoints];

160

newShape->y = new double[newShape->nbPoints];

161

#ifdef _DEBUG_FULL_

162

std::cout << "New shape: ";

163

#endif

164

// new shape # 1 from imin to imax

165

for ( j = 0, i = imin; i != ( imax + 1 ) % nbPoints; i = ( i + 1 ) % nbPoints, j++ )

166

{

167

newShape->x[j] = x[i];

168

newShape->y[j] = y[i];

169

#ifdef _DEBUG_FULL_

170

std::cout << x[i] << ";" << y[i] << std::endl;

171

#endif

172

}

173

// is the cutting point a new one ?

174

if ( fps != fpe )

175

{

176

// yes => so add it

177

newShape->x[j] = fptx;

178

newShape->y[j] = fpty;

179

#ifdef _DEBUG_FULL_

180

std::cout << fptx << ";" << fpty << std::endl;

181

#endif

182

}

183

184

#ifdef _DEBUG_FULL_

185

std::cout << "J = " << j << "/" << newShape->nbPoints << std::endl;

186

std::cout << std::endl;

187

std::cout << "This: " << this << std::endl;

188

#endif

189

190

return newShape;

191

}

192

193

194

void PointSet::splitPolygons( LinkedList<PointSet*> *shapes_toProcess,

195

LinkedList<PointSet*> *shapes_final,

196

double xrm, double yrm , char *uid )

197

{

198

#ifdef _DEBUG_

199

std::cout << "splitPolygons: " << uid << std::endl;

200

#endif

201

int i, j;

202

203

int nbp;

204

double *x;

205

double *y;

206

207

int *pts;

208

209

int *cHull;

210

int cHullSize;

211

212

double cp;

213

double bestcp = 0;

214

215

double bestArea = 0;

216

double area;

217

218

double base;

219

double b, c;

220

double s;

221

222

int ihs;

223

int ihn;

224

225

int ips;

226

int ipn;

227

228

int holeS = -1; // hole start and end points

229

int holeE = -1;

230

231

int retainedPt = -1;

232

int pt = 0;

233

234

double labelArea = xrm * yrm;

235

236

PointSet *shape;

237

238

while ( shapes_toProcess->size() > 0 )

239

{

240

#ifdef _DEBUG_FULL_

241

std::cout << "Shape popping()" << std::endl;

242

#endif

243

shape = shapes_toProcess->pop_front();

244

245

x = shape->x;

246

y = shape->y;

247

nbp = shape->nbPoints;

248

pts = new int[nbp];

249

250

#ifdef _DEBUG_FULL_

251

std::cout << "nbp: " << nbp << std::endl;

252

std::cout << " PtSet: ";

253

#endif

254

255

for ( i = 0; i < nbp; i++ )

256

{

257

pts[i] = i;

258

#ifdef _DEBUG_FULL_

259

std::cout << x[i] << ";" << y[i] << std::endl;

260

#endif

261

}

262

263

#ifdef _DEBUG_FULL_

264

std::cout << std::endl;

265

#endif

266

267

// conpute convex hull

268

shape->cHullSize = convexHullId( pts, x, y, nbp, shape->cHull );

269

270

cHull = shape->cHull;

271

cHullSize = shape->cHullSize;

272

273

274

#ifdef _DEBUG_FULL_

275

std::cout << " CHull: ";

276

for ( i = 0; i < cHullSize; i++ )

277

{

278

std::cout << cHull[i] << " ";

279

}

280

std::cout << std::endl;

281

#endif

282

283

bestArea = 0;

284

retainedPt = -1;

285

286

// lookup for a hole

287

for ( ihs = 0; ihs < cHullSize; ihs++ )

288

{

289

// ihs->ihn => cHull'seg

290

ihn = ( ihs + 1 ) % cHullSize;

291

292

ips = cHull[ihs];

293

ipn = ( ips + 1 ) % nbp;

294

if ( ipn != cHull[ihn] ) // next point on shape is not the next point on cHull => there is a hole here !

295

{

296

bestcp = 0;

297

pt = -1;

298

// lookup for the deepest point in the hole

299

for ( i = ips; i != cHull[ihn]; i = ( i + 1 ) % nbp )

300

{

301

cp = vabs( cross_product( x[cHull[ihs]], y[cHull[ihs]],

302

x[cHull[ihn]], y[cHull[ihn]],

303

x[i], y[i] ) );

304

if ( cp - bestcp > EPSILON )

305

{

306

bestcp = cp;

307

pt = i;

308

}

309

}

310

311

#ifdef _DEBUG_FULL_

312

std::cout << "Deeper POint: " << pt << " between " << ips << " and " << cHull[ihn] << std::endl;

313

#endif

314

if ( pt != -1 )

315

{

316

// compute the ihs->ihn->pt triangle's area

317

base = dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],

318

x[cHull[ihn]], y[cHull[ihn]] );

319

320

b = dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],

321

x[pt], y[pt] );

322

323

c = dist_euc2d( x[cHull[ihn]], y[cHull[ihn]],

324

x[pt], y[pt] );

325

326

s = ( base + b + c ) / 2; // s = half perimeter

327

area = s * ( s - base ) * ( s - b ) * ( s - c );

328

if ( area < 0 )

329

area = -area;

330

331

// retain the biggest area

332

if ( area - bestArea > EPSILON )

333

{

334

bestArea = area;

335

retainedPt = pt;

336

holeS = ihs;

337

holeE = ihn;

338

}

339

}

340

}

341

}

342

343

// we have a hole, its area, and the deppest point in hole

344

// we're going to find the second point to cup the shape

345

// holeS = hole starting point

346

// holeE = hole ending point

347

// retainedPt = deppest point in hole

348

// bestArea = area of triangle HoleS->holeE->retainedPoint

349

bestArea = sqrt( bestArea );

350

double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;

351

int ps = -1, pe = -1, fps = -1, fpe = -1;

352

if ( retainedPt >= 0 && bestArea > labelArea ) // there is a hole so we'll cut the shape in two new shape (only if hole area is bigger than twice labelArea)

353

{

354

#ifdef _DEBUG_FULL_

355

std::cout << "Trou: " << retainedPt << std::endl;

356

std::cout << "Hole: " << cHull[holeS] << " -> " << cHull[holeE] << std::endl;

357

std::cout << "iterate from " << ( cHull[holeE] + 1 ) % nbp << " to " << ( cHull[holeS] - 1 + nbp ) % nbp << std::endl;

358

#endif

359

c = DBL_MAX;

360

361

362

// iterate on all shape points except points which are in the hole

363

double isValid;

364

int k, l;

365

for ( i = ( cHull[holeE] + 1 ) % nbp; i != ( cHull[holeS] - 1 + nbp ) % nbp; i = j )

366

{

367

j = ( i + 1 ) % nbp; // i->j is shape segment not in hole

368

369

// compute distance between retainedPoint and segment

370

// whether perpendicular distance (if retaindPoint is fronting segment i->j)

371

// or distance between retainedPt and i or j (choose the nearest)

372

seg_length = dist_euc2d( x[i], y[i], x[j], y[j] );

373

cx = ( x[i] + x[j] ) / 2.0;

374

cy = ( y[i] + y[j] ) / 2.0;

375

dx = cy - y[i];

376

dy = cx - x[i];

377

378

ex = cx - dx;

379

ey = cy + dy;

380

fx = cx + dx;

381

fy = cy - dy;

382

#ifdef _DEBUG_FULL_

383

std::cout << "D: " << dx << " " << dy << std::endl;

384

std::cout << "seg_length: " << seg_length << std::endl;

385

#endif

386

if ( seg_length < EPSILON || vabs(( b = cross_product( ex, ey, fx, fy, x[retainedPt], y[retainedPt] ) / ( seg_length ) ) ) > ( seg_length / 2 ) ) // retainedPt is not fronting i->j

387

{

388

if (( ex = dist_euc2d_sq( x[i], y[i], x[retainedPt], y[retainedPt] ) ) < ( ey = dist_euc2d_sq( x[j], y[j], x[retainedPt], y[retainedPt] ) ) )

389

{

390

b = ex;

391

ps = i;

392

pe = i;

393

}

394

else

395

{

396

b = ey;

397

ps = j;

398

pe = j;

399

}

400

}

401

else // point fronting i->j => compute pependicular distance => create a new point

402

{

403

b = cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;

404

b *= b;

405

ps = i;

406

pe = j;

407

408

if ( !computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )

409

{

410

std::cout << "Oups ... il devrait par tomber la..." << std::endl;

411

}

412

#ifdef _DEBUG_FULL_

413

std::cout << "intersection : " << x[i] << " " << y[i] << " " << x[j] << " " << y[j] << " " << x[retainedPt] << " " << y[retainedPt] << " " << x[retainedPt] - dx << " " << y[retainedPt] + dy << std::endl;

414

std::cout << " => " << ptx << ";" << pty << std::endl;

415

std::cout << " cxy> " << cx << ";" << cy << std::endl;

416

std::cout << " dxy> " << dx << ";" << dy << std::endl;

417

#endif

418

}

419

420

isValid = true;

421

double pointX, pointY;

422

if ( ps == pe )

423

{

424

pointX = x[pe];

425

pointY = y[pe];

426

}

427

else

428

{

429

pointX = ptx;

430

pointY = pty;

431

}

432

433

for ( k = cHull[holeS]; k != cHull[holeE]; k = ( k + 1 ) % nbp )

434

{

435

l = ( k + 1 ) % nbp;

436

//std::cout << "test " << k << " " << l << std::endl;

437

if ( isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )

438

{

439

isValid = false;

440

//std::cout << "Invalid point" << pe << ps << std::endl;

441

break;

442

}

443

}

444

445

446

if ( isValid && b < c )

447

{

448

//std::cout << "new point: " << ps << " " << pe << std::endl;

449

c = b;

450

fps = ps;

451

fpe = pe;

452

fptx = ptx;

453

fpty = pty;

454

}

455

} // for point which are not in hole

456

457

#ifdef _DEBUG_FULL_

458

std::cout << " cut from " << retainedPt << " to ";

459

if ( fps == fpe )

460

std::cout << "point " << fps << std::endl;

461

else

462

{

463

std::cout << "new point (" << fptx << ";" << fpty << " between " << fps << " and " << fpe << std::endl;

464

}

465

#endif

466

467

// we will cut the shapeu in two new shapes, one from [retainedPoint] to [newPoint] and one form [newPoint] to [retainedPoint]

468

int imin = retainedPt;

469

int imax = ((( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? min( fps, fpe ) : max( fps, fpe ) );

470

471

int nbPtSh1, nbPtSh2; // how many points in new shapes ?

472

if ( imax > imin )

473

nbPtSh1 = imax - imin + 1 + ( fpe != fps );

474

else

475

nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );

476

477

if (( imax == fps ? fpe : fps ) < imin )

478

nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );

479

else

480

nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );

481

482

483

#ifdef _DEBUG_FULL_

484

std::cout << "imin: " << imin << " imax:" << imax << std::endl;

485

#endif

486

if ( retainedPt == -1 || fps == -1 || fpe == -1 )

487

{

488

#ifdef _DEBUG_

489

std::cout << std::endl << "Failed to split feature !!! (uid=" << uid << ")" << std::endl;

490

#endif

491

if ( shape->parent )

492

delete shape;

493

}

494

// check for useless spliting

495

else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )

496

{

497

shapes_final->push_back( shape );

498

}

499

else

500

{

501

502

PointSet *newShape = shape->extractShape( nbPtSh1, imin, imax, fps, fpe, fptx, fpty );

503

504

if ( shape->parent )

505

newShape->parent = shape->parent;

506

else

507

newShape->parent = shape;

508

509

#ifdef _DEBUG_FULL_

510

int i = 0;

511

std::cout << "push back:" << std::endl;

512

for ( i = 0; i < newShape->nbPoints; i++ )

513

{

514

std::cout << newShape->x[i] << ";" << newShape->y[i] << std::endl;

515

}

516

#endif

517

518

shapes_toProcess->push_back( newShape );

519

520

if ( imax == fps )

521

imax = fpe;

522

else

523

imax = fps;

524

525

newShape = shape->extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty );

526

527

if ( shape->parent )

528

newShape->parent = shape->parent;

529

else

530

newShape->parent = shape;

531

532

#ifdef _DEBUG_FULL_

533

std::cout << "push back:" << std::endl;

534

for ( i = 0; i < newShape->nbPoints; i++ )

535

{

536

std::cout << newShape->x[i] << ";" << newShape->y[i] << std::endl;

537

}

538

#endif

539

shapes_toProcess->push_back( newShape );

540

541

if ( shape->parent )

542

delete shape;

543

}

544

}

545

else

546

{

547

#ifdef _DEBUG_FULL_

548

std::cout << "Put shape into shapes_final" << std::endl;

549

#endif

550

shapes_final->push_back( shape );

551

}

552

delete[] pts;

553

}

554

}

555

556

557

PointSet* PointSet::createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside )

558

{

559

#ifdef _DEBUG_FULL_

560

std::cout << "CreateProblemSpecific:" << std::endl;

561

#endif

562

PointSet *shape = new PointSet();

563

shape->x = new double[nbPoints];

564

shape->y = new double[nbPoints];

565

shape->nbPoints = nbPoints;

566

shape->type = type;

567

568

shape->xmin = xmin;

569

shape->xmax = xmax;

570

shape->ymin = ymin;

571

shape->ymax = ymax;

572

573

*inside = true;

574

575

for ( int i = 0; i < nbPoints; i++ )

576

{

577

shape->x[i] = this->x[i];

578

shape->y[i] = this->y[i];

579

580

// check whether it's not outside

581

if ( x[i] < bbmin[0] || x[i] > bbmax[0] || y[i] < bbmin[1] || y[i] > bbmax[1] )

582

*inside = false;

583

}

584

585

shape->holeOf = NULL;

586

shape->parent = NULL;

587

588

return shape;

589

}

590

591

592

593

594

CHullBox * PointSet::compute_chull_bbox()

595

{

596

int i;

597

int j;

598

599

double bbox[4]; // xmin, ymin, xmax, ymax

600

601

double alpha;

602

int alpha_d;

603

604

double alpha_seg;

605

606

double dref;

607

double d1, d2;

608

609

double bb[16]; // {ax, ay, bx, by, cx, cy, dx, dy, ex, ey, fx, fy, gx, gy, hx, hy}}

610

611

double cp;

612

double best_cp;

613

double distNearestPoint;

614

int nearestPoint;

615

616

double area;

617

double width;

618

double length;

619

620

double best_area = DBL_MAX;

621

double best_alpha = -1;

622

double best_bb[16];

623

double best_length = 0;

624

double best_width = 0;

625

626

627

bbox[0] = DBL_MAX;

628

bbox[1] = DBL_MAX;

629

bbox[2] = - DBL_MAX;

630

bbox[3] = - DBL_MAX;

631

632

#ifdef _DEBUG_

633

std::cout << "Compute_chull_bbox" << std::endl;

634

#endif

635

636

for ( i = 0; i < cHullSize; i++ )

637

{

638

#ifdef _DEBUG_FULL_

639

std::cout << x[cHull[i]] << ";" << y[cHull[i]] << std::endl;

640

#endif

641

if ( x[cHull[i]] < bbox[0] )

642

bbox[0] = x[cHull[i]];

643

644

if ( x[cHull[i]] > bbox[2] )

645

bbox[2] = x[cHull[i]];

646

647

if ( y[cHull[i]] < bbox[1] )

648

bbox[1] = y[cHull[i]];

649

650

if ( y[cHull[i]] > bbox[3] )

651

bbox[3] = y[cHull[i]];

652

}

653

654

655

dref = bbox[2] - bbox[0];

656

657

for ( alpha_d = 0; alpha_d < 90; alpha_d++ )

658

{

659

alpha = alpha_d * M_PI / 180.0;

660

d1 = cos( alpha ) * dref;

661

d2 = sin( alpha ) * dref;

662

663

bb[0] = bbox[0];

664

bb[1] = bbox[3]; // ax, ay

665

666

bb[4] = bbox[0];

667

bb[5] = bbox[1]; // cx, cy

668

669

bb[8] = bbox[2];

670

bb[9] = bbox[1]; // ex, ey

671

672

bb[12] = bbox[2];

673

bb[13] = bbox[3]; // gx, gy

674

675

676

bb[2] = bb[0] + d1;

677

bb[3] = bb[1] + d2; // bx, by

678

bb[6] = bb[4] - d2;

679

bb[7] = bb[5] + d1; // dx, dy

680

bb[10] = bb[8] - d1;

681

bb[11] = bb[9] - d2; // fx, fy

682

bb[14] = bb[12] + d2;

683

bb[15] = bb[13] - d1; // hx, hy

684

685

// adjust all points

686

for ( i = 0; i < 16; i += 4 )

687

{

688

689

alpha_seg = (( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) * M_PI / 2 + alpha;

690

691

best_cp = DBL_MAX;

692

nearestPoint = -1;

693

for ( j = 0; j < nbPoints; j++ )

694

{

695

cp = cross_product( bb[i+2], bb[i+3], bb[i], bb[i+1], x[cHull[j]], y[cHull[j]] );

696

if ( cp < best_cp )

697

{

698

best_cp = cp;

699

nearestPoint = cHull[j];

700

}

701

}

702

703

distNearestPoint = best_cp / dref;

704

705

d1 = cos( alpha_seg ) * distNearestPoint;

706

d2 = sin( alpha_seg ) * distNearestPoint;

707

708

bb[i] += d1; // x

709

bb[i+1] += d2; // y

710

bb[i+2] += d1; // x

711

bb[i+3] += d2; // y

712

}

713

714

// compute and compare AREA

715

width = cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;

716

length = cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;

717

718

area = width * length;

719

720

if ( area < 0 )

721

area *= -1;

722

723

724

if ( best_area - area > EPSILON )

725

{

726

best_area = area;

727

best_length = length;

728

best_width = width;

729

best_alpha = alpha;

730

memcpy( best_bb, bb, sizeof( double ) *16 );

731

}

732

}

733

734

// best bbox is defined

735

736

CHullBox * finalBb = new CHullBox();

737

738

for ( i = 0; i < 16; i = i + 4 )

739

{

740

computeLineIntersection( best_bb[i], best_bb[i+1], best_bb[i+2], best_bb[i+3],

741

best_bb[( i+4 ) %16], best_bb[( i+5 ) %16], best_bb[( i+6 ) %16], best_bb[( i+7 ) %16],

742

&finalBb->x[int ( i/4 )], &finalBb->y[int ( i/4 )] );

743

}

744

745

finalBb->alpha = best_alpha;

746

finalBb->width = best_width;

747

finalBb->length = best_length;

748

749

#ifdef _DEBUG_FULL_

750

std::cout << "FINAL" << std::endl;

751

std::cout << "Length : " << best_length << std::endl;

752

std::cout << "Width : " << best_width << std::endl;

753

std::cout << "Alpha: " << best_alpha << " " << best_alpha*180 / M_PI << std::endl;

754

for ( i = 0; i < 4; i++ )

755

{

756

std::cout << finalBb->x[0] << " " << finalBb->y[0] << " ";

757

}

758

std::cout << std::endl;

759

#endif

760

761

return finalBb;

762

}

763

764

#if 0

765

double PointSet::getDistInside( double px, double py )

766

{

767

768

double dist[8];

769

double rpx[8];

770

double rpy[8];

771

bool ok[8];

772

773

if ( !isPointInPolygon( nbPoints, x, y, px, py ) )

774

{

775

double d = getDist( px, py );

776

//std::cout << "Outside : " << d << std::endl;

777

if ( d < 0 )

778

{

779

d = -( d * d * d * d );

780

}

781

else

782

{

783

d = d * d * d * d;

784

}

785

return d;

786

}

787

788

int i;

789

790

double width = ( xmax - xmin ) * 2;

791

double height = ( ymax - ymin ) * 2;

792

793

794

1 0 7

795

\ | /

796

2 --x -- 6

797

/ | \

798

3 4 5

799

800

801

// Compute references points

802

for ( i = 0; i < 8; i++ )

803

{

804

dist[i] = DBL_MAX;

805

ok[i] = false;

806

rpx[i] = px;

807

rpy[i] = py;

808

}

809

810

rpx[0] += 0;

811

rpy[0] += height;

812

813

rpx[1] -= width;

814

rpy[1] += height;

815

816

rpx[2] -= width;

817

rpy[2] += 0;

818

819

rpx[3] -= width;

820

rpy[3] -= height;

821

822

rpx[4] += 0;

823

rpy[4] -= height;

824

825

rpx[5] += width;

826

rpy[5] -= height;

827

828

rpx[6] += width;

829

rpy[6] += 0;

830

831

rpx[7] += width;

832

rpy[7] += height;

833

834

int j, k;

835

for ( i = 0; i < nbPoints; i++ )

836

{

837

j = ( i + 1 ) % nbPoints;

838

839

for ( k = 0; k < 8; k++ )

840

{

841

double ix, iy;

842

if ( computeSegIntersection( px, py, rpx[k], rpy[k], x[i], y[i], x[j], y[j], &ix, &iy ) )

843

{

844

double d = dist_euc2d_sq( px, py, ix, iy );

845

if ( dist[k] > d )

846

{

847

dist[k] = d;

848

ok[k] = true;

849

//std::cout << "new dist for " << k << ": " << dist[k] << std::endl;

850

}

851

}

852

}

853

}

854

855

double a, b, c, d;

856

857

858

for ( i = 0; i < 8; i++ )

859

{

860

if ( !ok[i] )

861

{

862

std::cout << "ERROR!!!!!!!!!!!!!!!!!" << std::endl;

863

dist[i] = 0;

864

}

865

//std::cout << "dist[" << i << "]: " << dist[i] << std::endl;

866

}

867

868

a = min( dist[0], dist[4] );

869

b = min( dist[1], dist[5] );

870

c = min( dist[2], dist[6] );

871

d = min( dist[3], dist[7] );

872

873

std::cout << "a: " << a << std::endl;

874

std::cout << "b: " << b << std::endl;

875

std::cout << "c: " << c << std::endl;

876

std::cout << "d: " << d << std::endl;

877

878

//a = (a+b+c+d)/4.0;

879

880

//a = min(a,b);

881

//c = min(c,d);

882

//return min(a,c);

883

884

885

return ( a*b*c*d );

886

}

887

#endif

888

889

double PointSet::getDist( double px, double py, double *rx, double *ry )

890

{

891

if ( nbPoints == 1 || type == GEOS_POINT )

892

{

893

if ( rx && ry )

894

{

895

*rx = x[0];

896

*ry = y[0];

897

}

898

return dist_euc2d_sq( x[0], y[0], px, py );

899

}

900

901

int a, b;

902

int nbP = ( type == GEOS_POLYGON ? nbPoints : nbPoints - 1 );

903

904

double best_dist = DBL_MAX;

905

double d;

906

907

for ( a = 0; a < nbP; a++ )

908

{

909

b = ( a + 1 ) % nbPoints;

910

911

double px2, py2;

912

px2 = px - y[b] + y[a];

913

py2 = py + x[b] - x[a];

914

double ix, iy;

915

916

// (px,py)->(px2,py2) is a line perpendicular to a->b

917

// Check the line p->p2 cross the segment a->b

918

if ( computeLineSegIntersection( px, py, px2, py2,

919

x[a], y[a], x[b], y[b],

920

&ix, &iy ) )

921

{

922

d = dist_euc2d_sq( px, py, ix, iy );

923

}

924

else

925

{

926

double d1 = dist_euc2d_sq( x[a], y[a], px, py );

927

double d2 = dist_euc2d_sq( x[b], y[b], px, py );

928

if ( d1 < d2 )

929

{

930

d = d1;

931

ix = x[a];

932

iy = y[a];

933

}

934

else

935

{

936

d = d2;

937

ix = x[b];

938

iy = y[b];

939

}

940

}

941

942

if ( d < best_dist )

943

{

944

best_dist = d;

945

if ( rx && ry )

946

{

947

*rx = ix;

948

*ry = iy;

949

}

950

}

951

} // end for (a in nbPoints)

952

953

return best_dist;

954

}

955

956

957

958

void PointSet::getCentroid( double &px, double &py )

959

{

960

// for explanation see this page:

961

// http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/

962

963

int i, j;

964

double cx = 0, cy = 0, A = 0, tmp;

965

for ( i = 0; i < nbPoints; i++ )

966

{

967

j = i + 1; if ( j == nbPoints ) j = 0;

968

tmp = ( x[i] * y[j] - x[j] * y[i] );

969

cx += ( x[i] + x[j] ) * tmp;

970

cy += ( y[i] + y[j] ) * tmp;

971

A += tmp;

972

}

973

974

px = cx / ( 3 * A );

975

py = cy / ( 3 * A );

976

}

977

978

} // end namespace

979

980

#endif

Older »