123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470 |
- // David Eberly, Geometric Tools, Redmond WA 98052
- // Copyright (c) 1998-2020
- // Distributed under the Boost Software License, Version 1.0.
- // https://www.boost.org/LICENSE_1_0.txt
- // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
- // Version: 4.0.2019.08.13
- #pragma once
- #include <Mathematics/CurveExtractor.h>
- // The level set extraction algorithm implemented here is described
- // in Section 3 of the document
- // https://www.geometrictools.com/Documentation/ExtractLevelCurves.pdf
- namespace WwiseGTE
- {
- // The image type T must be one of the integer types: int8_t, int16_t,
- // int32_t, uint8_t, uint16_t or uint32_t. Internal integer computations
- // are performed using int64_t. The type Real is for extraction to
- // floating-point vertices.
- template <typename T, typename Real>
- class CurveExtractorSquares : public CurveExtractor<T, Real>
- {
- public:
- // Convenience type definitions.
- typedef typename CurveExtractor<T, Real>::Vertex Vertex;
- typedef typename CurveExtractor<T, Real>::Edge Edge;
- // The input is a 2D image with lexicographically ordered pixels (x,y)
- // stored in a linear array. Pixel (x,y) is stored in the array at
- // location index = x + xBound * y. The inputs xBound and yBound must
- // each be 2 or larger so that there is at least one image square to
- // process. The inputPixels must be nonnull and point to contiguous
- // storage that contains at least xBound * yBound elements.
- CurveExtractorSquares(int xBound, int yBound, T const* inputPixels)
- :
- CurveExtractor<T, Real>(xBound, yBound, inputPixels)
- {
- }
- // Extract level curves and return rational vertices. Use the
- // base-class Extract if you want real-valued vertices.
- virtual void Extract(T level, std::vector<Vertex>& vertices,
- std::vector<Edge>& edges) override
- {
- // Adjust the image so that the level set is F(x,y) = 0.
- int64_t levelI64 = static_cast<int64_t>(level);
- for (size_t i = 0; i < this->mPixels.size(); ++i)
- {
- int64_t inputI64 = static_cast<int64_t>(this->mInputPixels[i]);
- this->mPixels[i] = inputI64 - levelI64;
- }
- vertices.clear();
- edges.clear();
- for (int y = 0, yp = 1; yp < this->mYBound; ++y, ++yp)
- {
- for (int x = 0, xp = 1; xp < this->mXBound; ++x, ++xp)
- {
- // Get the image values at the corners of the square.
- int i00 = x + this->mXBound * y;
- int i10 = i00 + 1;
- int i01 = i00 + this->mXBound;
- int i11 = i10 + this->mXBound;
- int64_t f00 = this->mPixels[i00];
- int64_t f10 = this->mPixels[i10];
- int64_t f01 = this->mPixels[i01];
- int64_t f11 = this->mPixels[i11];
- // Construct the vertices and edges of the level curve in
- // the square. The x, xp, y and yp values are implicitly
- // converted from int to int64_t (which is guaranteed to
- // be correct).
- ProcessSquare(vertices, edges, x, xp, y, yp, f00, f10, f11, f01);
- }
- }
- }
- protected:
- void ProcessSquare(std::vector<Vertex>& vertices, std::vector<Edge>& edges,
- int64_t x, int64_t xp, int64_t y, int64_t yp,
- int64_t f00, int64_t f10, int64_t f11, int64_t f01)
- {
- int64_t xn0, yn0, xn1, yn1, d0, d1, d2, d3, det;
- if (f00 != 0)
- {
- // convert to case "+***"
- if (f00 < 0)
- {
- f00 = -f00;
- f10 = -f10;
- f11 = -f11;
- f01 = -f01;
- }
- if (f10 > 0)
- {
- if (f11 > 0)
- {
- if (f01 > 0)
- {
- // ++++
- return;
- }
- else if (f01 < 0)
- {
- // +++-
- d0 = f11 - f01;
- xn0 = f11 * x - f01 * xp;
- d1 = f00 - f01;
- yn1 = f00 * yp - f01 * y;
- this->AddEdge(vertices, edges, xn0, d0, yp, 1, x, 1, yn1, d1);
- }
- else
- {
- // +++0
- this->AddVertex(vertices, x, 1, yp, 1);
- }
- }
- else if (f11 < 0)
- {
- d0 = f10 - f11;
- yn0 = f10 * yp - f11 * y;
- if (f01 > 0)
- {
- // ++-+
- d1 = f01 - f11;
- xn1 = f01 * xp - f11 * x;
- this->AddEdge(vertices, edges, xp, 1, yn0, d0, xn1, d1, yp, 1);
- }
- else if (f01 < 0)
- {
- // ++--
- d1 = f01 - f00;
- yn1 = f01 * y - f00 * yp;
- this->AddEdge(vertices, edges, x, 1, yn1, d1, xp, 1, yn0, d0);
- }
- else
- {
- // ++-0
- this->AddEdge(vertices, edges, x, 1, yp, 1, xp, 1, yn0, d0);
- }
- }
- else
- {
- if (f01 > 0)
- {
- // ++0+
- this->AddVertex(vertices, xp, 1, yp, 1);
- }
- else if (f01 < 0)
- {
- // ++0-
- d0 = f01 - f00;
- yn0 = f01 * y - f00 * yp;
- this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yn0, d0);
- }
- else
- {
- // ++00
- this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
- }
- }
- }
- else if (f10 < 0)
- {
- d0 = f00 - f10;
- xn0 = f00 * xp - f10 * x;
- if (f11 > 0)
- {
- d1 = f11 - f10;
- yn1 = f11 * y - f10 * yp;
- if (f01 > 0)
- {
- // +-++
- this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, yn1, d1);
- }
- else if (f01 < 0)
- {
- // +-+-
- d3 = f11 - f01;
- xn1 = f11 * x - f01 * xp;
- d2 = f01 - f00;
- yn0 = f01 * y - f00 * yp;
- if (d0*d3 > 0)
- {
- det = xn1 * d0 - xn0 * d3;
- }
- else
- {
- det = xn0 * d3 - xn1 * d0;
- }
- if (det > 0)
- {
- this->AddEdge(vertices, edges, xn1, d3, yp, 1, xp, 1, yn1, d1);
- this->AddEdge(vertices, edges, xn0, d0, y, 1, x, 1, yn0, d2);
- }
- else if (det < 0)
- {
- this->AddEdge(vertices, edges, xn1, d3, yp, 1, x, 1, yn0, d2);
- this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, yn1, d1);
- }
- else
- {
- this->AddEdge(vertices, edges, xn0, d0, yn0, d2, xn0, d0, y, 1);
- this->AddEdge(vertices, edges, xn0, d0, yn0, d2, xn0, d0, yp, 1);
- this->AddEdge(vertices, edges, xn0, d0, yn0, d2, x, 1, yn0, d2);
- this->AddEdge(vertices, edges, xn0, d0, yn0, d2, xp, 1, yn0, d2);
- }
- }
- else
- {
- // +-+0
- this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, yn1, d1);
- this->AddVertex(vertices, x, 1, yp, 1);
- }
- }
- else if (f11 < 0)
- {
- if (f01 > 0)
- {
- // +--+
- d1 = f11 - f01;
- xn1 = f11 * x - f01 * xp;
- this->AddEdge(vertices, edges, xn0, d0, y, 1, xn1, d1, yp, 1);
- }
- else if (f01 < 0)
- {
- // +---
- d1 = f01 - f00;
- yn1 = f01 * y - f00 * yp;
- this->AddEdge(vertices, edges, x, 1, yn1, d1, xn0, d0, y, 1);
- }
- else
- {
- // +--0
- this->AddEdge(vertices, edges, x, 1, yp, 1, xn0, d0, y, 1);
- }
- }
- else
- {
- if (f01 > 0)
- {
- // +-0+
- this->AddEdge(vertices, edges, xp, 1, yp, 1, xn0, d0, y, 1);
- }
- else if (f01 < 0)
- {
- // +-0-
- d1 = f01 - f00;
- yn1 = f01 * y - f00 * yp;
- this->AddEdge(vertices, edges, x, 1, yn1, d1, xn0, d0, y, 1);
- this->AddVertex(vertices, xp, 1, yp, 1);
- }
- else
- {
- // +-00
- this->AddEdge(vertices, edges, xp, 1, yp, 1, xn0, d0, yp, 1);
- this->AddEdge(vertices, edges, xn0, d0, yp, 1, x, 1, yp, 1);
- this->AddEdge(vertices, edges, xn0, d0, yp, 1, xn0, d0, y, 1);
- }
- }
- }
- else
- {
- if (f11 > 0)
- {
- if (f01 > 0)
- {
- // +0++
- this->AddVertex(vertices, xp, 1, y, 1);
- }
- else if (f01 < 0)
- {
- // +0+-
- d0 = f11 - f01;
- xn0 = f11 * x - f01 * xp;
- d1 = f00 - f01;
- yn1 = f00 * yp - f01 * y;
- this->AddEdge(vertices, edges, xn0, d0, yp, 1, x, 1, yn1, d1);
- this->AddVertex(vertices, xp, 1, y, 1);
- }
- else
- {
- // +0+0
- this->AddVertex(vertices, xp, 1, y, 1);
- this->AddVertex(vertices, x, 1, yp, 1);
- }
- }
- else if (f11 < 0)
- {
- if (f01 > 0)
- {
- // +0-+
- d0 = f11 - f01;
- xn0 = f11 * x - f01 * xp;
- this->AddEdge(vertices, edges, xp, 1, y, 1, xn0, d0, yp, 1);
- }
- else if (f01 < 0)
- {
- // +0--
- d0 = f01 - f00;
- yn0 = f01 * y - f00 * yp;
- this->AddEdge(vertices, edges, xp, 1, y, 1, x, 1, yn0, d0);
- }
- else
- {
- // +0-0
- this->AddEdge(vertices, edges, xp, 1, y, 1, x, 1, yp, 1);
- }
- }
- else
- {
- if (f01 > 0)
- {
- // +00+
- this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
- }
- else if (f01 < 0)
- {
- // +00-
- d0 = f00 - f01;
- yn0 = f00 * yp - f01 * y;
- this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yn0, d0);
- this->AddEdge(vertices, edges, xp, 1, yn0, d0, xp, 1, yp, 1);
- this->AddEdge(vertices, edges, xp, 1, yn0, d0, x, 1, yn0, d0);
- }
- else
- {
- // +000
- this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
- this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
- }
- }
- }
- }
- else if (f10 != 0)
- {
- // convert to case 0+**
- if (f10 < 0)
- {
- f10 = -f10;
- f11 = -f11;
- f01 = -f01;
- }
- if (f11 > 0)
- {
- if (f01 > 0)
- {
- // 0+++
- this->AddVertex(vertices, x, 1, y, 1);
- }
- else if (f01 < 0)
- {
- // 0++-
- d0 = f11 - f01;
- xn0 = f11 * x - f01 * xp;
- this->AddEdge(vertices, edges, x, 1, y, 1, xn0, d0, yp, 1);
- }
- else
- {
- // 0++0
- this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
- }
- }
- else if (f11 < 0)
- {
- if (f01 > 0)
- {
- // 0+-+
- d0 = f10 - f11;
- yn0 = f10 * yp - f11 * y;
- d1 = f01 - f11;
- xn1 = f01 * xp - f11 * x;
- this->AddEdge(vertices, edges, xp, 1, yn0, d0, xn1, d1, yp, 1);
- this->AddVertex(vertices, x, 1, y, 1);
- }
- else if (f01 < 0)
- {
- // 0+--
- d0 = f10 - f11;
- yn0 = f10 * yp - f11 * y;
- this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, yn0, d0);
- }
- else
- {
- // 0+-0
- d0 = f10 - f11;
- yn0 = f10 * yp - f11 * y;
- this->AddEdge(vertices, edges, x, 1, y, 1, x, 1, yn0, d0);
- this->AddEdge(vertices, edges, x, 1, yn0, d0, x, 1, yp, 1);
- this->AddEdge(vertices, edges, x, 1, yn0, d0, xp, 1, yn0, d0);
- }
- }
- else
- {
- if (f01 > 0)
- {
- // 0+0+
- this->AddVertex(vertices, x, 1, y, 1);
- this->AddVertex(vertices, xp, 1, yp, 1);
- }
- else if (f01 < 0)
- {
- // 0+0-
- this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, yp, 1);
- }
- else
- {
- // 0+00
- this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
- this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
- }
- }
- }
- else if (f11 != 0)
- {
- // convert to case 00+*
- if (f11 < 0)
- {
- f11 = -f11;
- f01 = -f01;
- }
- if (f01 > 0)
- {
- // 00++
- this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
- }
- else if (f01 < 0)
- {
- // 00+-
- d0 = f01 - f11;
- xn0 = f01 * xp - f11 * x;
- this->AddEdge(vertices, edges, x, 1, y, 1, xn0, d0, y, 1);
- this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, y, 1);
- this->AddEdge(vertices, edges, xn0, d0, y, 1, xn0, d0, yp, 1);
- }
- else
- {
- // 00+0
- this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
- this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
- }
- }
- else if (f01 != 0)
- {
- // cases 000+ or 000-
- this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
- this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
- }
- else
- {
- // case 0000
- this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
- this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
- this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
- this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
- }
- }
- };
- }
|