matplotlib/src/tri/_tri.cpp at main · commit-0/matplotlib · GitHub

History

2071 lines (1802 loc) · 69.3 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

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

308

309

310

311

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

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

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

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

691

692

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

735

736

737

738

739

740

741

742

743

744

745

746

747

748

749

750

751

752

753

754

755

756

757

758

759

760

761

762

763

764

765

766

767

768

769

770

771

772

773

774

775

776

777

778

779

780

781

782

783

784

785

786

787

788

789

790

791

792

793

794

795

796

797

798

799

800

801

802

803

804

805

806

807

808

809

810

811

812

813

814

815

816

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

866

867

868

869

870

871

872

873

874

875

876

877

878

879

880

881

882

883

884

885

886

887

888

889

890

891

892

893

894

895

896

897

898

899

900

901

902

903

904

905

906

907

908

909

910

911

912

913

914

915

916

917

918

919

920

921

922

923

924

925

926

927

928

929

930

931

932

933

934

935

936

937

938

939

940

941

942

943

944

945

946

947

948

949

950

951

952

953

954

955

956

957

958

959

960

961

962

963

964

965

966

967

968

969

970

971

972

973

974

975

976

977

978

979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

/* This file contains liberal use of asserts to assist code development and debugging.

* Standard Matplotlib builds disable asserts so they cause no performance reduction. To

* enable the asserts, you need to undefine the NDEBUG macro, which is achieved by

* passing ``b_ndebug=false`` to the Meson configuration.

#include "../mplutils.h"

#include "_tri.h"

#include <algorithm>

#include <random>

#include <set>

TriEdge::TriEdge()

: tri(-1), edge(-1)

{}

TriEdge::TriEdge(int tri_, int edge_)

: tri(tri_), edge(edge_)

{}

bool TriEdge::operator<(const TriEdge& other) const

{

if (tri != other.tri)

return tri < other.tri;

else

return edge < other.edge;

}

bool TriEdge::operator==(const TriEdge& other) const

{

return tri == other.tri && edge == other.edge;

}

bool TriEdge::operator!=(const TriEdge& other) const

{

return !operator==(other);

}

std::ostream& operator<<(std::ostream& os, const TriEdge& tri_edge)

{

return os << tri_edge.tri << ' ' << tri_edge.edge;

}

XY::XY()

{}

XY::XY(const double& x_, const double& y_)

: x(x_), y(y_)

{}

double XY::angle() const

{

return atan2(y, x);

}

double XY::cross_z(const XY& other) const

{

return x*other.y - y*other.x;

}

bool XY::is_right_of(const XY& other) const

{

if (x == other.x)

return y > other.y;

else

return x > other.x;

}

bool XY::operator==(const XY& other) const

{

return x == other.x && y == other.y;

}

bool XY::operator!=(const XY& other) const

{

return x != other.x || y != other.y;

}

XY XY::operator*(const double& multiplier) const

{

return XY(x*multiplier, y*multiplier);

}

const XY& XY::operator+=(const XY& other)

{

x += other.x;

y += other.y;

return *this;

}

const XY& XY::operator-=(const XY& other)

{

x -= other.x;

y -= other.y;

return *this;

}

XY XY::operator+(const XY& other) const

{

return XY(x + other.x, y + other.y);

}

XY XY::operator-(const XY& other) const

{

return XY(x - other.x, y - other.y);

}

std::ostream& operator<<(std::ostream& os, const XY& xy)

{

return os << '(' << xy.x << ' ' << xy.y << ')';

}

XYZ::XYZ(const double& x_, const double& y_, const double& z_)

: x(x_), y(y_), z(z_)

{}

XYZ XYZ::cross(const XYZ& other) const

{

return XYZ(y*other.z - z*other.y,

z*other.x - x*other.z,

x*other.y - y*other.x);

}

double XYZ::dot(const XYZ& other) const

{

return x*other.x + y*other.y + z*other.z;

}

XYZ XYZ::operator-(const XYZ& other) const

{

return XYZ(x - other.x, y - other.y, z - other.z);

}

std::ostream& operator<<(std::ostream& os, const XYZ& xyz)

{

return os << '(' << xyz.x << ' ' << xyz.y << ' ' << xyz.z << ')';

}

BoundingBox::BoundingBox()

: empty(true), lower(0.0, 0.0), upper(0.0, 0.0)

{}

void BoundingBox::add(const XY& point)

{

if (empty) {

empty = false;

lower = upper = point;

} else {

if (point.x < lower.x) lower.x = point.x;

else if (point.x > upper.x) upper.x = point.x;

if (point.y < lower.y) lower.y = point.y;

else if (point.y > upper.y) upper.y = point.y;

}

void BoundingBox::expand(const XY& delta)

{

if (!empty) {

lower -= delta;

upper += delta;

}

ContourLine::ContourLine()

: std::vector<XY>()

{}

void ContourLine::push_back(const XY& point)

{

if (empty() || point != back())

std::vector<XY>::push_back(point);

}

void ContourLine::write() const

{

std::cout << "ContourLine of " << size() << " points:";

for (const_iterator it = begin(); it != end(); ++it)

std::cout << ' ' << *it;

std::cout << std::endl;

}

void write_contour(const Contour& contour)

{

std::cout << "Contour of " << contour.size() << " lines." << std::endl;

for (Contour::const_iterator it = contour.begin(); it != contour.end(); ++it)

it->write();

}

Triangulation::Triangulation(const CoordinateArray& x,

const CoordinateArray& y,

const TriangleArray& triangles,

const MaskArray& mask,

const EdgeArray& edges,

const NeighborArray& neighbors,

bool correct_triangle_orientations)

: _x(x),

_y(y),

_triangles(triangles),

_mask(mask),

_edges(edges),

_neighbors(neighbors)

{

if (_x.ndim() != 1 || _y.ndim() != 1 || _x.shape(0) != _y.shape(0))

throw std::invalid_argument("x and y must be 1D arrays of the same length");

if (_triangles.ndim() != 2 || _triangles.shape(1) != 3)

throw std::invalid_argument("triangles must be a 2D array of shape (?,3)");

// Optional mask.

if (_mask.size() > 0 &&

(_mask.ndim() != 1 || _mask.shape(0) != _triangles.shape(0)))

throw std::invalid_argument(

"mask must be a 1D array with the same length as the triangles array");

// Optional edges.

if (_edges.size() > 0 &&

(_edges.ndim() != 2 || _edges.shape(1) != 2))

throw std::invalid_argument("edges must be a 2D array with shape (?,2)");

// Optional neighbors.

if (_neighbors.size() > 0 &&

(_neighbors.ndim() != 2 || _neighbors.shape() != _triangles.shape()))

throw std::invalid_argument(

"neighbors must be a 2D array with the same shape as the triangles array");

if (correct_triangle_orientations)

correct_triangles();

}

void Triangulation::calculate_boundaries()

{

get_neighbors(); // Ensure _neighbors has been created.

// Create set of all boundary TriEdges, which are those which do not

// have a neighbor triangle.

typedef std::set<TriEdge> BoundaryEdges;

BoundaryEdges boundary_edges;

for (int tri = 0; tri < get_ntri(); ++tri) {

if (!is_masked(tri)) {

for (int edge = 0; edge < 3; ++edge) {

if (get_neighbor(tri, edge) == -1) {

boundary_edges.insert(TriEdge(tri, edge));

}

// Take any boundary edge and follow the boundary until return to start

// point, removing edges from boundary_edges as they are used. At the same

// time, initialise the _tri_edge_to_boundary_map.

while (!boundary_edges.empty()) {

// Start of new boundary.

BoundaryEdges::iterator it = boundary_edges.begin();

int tri = it->tri;

int edge = it->edge;

_boundaries.push_back(Boundary());

Boundary& boundary = _boundaries.back();

while (true) {

boundary.push_back(TriEdge(tri, edge));

boundary_edges.erase(it);

_tri_edge_to_boundary_map[TriEdge(tri, edge)] =

BoundaryEdge(_boundaries.size()-1, boundary.size()-1);

// Move to next edge of current triangle.

edge = (edge+1) % 3;

// Find start point index of boundary edge.

int point = get_triangle_point(tri, edge);

// Find next TriEdge by traversing neighbors until find one

// without a neighbor.

while (get_neighbor(tri, edge) != -1) {

tri = get_neighbor(tri, edge);

edge = get_edge_in_triangle(tri, point);

}

if (TriEdge(tri,edge) == boundary.front())

break; // Reached beginning of this boundary, so finished it.

else

it = boundary_edges.find(TriEdge(tri, edge));

}

void Triangulation::calculate_edges()

{

assert(!has_edges() && "Expected empty edges array");

// Create set of all edges, storing them with start point index less than

// end point index.

typedef std::set<Edge> EdgeSet;

EdgeSet edge_set;

for (int tri = 0; tri < get_ntri(); ++tri) {

if (!is_masked(tri)) {

for (int edge = 0; edge < 3; edge++) {

int start = get_triangle_point(tri, edge);

int end = get_triangle_point(tri, (edge+1)%3);

edge_set.insert(start > end ? Edge(start,end) : Edge(end,start));

}

// Convert to python _edges array.

py::ssize_t dims[2] = {static_cast<py::ssize_t>(edge_set.size()), 2};

_edges = EdgeArray(dims);

auto edges = _edges.mutable_data();

int i = 0;

for (EdgeSet::const_iterator it = edge_set.begin(); it != edge_set.end(); ++it) {

edges[i++] = it->start;

edges[i++] = it->end;

}

void Triangulation::calculate_neighbors()

{

assert(!has_neighbors() && "Expected empty neighbors array");

// Create _neighbors array with shape (ntri,3) and initialise all to -1.

py::ssize_t dims[2] = {get_ntri(), 3};

_neighbors = NeighborArray(dims);

auto* neighbors = _neighbors.mutable_data();

int tri, edge;

std::fill(neighbors, neighbors+3*get_ntri(), -1);

// For each triangle edge (start to end point), find corresponding neighbor

// edge from end to start point. Do this by traversing all edges and

// storing them in a map from edge to TriEdge. If corresponding neighbor

// edge is already in the map, don't need to store new edge as neighbor

// already found.

typedef std::map<Edge, TriEdge> EdgeToTriEdgeMap;

EdgeToTriEdgeMap edge_to_tri_edge_map;

for (tri = 0; tri < get_ntri(); ++tri) {

if (!is_masked(tri)) {

for (edge = 0; edge < 3; ++edge) {

int start = get_triangle_point(tri, edge);

int end = get_triangle_point(tri, (edge+1)%3);

EdgeToTriEdgeMap::iterator it =

edge_to_tri_edge_map.find(Edge(end,start));

if (it == edge_to_tri_edge_map.end()) {

// No neighbor edge exists in the edge_to_tri_edge_map, so

// add this edge to it.

edge_to_tri_edge_map[Edge(start,end)] = TriEdge(tri,edge);

} else {

// Neighbor edge found, set the two elements of _neighbors

// and remove edge from edge_to_tri_edge_map.

neighbors[3*tri + edge] = it->second.tri;

neighbors[3*it->second.tri + it->second.edge] = tri;

edge_to_tri_edge_map.erase(it);

}

// Note that remaining edges in the edge_to_tri_edge_map correspond to

// boundary edges, but the boundaries are calculated separately elsewhere.

}

Triangulation::TwoCoordinateArray Triangulation::calculate_plane_coefficients(

const CoordinateArray& z)

{

if (z.ndim() != 1 || z.shape(0) != _x.shape(0))

throw std::invalid_argument(

"z must be a 1D array with the same length as the triangulation x and y arrays");

int dims[2] = {get_ntri(), 3};

Triangulation::TwoCoordinateArray planes_array(dims);

auto planes = planes_array.mutable_unchecked<2>();

auto triangles = _triangles.unchecked<2>();

auto x = _x.unchecked<1>();

auto y = _y.unchecked<1>();

auto z_ptr = z.unchecked<1>();

int point;

for (int tri = 0; tri < get_ntri(); ++tri) {

if (is_masked(tri)) {

planes(tri, 0) = 0.0;

planes(tri, 1) = 0.0;

planes(tri, 2) = 0.0;

}

else {

// Equation of plane for all points r on plane is r.normal = p

// where normal is vector normal to the plane, and p is a

// constant. Rewrite as

// r_x*normal_x + r_y*normal_y + r_z*normal_z = p

// and rearrange to give

// r_z = (-normal_x/normal_z)*r_x + (-normal_y/normal_z)*r_y +

// p/normal_z

point = triangles(tri, 0);

XYZ point0(x(point), y(point), z_ptr(point));

point = triangles(tri, 1);

XYZ side01 = XYZ(x(point), y(point), z_ptr(point)) - point0;

point = triangles(tri, 2);

XYZ side02 = XYZ(x(point), y(point), z_ptr(point)) - point0;

XYZ normal = side01.cross(side02);

if (normal.z == 0.0) {

// Normal is in x-y plane which means triangle consists of

// colinear points. To avoid dividing by zero, we use the

// Moore-Penrose pseudo-inverse.

double sum2 = (side01.x*side01.x + side01.y*side01.y +

side02.x*side02.x + side02.y*side02.y);

double a = (side01.x*side01.z + side02.x*side02.z) / sum2;

double b = (side01.y*side01.z + side02.y*side02.z) / sum2;

planes(tri, 0) = a;

planes(tri, 1) = b;

planes(tri, 2) = point0.z - a*point0.x - b*point0.y;

}

else {

planes(tri, 0) = -normal.x / normal.z; // x

planes(tri, 1) = -normal.y / normal.z; // y

planes(tri, 2) = normal.dot(point0) / normal.z; // constant

}

return planes_array;

}

void Triangulation::correct_triangles()

{

auto triangles = _triangles.mutable_data();

auto neighbors = _neighbors.mutable_data();

for (int tri = 0; tri < get_ntri(); ++tri) {

XY point0 = get_point_coords(triangles[3*tri]);

XY point1 = get_point_coords(triangles[3*tri+1]);

XY point2 = get_point_coords(triangles[3*tri+2]);

if ( (point1 - point0).cross_z(point2 - point0) < 0.0) {

// Triangle points are clockwise, so change them to anticlockwise.

std::swap(triangles[3*tri+1], triangles[3*tri+2]);

if (has_neighbors())

std::swap(neighbors[3*tri+1], neighbors[3*tri+2]);

}

const Triangulation::Boundaries& Triangulation::get_boundaries() const

{

if (_boundaries.empty())

const_cast<Triangulation*>(this)->calculate_boundaries();

return _boundaries;

}

void Triangulation::get_boundary_edge(const TriEdge& triEdge,

int& boundary,

int& edge) const

{

get_boundaries(); // Ensure _tri_edge_to_boundary_map has been created.

TriEdgeToBoundaryMap::const_iterator it =

_tri_edge_to_boundary_map.find(triEdge);

assert(it != _tri_edge_to_boundary_map.end() &&

"TriEdge is not on a boundary");

boundary = it->second.boundary;

edge = it->second.edge;

}

int Triangulation::get_edge_in_triangle(int tri, int point) const

{

assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds");

assert(point >= 0 && point < get_npoints() && "Point index out of bounds.");

auto triangles = _triangles.data();

for (int edge = 0; edge < 3; ++edge) {

if (triangles[3*tri + edge] == point)

return edge;

}

return -1; // point is not in triangle.

}

Triangulation::EdgeArray& Triangulation::get_edges()

{

if (!has_edges())

calculate_edges();

return _edges;

}

int Triangulation::get_neighbor(int tri, int edge) const

{

assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds");

assert(edge >= 0 && edge < 3 && "Edge index out of bounds");

if (!has_neighbors())

const_cast<Triangulation&>(*this).calculate_neighbors();

return _neighbors.data()[3*tri + edge];

}

TriEdge Triangulation::get_neighbor_edge(int tri, int edge) const

{

int neighbor_tri = get_neighbor(tri, edge);

if (neighbor_tri == -1)

return TriEdge(-1,-1);

else

return TriEdge(neighbor_tri,

get_edge_in_triangle(neighbor_tri,

get_triangle_point(tri,

(edge+1)%3)));

}

Triangulation::NeighborArray& Triangulation::get_neighbors()

{

if (!has_neighbors())

calculate_neighbors();

return _neighbors;

}

int Triangulation::get_npoints() const

{

return _x.shape(0);

}

int Triangulation::get_ntri() const

{

return _triangles.shape(0);

}

XY Triangulation::get_point_coords(int point) const

{

assert(point >= 0 && point < get_npoints() && "Point index out of bounds.");

return XY(_x.data()[point], _y.data()[point]);

}

int Triangulation::get_triangle_point(int tri, int edge) const

{

assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds");

assert(edge >= 0 && edge < 3 && "Edge index out of bounds");

return _triangles.data()[3*tri + edge];

}

int Triangulation::get_triangle_point(const TriEdge& tri_edge) const

{

return get_triangle_point(tri_edge.tri, tri_edge.edge);

}

bool Triangulation::has_edges() const

{

return _edges.size() > 0;

}

bool Triangulation::has_mask() const

{

return _mask.size() > 0;

}

bool Triangulation::has_neighbors() const

{

return _neighbors.size() > 0;

}

bool Triangulation::is_masked(int tri) const

{

assert(tri >= 0 && tri < get_ntri() && "Triangle index out of bounds.");

return has_mask() && _mask.data()[tri];

}

void Triangulation::set_mask(const MaskArray& mask)

{

if (mask.size() > 0 &&

(mask.ndim() != 1 || mask.shape(0) != _triangles.shape(0)))

throw std::invalid_argument(

"mask must be a 1D array with the same length as the triangles array");

_mask = mask;

// Clear derived fields so they are recalculated when needed.

_edges = EdgeArray();

_neighbors = NeighborArray();

_boundaries.clear();

}

void Triangulation::write_boundaries() const

{

const Boundaries& bs = get_boundaries();

std::cout << "Number of boundaries: " << bs.size() << std::endl;

for (Boundaries::const_iterator it = bs.begin(); it != bs.end(); ++it) {

const Boundary& b = *it;

std::cout << " Boundary of " << b.size() << " points: ";

for (Boundary::const_iterator itb = b.begin(); itb != b.end(); ++itb) {

std::cout << *itb << ", ";

}

std::cout << std::endl;

}

TriContourGenerator::TriContourGenerator(Triangulation& triangulation,

const CoordinateArray& z)

: _triangulation(triangulation),

_z(z),

_interior_visited(2*_triangulation.get_ntri()),

_boundaries_visited(0),

_boundaries_used(0)

{

if (_z.ndim() != 1 || _z.shape(0) != _triangulation.get_npoints())

throw std::invalid_argument(

"z must be a 1D array with the same length as the x and y arrays");

}

void TriContourGenerator::clear_visited_flags(bool include_boundaries)

{

// Clear _interiorVisited.

std::fill(_interior_visited.begin(), _interior_visited.end(), false);

if (include_boundaries) {

if (_boundaries_visited.empty()) {

const Boundaries& boundaries = get_boundaries();

// Initialise _boundaries_visited.

_boundaries_visited.reserve(boundaries.size());

for (Boundaries::const_iterator it = boundaries.begin();

it != boundaries.end(); ++it)

_boundaries_visited.push_back(BoundaryVisited(it->size()));

// Initialise _boundaries_used.

_boundaries_used = BoundariesUsed(boundaries.size());

}

// Clear _boundaries_visited.

for (BoundariesVisited::iterator it = _boundaries_visited.begin();

it != _boundaries_visited.end(); ++it)

std::fill(it->begin(), it->end(), false);

// Clear _boundaries_used.

std::fill(_boundaries_used.begin(), _boundaries_used.end(), false);

}

py::tuple TriContourGenerator::contour_line_to_segs_and_kinds(const Contour& contour)

{

// Convert all of the lines generated by a call to create_contour() into

// their Python equivalents for return to the calling function.

// A line is either a closed line loop (in which case the last point is

// identical to the first) or an open line strip. Two NumPy arrays are

// created for each line:

// vertices is a double array of shape (npoints, 2) containing the (x, y)

// coordinates of the points in the line

// codes is a uint8 array of shape (npoints,) containing the 'kind codes'

// which are defined in the Path class

// and they are appended to the Python lists vertices_list and codes_list

// respectively for return to the Python calling function.

py::list vertices_list(contour.size());

py::list codes_list(contour.size());

for (Contour::size_type i = 0; i < contour.size(); ++i) {

const ContourLine& contour_line = contour[i];

py::ssize_t npoints = static_cast<py::ssize_t>(contour_line.size());

py::ssize_t segs_dims[2] = {npoints, 2};

CoordinateArray segs(segs_dims);

double* segs_ptr = segs.mutable_data();

py::ssize_t codes_dims[1] = {npoints};

CodeArray codes(codes_dims);

unsigned char* codes_ptr = codes.mutable_data();

for (ContourLine::const_iterator it = contour_line.begin();

it != contour_line.end(); ++it) {

*segs_ptr++ = it->x;

*segs_ptr++ = it->y;

*codes_ptr++ = (it == contour_line.begin() ? MOVETO : LINETO);

}

// Closed line loop has identical first and last (x, y) points.

if (contour_line.size() > 1 &&

contour_line.front() == contour_line.back())

*(codes_ptr-1) = CLOSEPOLY;

vertices_list[i] = segs;

codes_list[i] = codes;

}

return py::make_tuple(vertices_list, codes_list);

}

py::tuple TriContourGenerator::contour_to_segs_and_kinds(const Contour& contour)

{

// Convert all of the polygons generated by a call to

// create_filled_contour() into their Python equivalents for return to the

// calling function. All of the polygons' points and kinds codes are

// combined into single NumPy arrays for each; this avoids having

// to determine which polygons are holes as this will be determined by the

// renderer. If there are ntotal points in all of the polygons, the two

// NumPy arrays created are:

// vertices is a double array of shape (ntotal, 2) containing the (x, y)

// coordinates of the points in the polygons

// codes is a uint8 array of shape (ntotal,) containing the 'kind codes'

// which are defined in the Path class

// and they are returned in the Python lists vertices_list and codes_list

// respectively.

Contour::const_iterator line;

ContourLine::const_iterator point;

// Find total number of points in all contour lines.

py::ssize_t n_points = 0;

for (line = contour.begin(); line != contour.end(); ++line)

n_points += static_cast<py::ssize_t>(line->size());

// Create segs array for point coordinates.

py::ssize_t segs_dims[2] = {n_points, 2};

TwoCoordinateArray segs(segs_dims);

double* segs_ptr = segs.mutable_data();

// Create kinds array for code types.

py::ssize_t codes_dims[1] = {n_points};

CodeArray codes(codes_dims);

unsigned char* codes_ptr = codes.mutable_data();

for (line = contour.begin(); line != contour.end(); ++line) {

for (point = line->begin(); point != line->end(); point++) {

*segs_ptr++ = point->x;

*segs_ptr++ = point->y;

*codes_ptr++ = (point == line->begin() ? MOVETO : LINETO);

}

if (line->size() > 1)

*(codes_ptr-1) = CLOSEPOLY;

}

py::list vertices_list(1);

vertices_list[0] = segs;

py::list codes_list(1);

codes_list[0] = codes;

return py::make_tuple(vertices_list, codes_list);

}

py::tuple TriContourGenerator::create_contour(const double& level)

{

clear_visited_flags(false);

Contour contour;

find_boundary_lines(contour, level);

find_interior_lines(contour, level, false);

return contour_line_to_segs_and_kinds(contour);

}

py::tuple TriContourGenerator::create_filled_contour(const double& lower_level,

const double& upper_level)

{

if (lower_level >= upper_level)

throw std::invalid_argument("filled contour levels must be increasing");

clear_visited_flags(true);

Contour contour;

find_boundary_lines_filled(contour, lower_level, upper_level);

find_interior_lines(contour, lower_level, false);

find_interior_lines(contour, upper_level, true);

return contour_to_segs_and_kinds(contour);

}

XY TriContourGenerator::edge_interp(int tri, int edge, const double& level)

{

return interp(_triangulation.get_triangle_point(tri, edge),

_triangulation.get_triangle_point(tri, (edge+1)%3),

level);

}

void TriContourGenerator::find_boundary_lines(Contour& contour,

const double& level)

{

// Traverse boundaries to find starting points for all contour lines that

// intersect the boundaries. For each starting point found, follow the

// line to its end before continuing.

const Triangulation& triang = _triangulation;

const Boundaries& boundaries = get_boundaries();

for (Boundaries::const_iterator it = boundaries.begin();

it != boundaries.end(); ++it) {

const Boundary& boundary = *it;

bool startAbove, endAbove = false;

for (Boundary::const_iterator itb = boundary.begin();

itb != boundary.end(); ++itb) {

if (itb == boundary.begin())

startAbove = get_z(triang.get_triangle_point(*itb)) >= level;

else

startAbove = endAbove;

endAbove = get_z(triang.get_triangle_point(itb->tri,

(itb->edge+1)%3)) >= level;

if (startAbove && !endAbove) {

// This boundary edge is the start point for a contour line,

// so follow the line.

contour.push_back(ContourLine());

ContourLine& contour_line = contour.back();

TriEdge tri_edge = *itb;

follow_interior(contour_line, tri_edge, true, level, false);

}

void TriContourGenerator::find_boundary_lines_filled(Contour& contour,

const double& lower_level,

const double& upper_level)

{

// Traverse boundaries to find starting points for all contour lines that

// intersect the boundaries. For each starting point found, follow the

// line to its end before continuing.

const Triangulation& triang = _triangulation;

const Boundaries& boundaries = get_boundaries();

for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {

const Boundary& boundary = boundaries[i];

for (Boundary::size_type j = 0; j < boundary.size(); ++j) {

if (!_boundaries_visited[i][j]) {

// z values of start and end of this boundary edge.

double z_start = get_z(triang.get_triangle_point(boundary[j]));

double z_end = get_z(triang.get_triangle_point(

boundary[j].tri, (boundary[j].edge+1)%3));

// Does this boundary edge's z increase through upper level

// and/or decrease through lower level?

bool incr_upper = (z_start < upper_level && z_end >= upper_level);

bool decr_lower = (z_start >= lower_level && z_end < lower_level);

if (decr_lower || incr_upper) {

// Start point for contour line, so follow it.

contour.push_back(ContourLine());

ContourLine& contour_line = contour.back();

TriEdge start_tri_edge = boundary[j];

TriEdge tri_edge = start_tri_edge;

// Traverse interior and boundaries until return to start.

bool on_upper = incr_upper;

do {

follow_interior(contour_line, tri_edge, true,

on_upper ? upper_level : lower_level, on_upper);

on_upper = follow_boundary(contour_line, tri_edge,

lower_level, upper_level, on_upper);

} while (tri_edge != start_tri_edge);

// Close polygon.

contour_line.push_back(contour_line.front());

}

// Add full boundaries that lie between the lower and upper levels. These

// are boundaries that have not been touched by an internal contour line

// which are stored in _boundaries_used.

for (Boundaries::size_type i = 0; i < boundaries.size(); ++i) {

if (!_boundaries_used[i]) {

const Boundary& boundary = boundaries[i];

double z = get_z(triang.get_triangle_point(boundary[0]));

if (z >= lower_level && z < upper_level) {

contour.push_back(ContourLine());

ContourLine& contour_line = contour.back();

for (Boundary::size_type j = 0; j < boundary.size(); ++j)

contour_line.push_back(triang.get_point_coords(

triang.get_triangle_point(boundary[j])));

// Close polygon.

contour_line.push_back(contour_line.front());

}

void TriContourGenerator::find_interior_lines(Contour& contour,

const double& level,

bool on_upper)

{

const Triangulation& triang = _triangulation;

int ntri = triang.get_ntri();

for (int tri = 0; tri < ntri; ++tri) {

int visited_index = (on_upper ? tri+ntri : tri);

if (_interior_visited[visited_index] || triang.is_masked(tri))

continue; // Triangle has already been visited or is masked.

_interior_visited[visited_index] = true;

// Determine edge via which to leave this triangle.

int edge = get_exit_edge(tri, level, on_upper);

assert(edge >= -1 && edge < 3 && "Invalid exit edge");

if (edge == -1)

continue; // Contour does not pass through this triangle.

// Found start of new contour line loop.

contour.push_back(ContourLine());

ContourLine& contour_line = contour.back();

TriEdge tri_edge = triang.get_neighbor_edge(tri, edge);

follow_interior(contour_line, tri_edge, false, level, on_upper);

// Close line loop

contour_line.push_back(contour_line.front());

}

bool TriContourGenerator::follow_boundary(ContourLine& contour_line,

TriEdge& tri_edge,

const double& lower_level,

const double& upper_level,

bool on_upper)

{

const Triangulation& triang = _triangulation;

const Boundaries& boundaries = get_boundaries();

// Have TriEdge to start at, need equivalent boundary edge.

int boundary, edge;

triang.get_boundary_edge(tri_edge, boundary, edge);

_boundaries_used[boundary] = true;

bool stop = false;

bool first_edge = true;

double z_start, z_end = 0;

while (!stop)

{

assert(!_boundaries_visited[boundary][edge] && "Boundary already visited");

_boundaries_visited[boundary][edge] = true;

// z values of start and end points of boundary edge.

if (first_edge)

z_start = get_z(triang.get_triangle_point(tri_edge));

else

z_start = z_end;

z_end = get_z(triang.get_triangle_point(tri_edge.tri,

(tri_edge.edge+1)%3));

if (z_end > z_start) { // z increasing.

if (!(!on_upper && first_edge) &&

z_end >= lower_level && z_start < lower_level) {

stop = true;

on_upper = false;

} else if (z_end >= upper_level && z_start < upper_level) {

stop = true;

on_upper = true;

}

} else { // z decreasing.

if (!(on_upper && first_edge) &&

z_start >= upper_level && z_end < upper_level) {

stop = true;

on_upper = true;

} else if (z_start >= lower_level && z_end < lower_level) {

stop = true;

on_upper = false;

}

first_edge = false;

if (!stop) {

// Move to next boundary edge, adding point to contour line.

edge = (edge+1) % (int)boundaries[boundary].size();

tri_edge = boundaries[boundary][edge];

contour_line.push_back(triang.get_point_coords(

triang.get_triangle_point(tri_edge)));

}

return on_upper;

}

void TriContourGenerator::follow_interior(ContourLine& contour_line,

TriEdge& tri_edge,

bool end_on_boundary,

const double& level,

bool on_upper)

{

int& tri = tri_edge.tri;

int& edge = tri_edge.edge;

// Initial point.

contour_line.push_back(edge_interp(tri, edge, level));

while (true) {

int visited_index = tri;

if (on_upper)

visited_index += _triangulation.get_ntri();

// Check for end not on boundary.

if (!end_on_boundary && _interior_visited[visited_index])

break; // Reached start point, so return.

// Determine edge by which to leave this triangle.

edge = get_exit_edge(tri, level, on_upper);

assert(edge >= 0 && edge < 3 && "Invalid exit edge");

View remainder of file in raw view

matplotlib/src/tri/_tri.cpp at main · commit-0/matplotlib · GitHub

Uh oh!

Expand file tree

Breadcrumbs

_tri.cpp

Latest commit

History

Breadcrumbs

_tri.cpp

File metadata and controls