diff --git a/docs/docs/reference/wasm/modules/image/graph/_category_.json b/docs/docs/reference/wasm/modules/image/graph/_category_.json new file mode 100644 index 000000000..ae8de2258 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/_category_.json @@ -0,0 +1,9 @@ +{ + "label": "graph.h", + "link": { + "type": "generated-index", + "title": "Graph and Node Representation", + "description": "Documentation for Graph and Node Data Structure in the Image WebAssembly (WASM) module in Img2Num.", + "slug": "/reference/wasm/modules/image/graph" + } +} diff --git a/docs/docs/reference/wasm/modules/image/graph/api.md b/docs/docs/reference/wasm/modules/image/graph/api.md new file mode 100644 index 000000000..74e42fad8 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/api.md @@ -0,0 +1,298 @@ +--- +id: api +title: Graph / Node API +sidebar_label: API / Usage +sidebar_position: 5 +--- + +# Graph / Node API + +Each `Node` is a collection of pixels. A `Node` holds a `unique_ptr` to a vectors of pixels with `RGBXY` structure. \ +Each pixel has its own color and position. + +`Node`s reference neigbors through node pointers (`shared_ptr`) + +```cpp title="Nodes reference neigbors through node shared pointers" +Node_ptr n_ptr = std::make_shared(, > pixels>); +``` + +A `Graph` takes ownership over a collection of Nodes. It does so by referencing a list of Node pointers. + +```cpp title="A Graph takes ownership over a collection of Nodes. It does so by referencing a list of Node pointers." +std::unique_ptr> node_ptr = + std::make_unique>(std::move(nodes)); +Graph G(node_ptr, width, height); +``` + +:::tip +In sum, `Graph`s manage a list of Nodes through their pointers. Each `Node` can reference neighboring nodes as edges also through their pointers. +Since multiple entities can reference the same `Node` we use `shared_ptr`. +::: + +# Usage + +This follows the step-by-step guide in the [explanation](explained.md). + +1. Graph creation from kmeans labels + +- Initialize nodes and region map +- Use floodfill to fill out the region map and construct nodes + +```cpp +std::vector region_labels; +std::vector nodes; + +region_labeling(image_data, kmeans_labels, region_labels, width, height, nodes); +``` + +In `region_labeling` each Node is assigned an id and a collections of pixels: + +```cpp +Node_ptr n_ptr = std::make_shared(r_lbl, p_ptr); +nodes.push_back(n_ptr); +``` + +Then initialize the `Graph` + +```cpp +std::unique_ptr> node_ptr = + std::make_unique>(std::move(nodes)); +Graph G(node_ptr, width, height); +``` + +Finally add edges between `Nodes`: + +```cpp +G.discover_edges(region_labels, width, height); +``` + +2. Merge small regions/nodes + +```cpp +G.merge_small_area_nodes(min_area); +``` + +3. Compute contours and manage gaps + +```cpp +G.compute_contours(); +``` + +In this function nodes are iterated over one at a time. +Pseudocode: + +``` +for node in G.nodes +{ + // Consider all neigbors + for neigbor in node.edges + { + // collect pixels for each neighbor + } + /* + 1. Create joint grid plot of all pixels in node and neighbors + 2. Find edge pixels + 3. Decide if edge pixel should be added to the `node`'s or `neigbor`'s edge_pixel collection to ensure contour overlap + */ +} + +for node in G.nodes +{ + // compute contour per node +} +``` + +4. Collect all contours for SVG export + +--- + +# Node Class Documentation + +## Member Variables + +### Protected Members (Internal State) + +| Variable Name | Type | Description | +| :-------------- | :------------------------------------ | :----------------------------------------------------------------------------------------------------------------------- | +| `m_id` | `int32_t` | Unique identifier for the node. | +| `m_pixels` | `std::unique_ptr>` | Exclusive ownership of the raw pixel data defining this region. | +| `m_edges` | `std::set` | Adjacency list containing pointers to neighboring `Node` objects. | +| `m_edge_pixels` | `std::set` | Auxiliary pixels used for contour tracing. These are distinct from `m_pixels` and do not affect color/area calculations. | + +### Public Members + +| Variable Name | Type | Description | +| :------------ | :---------------- | :---------------------------------------------------------------------------------------------- | +| `m_contours` | `ColoredContours` | Vector representation of the node boundaries. Populated only after calling `compute_contour()`. | + +--- + +## API Reference + +### 1. Lifecycle + +#### `Node(int32_t id, std::unique_ptr> &pixels)` + +Constructs a new Node. + +- **id:** The unique integer ID. +- **pixels:** Reference to a unique pointer containing pixel data. Ownership is transferred to the Node using `std::move`. + +#### `void clear_all()` + +Resets the node completely, clearing pixel data, edges, and internal buffers. + +--- + +### 2. Geometric & Visual Properties + +#### `XY centroid() const` + +Calculates the geometric center of mass (average X, Y) of the region. + +#### `ImageLib::RGBPixel color() const` + +Computes the representative color of the node (typically the average color of all pixels in `m_pixels`). + +#### `std::array bounding_box_xywh() const` + +Calculates the axis-aligned bounding box. + +- **Returns:** `[min_x, min_y, width, height]` + +#### `size_t area() const` + +Returns the total number of pixels currently contained in the node. + +--- + +### 3. Graph Topology Management + +Methods to manage the adjacency list (`m_edges`). + +- `void add_edge(const Node_ptr &node)`: Adds a connection to a neighbor. +- `void remove_edge(const Node_ptr &node)`: Removes a specific connection. +- `void remove_all_edges()`: Clears all connections (isolates the node). +- `const std::set &edges() const`: Returns a read-only reference to the neighbor set. +- `size_t num_edges() const`: Returns the degree of the node. + +--- + +### 4. Image & Contour Operations + +#### `std::array create_binary_image(std::vector &binary) const` + +Rasterizes the node into a binary mask. + +- **binary:** Output buffer where the mask is written. +- **Returns:** Array describing dimensions/offsets of the generated mask. + +#### `void compute_contour()` + +Calculates the vector contours of the node based on edge pixels and populates `m_contours`. + +#### `void add_edge_pixel(const XY edge_pixel)` + +Adds a coordinate to the set used specifically for boundary tracing. + +#### `void clear_edge_pixels()` + +Clears the temporary edge pixel buffer. + +--- + +### 5. Data Access & Modification + +- `int32_t id() const`: Getter for the Node ID. +- `const std::vector &get_pixels() const`: Read-only access to the raw pixel vector. +- `ColoredContours &get_contours()`: Mutable access to the contour data. +- `void add_pixels(const std::vector &new_pixels)`: Merges new pixels into the existing node. + +--- + +# Graph Class Documentation + +## Member Variables + +### Protected Members (Internal State) + +| Variable Name | Type | Description | +| --------------------- | ---------------------------------------- | -------------------------------------------------------------------- | +| `m_width`, `m_height` | `int` | Dimensions of the original source image. | +| `m_nodes` | `std::unique_ptr>` | The collection of all nodes in the graph. | +| `m_node_ids` | `std::unordered_map` | A lookup map linking `Node ID` to `Vector Index` for fast retrieval. | + +--- + +### 1. Initialization + +#### `Graph(std::unique_ptr> &nodes, int width, int height)` + +Constructs the Graph. + +- **nodes:** A unique pointer to a vector of Node pointers. +- **Behavior:** The constructor calls `std::move` on the `nodes` argument, taking full ownership of the data. It also triggers `hash_node_ids()` to build the internal lookup map. + +--- + +### 2. Topology Analysis (Edge Discovery) + +#### `void discover_edges(const std::vector ®ion_labels, int32_t width, int32_t height)` + +Iterates through a raster label image to find adjacent regions. + +- **region_labels:** A flattened vector where each value represents the Node ID that the pixel belongs to. +- **Behavior:** Scans neighbors (8-connected) in the label map. If two adjacent pixels have different labels, an edge is added between the corresponding Nodes. + +#### `bool add_edge(int32_t node_id1, int32_t node_id2)` + +Manually creates a connection between two nodes identified by their IDs. +Calls `add_edge` for both Nodes + +- **Returns:** `true` if the edge was successfully added, `false` if nodes were not found. + +--- + +### 3. Graph Simplification (Merging & Pruning) + +#### `bool merge_nodes(const Node_ptr &node_to_keep, const Node_ptr &node_to_remove)` + +Combines two nodes into one. Called by `merge_small_area_nodes`. + +- **Behavior:** + +1. Transfers pixels and edges from `node_to_remove` to `node_to_keep`. +2. Updates the topology of neighbors. +3. Removes `node_to_remove` from the active graph. + +- **Returns:** `true` if merge successful. + +#### `void merge_small_area_nodes(int32_t min_area)` + +Iteratively merges nodes smaller than `min_area` into their largest neighbors. This is used to clean up "speckle" noise or insignificant regions. + +#### `void clear_unconnected_nodes()` + +Removes nodes that have no edges (orphaned regions) from the internal list. + +--- + +### 4. Data Processing & Access + +#### `void compute_contours()` + +Iterates through all nodes in the graph and triggers their individual `compute_contour()` methods. + +#### `const std::vector &get_nodes() const` + +Returns a read-only reference to the underlying vector of nodes. + +#### `size_t size()` + +Returns the number of nodes currently in the graph. + +#### `bool all_areas_bigger_than(int32_t min_area)` + +Utility check to verify if the graph simplification process (merging small nodes) is complete. + +- **Returns:** `true` if every node in the graph has an area greater than `min_area`. diff --git a/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide1.SVG b/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide1.SVG new file mode 100644 index 000000000..0d964a51c --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide1.SVG @@ -0,0 +1 @@ +0123456701234567Nodestd::vector<Node_ptr>Node_ptr= std::shared_ptr<Node>GraphGraph ownership viastd::unique_ptr<std::vector<Node_ptr>>01234567std::make_shared<Node>(id,pixels) \ No newline at end of file diff --git a/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide2.SVG b/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide2.SVG new file mode 100644 index 000000000..9dc8cc32e --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide2.SVG @@ -0,0 +1 @@ +0132456701231033012654ab \ No newline at end of file diff --git a/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide3.SVG b/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide3.SVG new file mode 100644 index 000000000..81366ea9c --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/diagrams/Slide3.SVG @@ -0,0 +1 @@ +01324567013245670132456713245671345672abcde \ No newline at end of file diff --git a/docs/docs/reference/wasm/modules/image/graph/explained.md b/docs/docs/reference/wasm/modules/image/graph/explained.md new file mode 100644 index 000000000..51e6366c2 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/explained.md @@ -0,0 +1,109 @@ +--- +id: explained +title: Graph and Node Explained +sidebar_position: 3 +--- + +To explain Graph creation from image regions we show a diagram representation followed by a concrete step-by-step example. + +## Mapping Images to Graphs + +Assume a region partition image as shown. Each partition (and the pixels contained within) are represented as a Node. Nodes exist in the heap. A Graph is a collection of shared pointers to these Nodes (`std::unique_ptr>>`). A Graph has unique ownership over these Nodes pointers. + +![slide1](./diagrams/Slide1.SVG) +_Figure 1: Partitioned Image_ + +Below is a visualization of the Graph with its connections forming a undirected graph (a). In reality each `Node` contains a list of pointers to their neighbors, represented as edges (b). Note that adjacent nodes point to each other. Edge management is handled by the Graph. +![slide2](./diagrams/Slide2.SVG) +_Figure 2: Visualizing a Graph. Note in (b) only 3 nodes are displayed 0,1,3. The remaining nodes have a similar structure_ + +## Small Area Node Merging (Small Region Pruning) + +One of the Graph's main operations is Node merging. Suppose we want to absorb/merge region 0 (Node 0) into region 2 (Node 2). Node 2 will assume ownership of Nodes 0's pixels and neighbors. The following figure shows the step-by-step process that happens in Graph. + +![slide3](./diagrams/Slide3.SVG) +_Figure 3: Step-by-step visualization of node merging_ + +**a**. Consider merging Node 0 into Node 2 + +**b**. Disconnect edges to Node 0 neighbors. This requires iterating over Nodes 1, 2, and 3, and removing Node 0 from their edge set. + +**c**. Transferring edges to absorbing node (Node 2). Again, iterate over Nodes 1 and 3 to assign Node 2 as an edge. Similarily Node 2 adds Node 1 and 3 as edges. In this case Node 2 and Node 3 already share an edge, but Node 1 gets a new edge. + +**d**. Pixels owned by Node 0 are passed to Node 2. Node 0 is removed from the graph. + +**e**. In image space, Node 2 (region 2) now contains the area that used to be Node 0's. + +## Example + +1. Consider starting from an image (after bilateral-filtering). + +![Image](./img/image.png) + +We use KMeans or a similar method to generate initial regions. But these regions do not have unique identifiers. \ +Notice also that there are many tiny regions caused by quantization noise even after filtering. + +![KMeans](./img/kmeans.png) + +We want a smart approach to merge small regions into neighboring regions. To do this build a `Graph`. + +2. Building graphs consists of 2 steps: + - Uniquely label neighboring regions. Floodfill is used for this. This discovers connected components - each connected component becomes a `Node` + + ![Regions](./img/regions.png) + - Tracking neighbors - the region labels are parsed to check for neighbors by looking at 3x3 neighborhoods. Each `Node` is updated with a list of neighbors forming an undirected graph. + + ![Graph](./img/region_graph.png) + +:::note Unexpected Centroid Locations +Since many regions are concave, their centroids may appear to be lying far away, but this is correct. Centroids here are just used for visualizing the graph. +::: + +3. Merging small regions + +The `Graph` is parsed to locate small area nodes. Each small node is merged with a larger neighbor node by transferring its pixels and edges to the absorbing neighbor node. + +![Graph2](./img/region_graph2.png) +![Regions2](./img/regions2.png) + +4. Contour creation + +Now that small nodes have been pruned, each `Node` can compute its own contour using the Suzuki-Abe method. + +![Contours](./img/contours.png) + +:::caution SVG space offset + +While Suzuki-Abe contour tracing correctly captures topological data, adjacent nodes end up with contours offset by ~1 pixel when mapped to SVG coordinates. + +
+See [this PR discussion](https://github.com/Ryan-Millard/Img2Num/pull/245#issuecomment-3807553622) for details on why this happens. +![ContoursZ](./img/contours_zoom.png) +
+ +::: + +5. Contour Management + +To solve this problem `Graph` has to overlap neighboring contours. For this `Nodes` have a special `edge_pixel` property to keep track of additional pixels to be considered for the contour. This forces neighboring contours to perfectly overlap creating no gaps or holes, which is important for SVGs. + +Before gap management + +
+Pull the code locally from PR #238 +```bash +# 1. Clone the repo (if you haven't already) +git clone https://github.com/Ryan-Millard/Img2Num.git +cd Img2Num +# 2. Fetch the specific commit from the PR +git fetch origin 9eb23f9a56edaeec95e2dfcfc8389b11bfd777b6 +# 3. Create a local branch pointing at it +git checkout -b try-pr-238 9eb23f9a56edaeec95e2dfcfc8389b11bfd777b6 +``` +
+ +![gaps](./img/cow_contours_gap.png) + +After gap management + +![nogaps](./img/cow_contours.png) diff --git a/docs/docs/reference/wasm/modules/image/graph/explained2.md b/docs/docs/reference/wasm/modules/image/graph/explained2.md new file mode 100644 index 000000000..091c77962 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/explained2.md @@ -0,0 +1,155 @@ +--- +id: explained2 +title: Important Functions Explained +sidebar_position: 4 +--- + +Nodes are responsible in computing their own contours using the Suzuki-Abe method. + +# Suzuki-Abe Contour Tracing Algorithm + +**Reference:** _Suzuki, S. and Abe, K., "Topological Structural Analysis of Digitized Binary Images by Border Following", CVGIP 30 1, pp 32-46 (1985)_ + +The Suzuki-Abe algorithm is a border-following technique designed to extract the topological structure of binary images. Unlike simple boundary tracing, it produces two distinct outputs: + +1. **Contours:** Vector lists of point coordinates representing the borders. +2. **Hierarchy:** A tree structure indicating the relationship between borders (e.g., is a border an external outline or a hole inside another shape?). + +--- + +## Key Concepts + +### 1. Border Types + +The algorithm distinguishes between two types of borders based on the transition between 0 (background) and 1 (foreground) pixels: + +- **Outer Border:** The boundary between a background region and a foreground component (surrounded by background). +- **Hole Border:** The boundary between a foreground component and an internal background hole (surrounded by foreground). + +### 2. Tracking Variables + +- **`f[i,j]`**: The value of the pixel at row `i`, column `j`. Initially, the image contains only `0` (background) and `1` (foreground). As the algorithm runs, `1`s are replaced by unique Border IDs. +- **`NBD` (Number of Border)**: A counter representing the current Border ID being assigned. Starts at `2` (since `1` is used for raw foreground). +- **`LNBD` (Last Number of Border)**: Tracks the Border ID of the most recently visited border during the raster scan. This acts as the "Parent" tracker. + +--- + +## The Algorithm Steps + +The algorithm combines a **Raster Scan** (to find new borders) with a **Border Following** routine (to trace and mark the entire boundary once found). + +### Phase 1: Raster Scan + +Iterate through the image row by row, from top-left to bottom-right. + +1. **Reset `LNBD`:** At the start of every row, reset `LNBD = 1` (frame ID). +2. **Check Pixel `f[i,j]`:** + - **Case A (Outer Border Start):** If `f[i,j] == 1` and `f[i, j-1] == 0` (transition from empty space to object): + - Increment `NBD`. + - This is a new **Outer Border**. + - Update Hierarchy: `LNBD` is the parent of `NBD`. + - **Trigger Phase 2 (Trace Border)** starting at `(i,j)`. + - **Case B (Hole Border Start):** If `f[i,j] >= 1` and `f[i, j+1] == 0` (transition from object to empty space): + - Increment `NBD`. + - This is a new **Hole Border**. + - Update Hierarchy: `LNBD` is the parent of `NBD`. (If `f[i,j] > 1`, set `LNBD = f[i,j]` first). + - **Trigger Phase 2 (Trace Border)** starting at `(i,j)`. + - **Case C (Non-Border):** If neither A nor B, update `LNBD` if `f[i,j] != 0` and continue scanning. + +### Phase 2: Border Following (The "Turtle" Logic) + +Once a starting pixel `(x, y)` is found, trace the connected edge pixels until returning to the start. + +1. **Search Neighborhood:** Starting from the previous pixel (or a default direction), check the 8-connected (or 4-connected) neighbors in a **Clockwise** (or Counter-Clockwise) direction. +2. **Find Next Pixel:** The first non-zero pixel found becomes the next current pixel. +3. **Mark Pixel:** Change the value of the current pixel in the image to `NBD` (or `-NBD` in specific cases to mark visited edges without destroying topology). + - _Note: This modification prevents the Raster Scan from re-detecting the same border later._ +4. **Record Coordinate:** Add the pixel `(x, y)` to the current contour vector. +5. **Termination:** Stop when the tracer returns to the **Starting Point** AND matches the **Starting Direction**. + +--- + +## Hierarchy Tree Structure + +The algorithm maintains a hierarchy table (often stored as an `std::array` : `[Next, Previous, First_Child, Parent]`). + +| Topology | Description | +| :--------- | :----------------------------------------------------- | +| **Root** | The image frame (background). | +| **Parent** | The border immediately surrounding the current one. | +| **Child** | A border immediately contained within the current one. | + +**Example Hierarchy:** + +1. **Outer Box (ID 2)** -> Parent: Frame. +2. **Inner Hole (ID 3)** -> Parent: ID 2. +3. **Island inside Hole (ID 4)** -> Parent: ID 3. + +--- + +# Contour Data Structures + +This file documents the core data structures used to represent vector boundaries, Bezier curves, and the topological hierarchy of image regions. + +## 1. Geometric Primitives + +### `struct Point` + +A fundamental 2D coordinate representing a position in the image space. + +| Member | Type | Description | +| :----- | :------ | :-------------------------------------------------------------------------------------------------------- | +| `x` | `float` | Horizontal coordinate (column). Uses `float` to support sub-pixel precision or smooth vector coordinates. | +| `y` | `float` | Vertical coordinate (row). | + +### `struct QuadBezier` + +Represents a Quadratic Bezier curve segment. This is used when the raw pixel contours are approximated or smoothed into vector paths. + +| Member | Type | Description | +| :----- | :------ | :--------------------------------------------------------------------------------------------------------------------------------- | +| `p0` | `Point` | **Start Point:** The anchor point where the curve begins. | +| `p1` | `Point` | **Control Point:** The handle that determines the curve's tangent and shape. The curve generally does not pass through this point. | +| `p2` | `Point` | **End Point:** The anchor point where the curve ends. | + +--- + +## 2. Contour Results + +### `struct ContoursResult` + +The primary container for the output of the contour extraction algorithm (e.g., Suzuki-Abe). It separates the raw pixel data from the topological relationship data. + +#### Member Variables + +| Member Variable | Type | Description | +| :-------------- | :--------------------------- | :-------------------------------------------------------------------------------------------------------------------- | +| `contours` | `vector>` | A list of contours. `contours[k]` is a vector of `Point`s tracing the boundary of the $k$-th region. | +| `curves` | `vector>` | A vectorized representation of `contours`. `curves[k]` contains the Bezier segments approximating the $k$-th contour. | +| `hierarchy` | `vector>` | Topological tree structure describing how contours nest within each other (see **Hierarchy Structure** below). | +| `is_hole` | `vector` | Flags the type of border. `true` if `contours[k]` is an internal hole; `false` if it is an external boundary. | + +#### Hierarchy Structure (`std::array`) + +The `hierarchy` vector follows the standard structure (compatible with OpenCV) to represent the nesting tree. For the $k$-th contour, `hierarchy[k]` contains: + +| Index | Name | Description | +| :---- | :--------------- | :------------------------------------------------------------------------------------ | +| `0` | **Next Sibling** | Index of the next contour at the same tree level. `-1` if none. | +| `1` | **Prev Sibling** | Index of the previous contour at the same tree level. `-1` if none. | +| `2` | **First Child** | Index of the first contour nested _inside_ this contour. `-1` if no children. | +| `3` | **Parent** | Index of the contour that surrounds this contour. `-1` if it is a root/frame contour. | + +--- + +## 3. Visual Extension + +### `struct ColoredContours` + +**Inherits from:** `ContoursResult` + +Extends the geometric data with visual attributes, specifically assigning a color to each contour. This is useful for visualization or when contours inherit the color properties of the underlying image nodes. + +| Member Variable | Type | Description | +| :-------------- | :--------------------------- | :---------------------------------------------------------------------------------------------------- | +| `colors` | `vector>` | A parallel vector to `contours`. `colors[k]` holds the RGBA color associated with the $k$-th contour. | diff --git a/docs/docs/reference/wasm/modules/image/graph/img/contours.png b/docs/docs/reference/wasm/modules/image/graph/img/contours.png new file mode 100644 index 000000000..97950841f Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/contours.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/contours_zoom.png b/docs/docs/reference/wasm/modules/image/graph/img/contours_zoom.png new file mode 100644 index 000000000..de1665bb8 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/contours_zoom.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/cow_contours.png b/docs/docs/reference/wasm/modules/image/graph/img/cow_contours.png new file mode 100644 index 000000000..d262d5705 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/cow_contours.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/cow_contours_gap.png b/docs/docs/reference/wasm/modules/image/graph/img/cow_contours_gap.png new file mode 100644 index 000000000..6951ff4ce Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/cow_contours_gap.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/image.png b/docs/docs/reference/wasm/modules/image/graph/img/image.png new file mode 100644 index 000000000..948a3f16e Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/image.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/kmeans.png b/docs/docs/reference/wasm/modules/image/graph/img/kmeans.png new file mode 100644 index 000000000..53ddb39e0 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/kmeans.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/region_graph.png b/docs/docs/reference/wasm/modules/image/graph/img/region_graph.png new file mode 100644 index 000000000..a224d52b3 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/region_graph.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/region_graph2.png b/docs/docs/reference/wasm/modules/image/graph/img/region_graph2.png new file mode 100644 index 000000000..ea8c4d1a7 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/region_graph2.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/regions.png b/docs/docs/reference/wasm/modules/image/graph/img/regions.png new file mode 100644 index 000000000..10cf09b86 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/regions.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/img/regions2.png b/docs/docs/reference/wasm/modules/image/graph/img/regions2.png new file mode 100644 index 000000000..9991c7678 Binary files /dev/null and b/docs/docs/reference/wasm/modules/image/graph/img/regions2.png differ diff --git a/docs/docs/reference/wasm/modules/image/graph/overview.md b/docs/docs/reference/wasm/modules/image/graph/overview.md new file mode 100644 index 000000000..cae1a3722 --- /dev/null +++ b/docs/docs/reference/wasm/modules/image/graph/overview.md @@ -0,0 +1,11 @@ +--- +id: overview +title: Graph Data Structure Overview +sidebar_label: Overview +sidebar_position: 2 +--- + +# Graph Data Structure Overview + +Images are pixels on a 2D grid. This is a type of densely connected graph, where each pixel in a quantized local region is a node and each neighbor is an edge. We can break images into regions using KMeans or SLIC. In this case we want to track neighboring relations between these regions. We convert images into a `Graph` datastructure consisting of `Node` nodes. Each `Node` is a collection of pixels representing a unique region. Each `Node` tracks its immediate neighbors. The `Graph` manages all `Node`s. +![illustration](diagrams/Slide1.SVG) diff --git a/src/components/WasmImageProcessor.jsx b/src/components/WasmImageProcessor.jsx index 52bee44c8..e05f7fe9a 100644 --- a/src/components/WasmImageProcessor.jsx +++ b/src/components/WasmImageProcessor.jsx @@ -103,7 +103,7 @@ const WasmImageProcessor = () => { step(95); const { svg } = await findContours({ - pixels: kmeansed, + pixels: imgBilateralFiltered, labels, width, height, diff --git a/src/hooks/useWasmWorker.js b/src/hooks/useWasmWorker.js index 4fd718e4f..428396436 100644 --- a/src/hooks/useWasmWorker.js +++ b/src/hooks/useWasmWorker.js @@ -65,7 +65,7 @@ export function useWasmWorker() { sigma_spatial = 3, sigma_range = 50, color_space = 0, - n_threads = 1, // Default because headers that enable threads aren't currently supported on GH Pages + n_threads = 8, }) => { const result = await call({ funcName: 'bilateral_filter', @@ -93,7 +93,7 @@ export function useWasmWorker() { num_colors, max_iter = 100, color_space = 0, - n_threads = 1, // Default because headers that enable threads aren't currently supported on GH Pages + n_threads = 8, }) => { const result = await call({ funcName: 'kmeans', diff --git a/src/wasm/modules/image/CMakeLists.txt b/src/wasm/modules/image/CMakeLists.txt index 4a32ab17c..286238021 100644 --- a/src/wasm/modules/image/CMakeLists.txt +++ b/src/wasm/modules/image/CMakeLists.txt @@ -30,6 +30,8 @@ target_include_directories(${MODULE_NAME}_wasm PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include ) +set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lidbfs.js") + # Shared Emscripten options set(COMMON_FLAGS "SHELL:-s MODULARIZE=1" diff --git a/src/wasm/modules/image/include/LABPixel.h b/src/wasm/modules/image/include/LABPixel.h index b16ac913d..d8552f1e9 100644 --- a/src/wasm/modules/image/include/LABPixel.h +++ b/src/wasm/modules/image/include/LABPixel.h @@ -2,6 +2,7 @@ #define LABPIXEL_H #include "Pixel.h" +#include /* can support signed data types @@ -30,6 +31,17 @@ template struct LABPixel : public Pixel { a = b = 0; } + static inline float colorDistance(const LABPixel &a, + const LABPixel &b) { + + LABPixel af{static_cast(a.l), static_cast(a.a), + static_cast(a.b)}; + LABPixel bf{static_cast(b.l), static_cast(b.a), + static_cast(b.b)}; + return std::sqrt((a.l - b.l) * (a.l - b.l) + (a.a - b.a) * (a.a - b.a) + + (a.b - b.b) * (a.b - b.b)); + } + } __attribute__((packed)); } // namespace ImageLib diff --git a/src/wasm/modules/image/include/Point.h b/src/wasm/modules/image/include/Point.h new file mode 100644 index 000000000..2526e3e7c --- /dev/null +++ b/src/wasm/modules/image/include/Point.h @@ -0,0 +1,41 @@ +#ifndef POINT_H +#define POINT_H + +// will start as integer values but can be adjusted to subpixel positions +struct Point { + float x = 0; + float y = 0; + + Point operator+(const Point &other) const { + return Point{x + other.x, y + other.y}; + } + + // Overload the - operator (subtraction) + Point operator-(const Point &other) const { + return Point{x - other.x, y - other.y}; + } + + // Overload the * operator (multiplication by a scalar) + // The left operand is the class object (this), the right is a double. + Point operator*(float scalar) const { return Point{x * scalar, y * scalar}; } + + Point operator/(float scalar) const { return Point{x / scalar, y / scalar}; } + + // Friend function to allow float * Vector + friend Point operator*(float scalar, const Point &v) { + return Point{v.x * scalar, v.y * scalar}; + ; // Calls the member function for the actual logic + } + + Point &operator+=(const Point &other) { + x += other.x; + y += other.y; + return *this; + } + + static float distSq(Point a, Point b) { + return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y); + } +}; + +#endif \ No newline at end of file diff --git a/src/wasm/modules/image/include/RGBPixel.h b/src/wasm/modules/image/include/RGBPixel.h index 7f8423737..3fee0126e 100644 --- a/src/wasm/modules/image/include/RGBPixel.h +++ b/src/wasm/modules/image/include/RGBPixel.h @@ -2,6 +2,7 @@ #define RGBPIXEL_H #include "Pixel.h" +#include namespace ImageLib { template struct RGBPixel : public Pixel { @@ -21,6 +22,19 @@ template struct RGBPixel : public Pixel { // ----- Utilities ----- inline void setGray(NumberT gray) { red = green = blue = gray; } + + static inline float colorDistance(const RGBPixel &a, + const RGBPixel &b) { + + RGBPixel af{static_cast(a.red), static_cast(a.green), + static_cast(a.blue)}; + RGBPixel bf{static_cast(b.red), static_cast(b.green), + static_cast(b.blue)}; + return std::sqrt((af.red - bf.red) * (af.red - bf.red) + + (af.green - bf.green) * (af.green - bf.green) + + (af.blue - bf.blue) * (af.blue - bf.blue)); + } + } __attribute__((packed)); } // namespace ImageLib diff --git a/src/wasm/modules/image/include/SavitskyGolay.h b/src/wasm/modules/image/include/SavitskyGolay.h new file mode 100644 index 000000000..7f940447a --- /dev/null +++ b/src/wasm/modules/image/include/SavitskyGolay.h @@ -0,0 +1,31 @@ +#ifndef SOLVER_H +#define SOLVER_H + +#include "Point.h" +#include +#include +#include +#include +#include + +struct Point; + +class SavitzkyGolay { +private: + int window_size_; + int m_; // half window size + int poly_order_; + std::vector coeffs_; + + std::vector> + invert_matrix(std::vector> A); + void compute_coefficients(); + +public: + SavitzkyGolay(int radius, int poly_order); + std::vector filter(const std::vector &data); + std::vector filter_wrap(const std::vector &data); + std::vector get_coeffs() const { return coeffs_; } +}; + +#endif \ No newline at end of file diff --git a/src/wasm/modules/image/include/bezier.h b/src/wasm/modules/image/include/bezier.h new file mode 100644 index 000000000..b35937143 --- /dev/null +++ b/src/wasm/modules/image/include/bezier.h @@ -0,0 +1,9 @@ +#ifndef BEZIER_H +#define BEZIER_H + +#include "contours.h" + +void fit_curve_reduction(const std::vector> &chains, + std::vector> &results, + float tolerance); +#endif \ No newline at end of file diff --git a/src/wasm/modules/image/include/contours.h b/src/wasm/modules/image/include/contours.h index 305011a23..3a575ab4c 100644 --- a/src/wasm/modules/image/include/contours.h +++ b/src/wasm/modules/image/include/contours.h @@ -1,22 +1,32 @@ #ifndef CONTOURS_H #define CONTOURS_H +#include "Image.h" +#include "PixelConverters.h" +#include "Point.h" +#include "RGBAPixel.h" +#include "SavitskyGolay.h" #include -#include #include #include #include #include -struct Point { - int x = 0; - int y = 0; +struct QuadBezier { + Point p0{0, 0}; // Start + Point p1{0, 0}; // Control + Point p2{0, 0}; // End +}; + +struct Rect { + float x, y, width, height; }; struct ContoursResult { // contours[k] is a sequence of boundary pixels (x,y) in image coordinates // (0..w-1, 0..h-1) std::vector> contours; + std::vector> curves; // hierarchy[k] = { next_sibling, prev_sibling, first_child, parent } // -1 means "none" @@ -26,8 +36,22 @@ struct ContoursResult { std::vector is_hole; }; +struct ColoredContours : ContoursResult { + // inherits: contours, hierarchy, is_hole + + std::vector> colors; +}; + namespace contours { ContoursResult find_contours(const std::vector &binary, int width, int height); -} + +void stitch_smooth(std::vector &vecA, std::vector &vecB); +void coupled_smooth(std::vector> &contours, Rect bounds); + +void pack_with_boundary_constraints(std::vector> &contours, + Rect bounds, int iterations = 15); + +} // namespace contours + #endif \ No newline at end of file diff --git a/src/wasm/modules/image/include/graph.h b/src/wasm/modules/image/include/graph.h index 23ab345cc..c3bdc0ae7 100644 --- a/src/wasm/modules/image/include/graph.h +++ b/src/wasm/modules/image/include/graph.h @@ -31,14 +31,16 @@ discover_edges(G, region_labels, width, height); class Graph { protected: + int m_width, m_height; std::unique_ptr> m_nodes; std::unordered_map m_node_ids; void hash_node_ids(void); public: - inline Graph(std::unique_ptr> &nodes) - : m_nodes(std::move(nodes)) { + inline Graph(std::unique_ptr> &nodes, int width, + int height) + : m_nodes(std::move(nodes)), m_width(width), m_height(height) { hash_node_ids(); } @@ -56,6 +58,7 @@ class Graph { void discover_edges(const std::vector ®ion_labels, const int32_t width, const int32_t height); void merge_small_area_nodes(const int32_t min_area); + void compute_contours(); }; #endif diff --git a/src/wasm/modules/image/include/node.h b/src/wasm/modules/image/include/node.h index ece8dfb68..286fd8d41 100644 --- a/src/wasm/modules/image/include/node.h +++ b/src/wasm/modules/image/include/node.h @@ -2,6 +2,7 @@ #define NODE_H #include "RGBPixel.h" +#include "contours.h" #include #include #include @@ -35,6 +36,9 @@ discover_edges(G, region_labels, width, height); struct XY { int32_t x, y; + std::pair xy; + XY(int32_t x_, int32_t y_) : x(x_), y(y_) { xy = std::make_pair(x, y); }; + bool operator<(const XY &rhs) const { return xy < rhs.xy; }; }; struct RGBXY { @@ -59,6 +63,10 @@ class Node { std::unique_ptr> m_pixels; std::set m_edges{}; + // pixels considered for contour tracing but not influencing other + // node properties such as color + std::set m_edge_pixels{}; + public: inline Node(int32_t id, std::unique_ptr> &pixels) : m_id(id), m_pixels(std::move(pixels)) {} @@ -68,15 +76,25 @@ class Node { std::array bounding_box_xywh() const; std::array create_binary_image(std::vector &binary) const; + // keep track of its own contour points + // only filled in when compute_contour() is called + // though these are public only Graph should access them + ColoredContours m_contours; + void clear_contour(); + void compute_contour(); + /* access member variables */ inline int32_t id() const { return m_id; }; inline size_t area() const { return m_pixels->size(); }; inline const std::set &edges() const { return m_edges; } inline size_t num_edges() const { return m_edges.size(); } inline const std::vector &get_pixels() const { return *m_pixels; } + inline ColoredContours &get_contours() { return m_contours; } /* modify member variables */ void add_pixels(const std::vector &new_pixels); + void add_edge_pixel(const XY edge_pixel); + void clear_edge_pixels(); void clear_all(); diff --git a/src/wasm/modules/image/include/utils.h b/src/wasm/modules/image/include/utils.h new file mode 100644 index 000000000..2fe8add18 --- /dev/null +++ b/src/wasm/modules/image/include/utils.h @@ -0,0 +1,103 @@ +#ifndef UTILS_H +#define UTILS_H + +#include +#include +#include +#include +#include +#include + +void EMSCRIPTEN_KEEPALIVE +saveSVG(const std::string &filename, int width, int height, + const std::vector> &rawChains, + const std::vector> &closedLoops) { + + EM_ASM({ + // Make a directory and mount IDBFS + FS.mkdir('/offline'); + FS.mount(IDBFS, {}, '/offline'); + // Then sync to load existing data and prepare for saving + FS.syncfs( + true, function(err) { + if (err) + console.log("Sync error: " + err); + else + console.log("IDBFS synced."); + }); + }); + + std::ofstream f(std::strcat("/offline/", filename.c_str())); + if (!f.is_open()) { + std::cerr << "Error opening file: " << filename << "\n"; + return; + } + + // SVG Header + f << "\n"; + + // Background (White) + f << "\n"; + + // Style Definitions + f << "\n"; + + // 1. Draw Raw Chains (Background Layer - Light Gray) + for (const auto &chain : rawChains) { + f << "\n"; + } + + // 2. Draw Closed Loops (Foreground Layer - Colors) + std::string colors[] = {"#FF0000", "#00AA00", "#0000FF", + "#FF00FF", "#FFAA00", "#00AAAA"}; + int colorIdx = 0; + + for (const auto &loop : closedLoops) { + std::string color = colors[colorIdx % 6]; + colorIdx++; + + f << "\n"; + } + + f << "\n"; + f.close(); + std::cout << "Saved visualization to " << filename << "\n"; + + // Sync again to save changes from MEMFS to IndexedDB + EM_ASM({ + FS.syncfs( + false, function(err) { + if (err) + console.log("Sync error: " + err); + else + console.log("File saved to IndexedDB."); + }); + }); + // 2. Trigger browser download using JS + /*EM_ASM({ + var a = document.createElement('a'); + a.href = 'output.svg'; + a.download = 'output.svg'; + document.body.appendChild(a); + a.click(); + document.body.removeChild(a); + });*/ +} + +#endif \ No newline at end of file diff --git a/src/wasm/modules/image/src/SavitskyGolay.cpp b/src/wasm/modules/image/src/SavitskyGolay.cpp new file mode 100644 index 000000000..f9ded82d7 --- /dev/null +++ b/src/wasm/modules/image/src/SavitskyGolay.cpp @@ -0,0 +1,137 @@ +#include "SavitskyGolay.h" +#include + +SavitzkyGolay::SavitzkyGolay(int radius, int poly_order) + : m_(radius), window_size_(2 * radius + 1), poly_order_(poly_order) { + + assert(radius >= 0); + assert(window_size_ > poly_order_); + + compute_coefficients(); +} + +std::vector SavitzkyGolay::filter(const std::vector &data) { + if (data.size() < window_size_) { + return data; // Data too short to filter + } + + std::vector result(data.size()); + + // 1. Convolution for the valid range + for (size_t i = m_; i < data.size() - m_; ++i) { + Point val{0.0, 0.0}; + for (int j = -m_; j <= m_; ++j) { + val += data[i + j] * coeffs_[j + m_]; + } + result[i] = val; + } + + // 2. Handle Edges (Simple Repeat/Nearest padding strategy) + // For a robust production app, you might calculate asymmetric kernels here. + for (int i = 0; i < m_; ++i) + result[i] = data[i]; + for (size_t i = data.size() - m_; i < data.size(); ++i) + result[i] = data[i]; + + return result; +} + +std::vector SavitzkyGolay::filter_wrap(const std::vector &data) { + // wrap around + if (data.size() < window_size_) { + return data; // Data too short to filter + } + + std::vector result(data.size()); + std::copy(data.begin(), data.end(), result.begin()); + + for (size_t i = 0; i < data.size(); ++i) { + Point val{0.0, 0.0}; + for (int j = -m_; j <= m_; ++j) { + int k = i + j; + if (k < 0) { + k = data.size() + k; + } else if (k >= data.size()) { + k = k - data.size(); + } + val = val + data[k] * coeffs_[j + m_]; + } + result[i] = val; + } + + return result; +} + +// Helper: Invert a matrix using Gauss-Jordan Elimination +// A is (N x N), returns A_inv +std::vector> +SavitzkyGolay::invert_matrix(std::vector> A) { + int n = A.size(); + std::vector> inv(n, std::vector(n, 0.0)); + + // Initialize inverse as identity + for (int i = 0; i < n; ++i) + inv[i][i] = 1.0; + + for (int i = 0; i < n; ++i) { + // Find pivot + float pivot = A[i][i]; + // (Simple pivot check, typically you'd swap rows for stability) + if (std::abs(pivot) < 1e-10) + throw std::runtime_error("Matrix singular, cannot invert."); + + // Normalize row + for (int j = 0; j < n; ++j) { + A[i][j] /= pivot; + inv[i][j] /= pivot; + } + + // Eliminate other rows + for (int k = 0; k < n; ++k) { + if (k != i) { + float factor = A[k][i]; + for (int j = 0; j < n; ++j) { + A[k][j] -= factor * A[i][j]; + inv[k][j] -= factor * inv[i][j]; + } + } + } + } + return inv; +} + +void SavitzkyGolay::compute_coefficients() { + // 1. Create the matrix J = (A^T * A) + // Size is (poly_order + 1) x (poly_order + 1) + // Element J[i][j] is the sum of k^(i+j) for k in -m..m + + int rows = poly_order_ + 1; + std::vector> J(rows, std::vector(rows)); + + for (int i = 0; i < rows; ++i) { + for (int j = 0; j < rows; ++j) { + float sum = 0; + for (int k = -m_; k <= m_; ++k) { + sum += std::pow(k, i + j); + } + J[i][j] = sum; + } + } + + // 2. Invert J to solve the normal equations + auto J_inv = invert_matrix(J); + + // 3. Compute the weights + // The smoothed value is the coefficient c_0 of the polynomial. + // c_0 = sum( weight_k * y_k ) + // weight_k = sum( J_inv[0][j] * k^j ) for j=0..order + + coeffs_.resize(window_size_); + for (int k = -m_; k <= m_; ++k) { + float weight = 0.0; + for (int j = 0; j < rows; ++j) { + weight += J_inv[0][j] * std::pow(k, j); + } + coeffs_[k + m_] = weight; + } +} \ No newline at end of file diff --git a/src/wasm/modules/image/src/bezier.cpp b/src/wasm/modules/image/src/bezier.cpp new file mode 100644 index 000000000..67cc74d83 --- /dev/null +++ b/src/wasm/modules/image/src/bezier.cpp @@ -0,0 +1,153 @@ +#include "bezier.h" +#include +#include +#include + +// --- Vector Math Helpers --- +inline float dot(Point a, Point b) { return a.x * b.x + a.y * b.y; } +inline float len(Point a, Point b) { + Point c = a - b; + return std::sqrt(c.x * c.x + c.y * c.y); +} + +// --- Evaluate Quadratic Bezier at t --- +Point evalBezier(const QuadBezier &b, float t) { + // B(t) = (1-t)^2 * P0 + 2*t*(1-t) * P1 + t^2 * P2 + float invT = 1.0 - t; + float c0 = invT * invT; + float c1 = 2.0f * t * invT; + float c2 = t * t; + + return c0 * b.p0 + c1 * b.p1 + c2 * b.p2; +} + +// --- Chord Length Parameterization --- +// Assigns a 't' value (0.0 to 1.0) to each point based on distance +std::vector chordLengthParameterize(const std::vector &points) { + std::vector u; + u.reserve(points.size()); + u.push_back(0.0f); + + for (int i = 1; i < points.size(); ++i) { + float dist = len(points[i], points[i - 1]); + u.push_back(u.back() + dist); + } + + float totalLen = u.back(); + if (totalLen == 0) + return u; // Should not happen for valid ranges + + for (float &val : u) { + val /= totalLen; + } + return u; +} + +// --- Least Squares Fit for Control Point Q1 --- +// We know Q0 (Start) and Q2 (End). We need to find Q1 that minimizes error. +// Based on equation: P(t) = (1-t)^2 Q0 + 2t(1-t) Q1 + t^2 Q2 +// Rearranged: Q1 * [2t(1-t)] = P(t) - (1-t)^2 Q0 - t^2 Q2 +Point generateQuadBezier(const std::vector &points, + const std::vector &u) { + Point Q0 = points.front(); + Point Q2 = points.back(); + + float numX = 0.0, numY = 0.0; + float den = 0.0; + + for (int i = 0; i < u.size(); ++i) { + float t = u[i]; + float invT = 1.0f - t; + + // A = 2t(1-t) + float A = 2.0f * t * invT; + + // V = P_actual - (Contribution of Q0 and Q2) + // V = P[i] - (1-t)^2 * Q0 - t^2 * Q2 + float B0 = invT * invT; + float B2 = t * t; + Point V = points[i] - (Q0 * B0 + Q2 * B2); + + // Least Squares Sums + numX += A * V.x; + numY += A * V.y; + den += A * A; + } + + if (den < 1e-9) { + // Fallback for straight lines (den is 0 if all t are 0 or 1) + return Q0 + (Q2 - Q0) * 0.5; + } + + return {numX / den, numY / den}; +} + +// --- Recursive Fit Function --- +void fitRecursive(const std::vector &points, float errorLimit, + std::vector &outCurves) { + + // Base Case: Not enough points, just connect them + if (points.size() <= 2) { + // Just a line segment + Point mid = points.front() + (points.back() - points.front()) * 0.5; + outCurves.push_back({points.front(), mid, points.back()}); + return; + } + + // 1. Parameterize Points + std::vector u = chordLengthParameterize(points); + + // 2. Find Optimal Control Point (Q1) + Point Q1 = generateQuadBezier(points, u); + QuadBezier curve = {points.front(), Q1, points.back()}; + + // 3. Calculate Maximum Error + float maxDistSq = 0.0f; + int splitPoint = 0; + + // Check distance of every intermediate point to the curve + // Note: Technically we should find the nearest point on curve, + // but evaluating at parameter 't' is a standard approximation for speed. + for (int i = 0; i < u.size(); ++i) { + Point P = points[i]; + Point CurveP = evalBezier(curve, u[i]); + float d2 = Point::distSq(P, CurveP); + + if (d2 > maxDistSq) { + maxDistSq = d2; + splitPoint = i; + } + } + + // 4. Check Error Threshold + if (maxDistSq < (errorLimit * errorLimit)) { + outCurves.push_back(curve); // Fit is good! + } else { + // Fit is bad, split at the point of maximum error + // Important: Prevent infinite recursion if split doesn't advance + if (splitPoint == 0 || splitPoint == points.size() - 1) { + // Fallback: simply bisect indices if geometric split fails + splitPoint = points.size() / 2; + } + + std::vector p1(points.begin(), points.begin() + splitPoint + 1); + std::vector p2(points.begin() + splitPoint, points.end()); + + fitRecursive(p1, errorLimit, outCurves); + fitRecursive(p2, errorLimit, outCurves); + } +} + +// --- Main Wrapper --- +void fit_curve_reduction(const std::vector> &chains, + std::vector> &results, + float tolerance) { + // if (chain.empty()) return result; + // results.resize(chains.size()); + for (int i = 0; i < chains.size(); ++i) { + // Start recursion on the whole chain + std::vector result; + fitRecursive(chains[i], tolerance, result); + results.push_back(result); + } +} \ No newline at end of file diff --git a/src/wasm/modules/image/src/contours.cpp b/src/wasm/modules/image/src/contours.cpp index a5dfe0ec3..035839bc0 100644 --- a/src/wasm/modules/image/src/contours.cpp +++ b/src/wasm/modules/image/src/contours.cpp @@ -1,5 +1,11 @@ #include "contours.h" +#include +#include +#include +#include +#include +#include namespace contours { @@ -71,7 +77,8 @@ static std::vector traceBorder(std::vector &f, int paddedW, int sy, // If isolated pixel: set f(sy,sx) = -NBD and return single-point contour if (!found) { set(sy, sx) = -nbd; - pts.push_back(Point{sx - 1, sy - 1}); + pts.push_back( + Point{static_cast(sx - 1), static_cast(sy - 1)}); return pts; } @@ -119,7 +126,8 @@ static std::vector traceBorder(std::vector &f, int paddedW, int sy, } // Record current point in unpadded coordinates - pts.push_back(Point{x3 - 1, y3 - 1}); + pts.push_back( + Point{static_cast(x3 - 1), static_cast(y3 - 1)}); // (3.5) Termination check if (y4 == sy && x4 == sx && y3 == firstY && x3 == firstX) { @@ -285,4 +293,619 @@ ContoursResult find_contours(const std::vector &binary, int width, return out; } -} // namespace contours +// Helper: 2D integer coordinate for map keys +struct Coord { + int x, y; + bool operator<(const Coord &other) const { + return std::tie(x, y) < std::tie(other.x, other.y); + } +}; + +// Helper: Normalize a vector +Point normalize(Point p) { + float len = std::sqrt(p.x * p.x + p.y * p.y); + if (len == 0) + return {0, 0}; + return {p.x / len, p.y / len}; +} + +// Helper: Calculate Tangent of A at index i using neighbors +Point getTangent(const std::vector &vec, int i) { + int n = vec.size(); + if (n < 2) + return {0, 0}; + + // Use previous and next points to determine the "flow" of the line + int prev = (i == 0) ? 0 : i - 1; + int next = (i == n - 1) ? n - 1 : i + 1; + + return normalize(vec[next] - vec[prev]); +} + +// Calculate the closest point on the segment V -> W from point P +// Returns {ClosestPoint, DistanceSquared} +std::pair getClosestPointOnSegment(Point p, Point v, Point w) { + float l2 = Point::distSq(v, w); + if (l2 == 0.0) + return {v, Point::distSq(p, v)}; + + // Project p onto line v-w + // t is the parameterized distance along the line (0.0 to 1.0) + float t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2; + + // Clamp to segment + t = std::max(0.0f, std::min(1.0f, t)); + + Point projection = {v.x + t * (w.x - v.x), v.y + t * (w.y - v.y)}; + return {projection, Point::distSq(p, projection)}; +} + +// Identifies a specific point: {ContourIndex, PointIndex} +struct PointID { + int cIdx; + int pIdx; +}; + +/** + * @brief Computes a smoothed target point using a local quadratic + * Savitzky-Golay filter. + * + * This function smooths a point in a sequence of 2D points by fitting a + * quadratic polynomial to a 5-point window (if possible) and evaluating it at + * the center. Essentially, it is a weighted local average where the reference + * point has the largest weight and neighbors have progressively smaller weights + * the farther they are from the reference. + * + * For points near the boundaries (indices 0, 1, n-2, n-1), where a full 5-point + * window is not available, a 3-point linear smoothing is applied instead. + * + * The 5-point quadratic filter uses the following coefficients: + * [-3, 12, 17, 12, -3] / 35 + * which preserves local peaks and slopes while reducing noise. + * + * @param pts A vector of 2D points (Point struct with x, y floats) to smooth. + * @param i The index of the point to smooth. + * @return The smoothed Point at index i. + * + * @note The function assumes pts has at least 3 points for linear smoothing and + * at least 5 points for quadratic smoothing. + */ + +Point getQuadraticTarget(const std::vector &pts, int i) { + int n = pts.size(); + + // 1. BOUNDARY FALLBACK: + // If we are too close to the end (indices 1 or n-2), we don't have 5 points. + // Fall back to standard Linear Laplacian (0.25, 0.5, 0.25). + if (i < 2 || i >= n - 2) { + Point prev = pts[i - 1]; + Point curr = pts[i]; + Point next = pts[i + 1]; + return 0.25f * prev + 0.5f * curr + 0.25f * next; + } + + // 2. QUADRATIC FILTER (Savitzky-Golay Window 5, Degree 2): + // Coefficients: [-3, 12, 17, 12, -3] / 35 + // This fits a local parabola and evaluates it at the center. + const Point &p2L = pts[i - 2]; // 2 Left + const Point &p1L = pts[i - 1]; // 1 Left + const Point &p = pts[i]; // Center + const Point &p1R = pts[i + 1]; // 1 Right + const Point &p2R = pts[i + 2]; // 2 Right + + return (-3.0f * p2L + 12.0f * p1L + 17.0f * p + 12.0f * p1R - 3.0f * p2R) / + 35.0f; +} + +void coupledSmooth(std::vector &contourA, std::vector &contourB) { + + std::vector> contours = {contourA, contourB}; + // 1. Build Spatial Grid for O(1) partner lookup + std::map> grid; + for (int c = 0; c < 2; ++c) { + for (int p = 0; p < (int)contours[c].size(); ++p) { + grid[{(int)std::round(contours[c][p].x), + (int)std::round(contours[c][p].y)}] + .push_back({c, p}); + } + } + + std::vector> targetPos = {contourA, contourB}; + + // Increased radius slightly to ensure we catch partners even on curves + float pairRadiusSq = 2.0 * 2.0; + + for (int c = 0; c < 2; ++c) { + // Skip endpoints (p=0 and p=last are usually locked anchors) + for (int p = 1; p < (int)contours[c].size() - 1; ++p) { + + Point myPos = contours[c][p]; + + // --- STEP A: Calculate My Ideal Position (Quadratic) --- + Point myTarget = getQuadraticTarget(contours[c], p); + + // --- STEP B: Find Partners & Calculate Their Ideal Positions --- + Point sumPartnerTargets = {0, 0}; + int partnerCount = 0; + + int gx = (int)std::round(myPos.x); + int gy = (int)std::round(myPos.y); + + // 3x3 Neighbor Search + for (int dy = -2; dy <= 2; ++dy) { + for (int dx = -2; dx <= 2; ++dx) { + auto it = grid.find({gx + dx, gy + dy}); + if (it == grid.end()) + continue; + + for (const auto &neighbor : it->second) { + if (neighbor.cIdx == c) + continue; // Ignore self + + Point otherPos = contours[neighbor.cIdx][neighbor.pIdx]; + + // If close enough to be a "Partner" + if (Point::distSq(myPos, otherPos) < pairRadiusSq) { + + // Check if partner is constrained + int op = neighbor.pIdx; + const auto &otherContour = contours[neighbor.cIdx]; + + Point oTarget; + if (op > 0 && op < (int)otherContour.size() - 1) { + + // Partner is free: Calculate THEIR Quadratic Target + oTarget = getQuadraticTarget(otherContour, op); + } else { + // Partner is locked: They want to stay put + oTarget = otherPos; + } + + sumPartnerTargets += oTarget; + partnerCount++; + } + } + } + } + + // --- STEP C: Consensus Averaging --- + // "I want to be a parabola" vs "My partner wants to be a parabola" + // We average the two perfect quadratic fits. + if (partnerCount > 0) { + targetPos[c][p] = (myTarget + sumPartnerTargets) / (1.0 + partnerCount); + } else { + targetPos[c][p] = myTarget; + } + } + } + + std::copy(targetPos[0].begin(), targetPos[0].end(), contourA.begin()); + std::copy(targetPos[1].begin(), targetPos[1].end(), contourB.begin()); +} + +void stitch_smooth(std::vector &vecA, std::vector &vecB) { + // 1. Map Vector B indices to Grid (Optimization) + // We map a coordinate to the INDEX in vector B + std::map mapB; + for (size_t i = 0; i < vecB.size(); ++i) { + mapB[{(int)std::round(vecB[i].x), (int)std::round(vecB[i].y)}] = i; + } + + // We store the calculated "Target Positions" here. + // We do NOT update in place immediately, or the math for the next point will + // be wrong. + struct Update { + int index; + Point newPos; + }; + std::vector updatesA; + std::vector updatesB; + + // --- Process Vector A (Snap to B's Geometry) --- + for (size_t i = 0; i < vecA.size(); ++i) { + int ax = (int)std::round(vecA[i].x); + int ay = (int)std::round(vecA[i].y); + + float minDst = std::numeric_limits::max(); + Point bestTarget = vecA[i]; + bool foundMatch = false; + + // Search 3x3 Neighborhood + for (int dy = -1; dy <= 1; ++dy) { + for (int dx = -1; dx <= 1; ++dx) { + auto it = mapB.find({ax + dx, ay + dy}); + if (it != mapB.end()) { + int bIdx = it->second; + + // CHECK FORWARD SEGMENT: B[bIdx] -> B[bIdx+1] + if (bIdx < (int)vecB.size() - 1) { + auto res = + getClosestPointOnSegment(vecA[i], vecB[bIdx], vecB[bIdx + 1]); + if (res.second < minDst) { + minDst = res.second; + bestTarget = res.first; + foundMatch = true; + } + } + + // CHECK BACKWARD SEGMENT: B[bIdx-1] -> B[bIdx] + if (bIdx > 0) { + auto res = + getClosestPointOnSegment(vecA[i], vecB[bIdx - 1], vecB[bIdx]); + if (res.second < minDst) { + minDst = res.second; + bestTarget = res.first; + foundMatch = true; + } + } + } + } + } + + if (foundMatch) { + // Move A to the midpoint between itself and the closest spot on B's line + Point mid = (vecA[i] + bestTarget) * 0.5f; + updatesA.push_back({(int)i, mid}); + } + } + + // --- Process Vector B (Snap to A's Geometry) --- + // (We need a map for A now to do the reverse) + std::map mapA; + for (size_t i = 0; i < vecA.size(); ++i) { + mapA[{(int)std::round(vecA[i].x), (int)std::round(vecA[i].y)}] = i; + } + + for (size_t i = 0; i < vecB.size(); ++i) { + int bx = (int)std::round(vecB[i].x); + int by = (int)std::round(vecB[i].y); + + float minDst = std::numeric_limits::max(); + Point bestTarget = vecB[i]; + bool foundMatch = false; + + for (int dy = -1; dy <= 1; ++dy) { + for (int dx = -1; dx <= 1; ++dx) { + auto it = mapA.find({bx + dx, by + dy}); + if (it != mapA.end()) { + int aIdx = it->second; + + // Check Forward Segment A + if (aIdx < (int)vecA.size() - 1) { + auto res = + getClosestPointOnSegment(vecB[i], vecA[aIdx], vecA[aIdx + 1]); + if (res.second < minDst) { + minDst = res.second; + bestTarget = res.first; + foundMatch = true; + } + } + // Check Backward Segment A + if (aIdx > 0) { + auto res = + getClosestPointOnSegment(vecB[i], vecA[aIdx - 1], vecA[aIdx]); + if (res.second < minDst) { + minDst = res.second; + bestTarget = res.first; + foundMatch = true; + } + } + } + } + } + + if (foundMatch) { + Point mid = (vecB[i] + bestTarget) * 0.5f; + updatesB.push_back({(int)i, mid}); + } + } + + if ((updatesA.size() == 0) | (updatesB.size() == 0)) { + return; + } + + // --- Apply Updates --- + for (const auto &u : updatesA) + vecA[u.index] = u.newPos; + for (const auto &u : updatesB) + vecB[u.index] = u.newPos; + + // --- OPTIONAL FINAL POLISH: Laplacian Smooth --- + // This removes any remaining high-frequency noise from the seam + // Only run this on the points that were actually touched/updated + // Formula: P_i = 0.25*P_prev + 0.5*P_curr + 0.25*P_next + + auto smoothVector = [](std::vector &pts, + const std::vector &updates) { + std::vector original = pts; + for (const auto &u : updates) { + int i = u.index; + if (i > 0 && i < (int)pts.size() - 1) { + pts[i] = 0.25f * original[i - 1] + 0.5f * original[i] + + 0.25f * original[i + 1]; + } + } + }; + + auto laplacianSmooth = [](std::vector &pts, + const std::vector &updates) { + std::vector original = pts; + for (const auto &u : updates) { + int i = u.index; + if (i > 0 && i < (int)pts.size() - 1) { + // Heavier weight on self (0.6) to preserve shape, but smooth noise (0.2 + // neighbors) + pts[i] = 0.2f * original[i - 1] + 0.6f * original[i] + + 0.2f * original[i + 1]; + } + } + }; + + smoothVector(vecA, updatesA); + smoothVector(vecB, updatesB); +} + +// Project p onto segment v-w. Returns closest point. +Point getClosestPoint(Point p, Point v, Point w) { + float l2 = Point::distSq(v, w); + if (l2 == 0.0f) + return v; + + float t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2; + // Allow slight extension (0.01) to help corners "kiss" + t = std::max(-0.1f, std::min(1.1f, t)); + + return {v.x + t * (w.x - v.x), v.y + t * (w.y - v.y)}; +} + +// --- Step 1: Boundary Locking --- +// Returns a mask: true = Locked (On Boundary), false = Free to move +std::vector> +createBoundaryMask(const std::vector> &contours, + Rect bounds) { + std::vector> locked(contours.size()); + float eps = 0.1; // Tolerance for "on the boundary" + + for (size_t c = 0; c < contours.size(); ++c) { + locked[c].resize(contours[c].size(), false); + for (size_t p = 0; p < contours[c].size(); ++p) { + Point pt = contours[c][p]; + // Check Left, Right, Top, Bottom + if (std::abs(pt.x - bounds.x) < eps || + std::abs(pt.x - (bounds.x + bounds.width - 1.0f)) < eps || + std::abs(pt.y - bounds.y) < eps || + std::abs(pt.y - (bounds.y + bounds.height - 1.0f)) < eps) { + + locked[c][p] = true; + } + } + } + return locked; +} + +// --- Corner Detection (Feature Preservation) --- +std::vector detectCorners(const std::vector &pts, + float angleThresholdDeg = 150.0) { + std::vector isCorner(pts.size(), false); + if (pts.size() < 3) + return isCorner; + float threshold = std::cos(angleThresholdDeg * M_PI / 180.0f); + + for (size_t i = 1; i < pts.size() - 1; ++i) { + Point v1 = normalize(pts[i] - pts[i - 1]); + Point v2 = normalize(pts[i + 1] - pts[i]); + if (v1.x * v2.x + v1.y * v2.y < threshold) + isCorner[i] = true; + } + // Endpoints are corners + isCorner[0] = true; + isCorner.back() = true; + return isCorner; +} + +// --- Selective Smoothing --- +void selectiveSmooth(std::vector &pts, + const std::vector &isLocked) { + std::vector original = pts; + for (size_t i = 1; i < pts.size() - 1; ++i) { + // DO NOT move if it's a Corner OR if it's Locked on the boundary + if (isLocked[i]) + continue; + + pts[i] = + 0.25f * original[i - 1] + 0.5f * original[i] + 0.25f * original[i + 1]; + } +} + +void coupledSmooth(std::vector> &contours, + const std::vector> &lockedMasks, + float pairRadiusSq = 2.25f) { + + SavitzkyGolay sg(2, 2); // radius, polynomial order + + // first fit + std::vector> smoothedContours; + for (int c = 0; c < (int)contours.size(); ++c) { + std::vector sc = sg.filter_wrap(contours[c]); + smoothedContours.push_back(sc); + } + + // 1. Build Spatial Grid to find partners quickly + std::map> grid; + for (int c = 0; c < (int)contours.size(); ++c) { + for (int p = 0; p < (int)contours[c].size(); ++p) { + grid[{static_cast(contours[c][p].x), + static_cast(contours[c][p].y)}] + .push_back({c, p}); + } + } + + // We calculate ALL targets before applying any updates to maintain stability + std::vector> targetPos = contours; + // float pairRadiusSq = 1.5f * 1.5f; // Radius to define "Connected/Paired" + + for (int c = 0; c < (int)contours.size(); ++c) { + for (int p = 1; p < (int)contours[c].size() - 1; ++p) { + // SKIP Constraints + if (lockedMasks[c][p]) + continue; + + Point myPos = contours[c][p]; + Point prev = contours[c][p - 1]; + Point next = contours[c][p + 1]; + + // 1. Calculate My Laplacian Target (Where I want to go to be smooth) + Point myTarget = smoothedContours[c][p]; + + // 2. Find Partners in OTHER contours + Point sumPartnerTargets = {0, 0}; + int partnerCount = 0; + + int gx{static_cast(myPos.x)}; + int gy{static_cast(myPos.y)}; + + // Check 3x3 grid + for (int dy = -1; dy <= 1; ++dy) { + for (int dx = -1; dx <= 1; ++dx) { + auto it = grid.find({gx + dx, gy + dy}); + if (it == grid.end()) + continue; + + for (const auto &neighbor : it->second) { + if (neighbor.cIdx == c) + continue; // Ignore self + + Point otherPos = contours[neighbor.cIdx][neighbor.pIdx]; + if (Point::distSq(myPos, otherPos) < pairRadiusSq) { + // Found a partner! + // Calculate where the PARTNER wants to go + // (We need safe access to partner's neighbors) + const auto &otherContour = contours[neighbor.cIdx]; + int op = neighbor.pIdx; + + // Only calculate partner target if they are not constrained + if (op > 0 && op < (int)otherContour.size() - 1 && + // !cornerMasks[neighbor.cIdx][op] && + !lockedMasks[neighbor.cIdx][op]) { + + Point oPrev = otherContour[op - 1]; + Point oNext = otherContour[op + 1]; + + Point oTarget = smoothedContours[neighbor.cIdx][op]; + + sumPartnerTargets += oTarget; + partnerCount++; + } else { + // If partner is constrained (e.g. corner), + // we should probably snap to THEM, not smooth them. + // For simplicity, we treat their current pos as their target. + sumPartnerTargets += otherPos; + partnerCount++; + } + } + } + } + } + + // 3. Average "My Desire" with "Partners' Desires" + if (partnerCount > 0) { + targetPos[c][p] = + (myTarget + sumPartnerTargets) / (1.0f + partnerCount); + } else { + // No partners, just smooth myself + targetPos[c][p] = myTarget; + } + } + } + + contours = targetPos; +} + +void coupled_smooth(std::vector> &contours, Rect bounds) { + auto lockedMasks = createBoundaryMask(contours, bounds); + coupledSmooth(contours, lockedMasks, 1.0f); +} + +// --- Main Solver --- +void pack_with_boundary_constraints(std::vector> &contours, + Rect bounds, int iterations) { + + // 1. Identify Locked Points (Boundary Constraint) + auto lockedMasks = createBoundaryMask(contours, bounds); + + // 2. Identify Feature Corners (Shape Constraint) + // std::vector> cornerMasks; + // for (const auto& c : contours) cornerMasks.push_back(detectCorners(c)); + + constexpr float radiusSq = 3.0f * 3.0f; // Search radius + + for (int iter = 0; iter < iterations; ++iter) { + + // Build Grid + std::map> grid; + for (int c = 0; c < (int)contours.size(); ++c) { + for (int p = 0; p < (int)contours[c].size(); ++p) { + grid[{static_cast(contours[c][p].x), + static_cast(contours[c][p].y)}] + .push_back({c, p}); + } + } + + std::vector> nextContours = contours; + + // Apply Forces + for (int c = 0; c < (int)contours.size(); ++c) { + for (int p = 0; p < (int)contours[c].size(); ++p) { + + // CRITICAL CHECK: If on boundary, skip all movement logic + if (lockedMasks[c][p]) + continue; + + Point currentPos = contours[c][p]; + Point sumTargets = {0, 0}; + int matchCount = 0; + + // Scan Neighborhood + int gx = static_cast(currentPos.x); + int gy = static_cast(currentPos.y); + + for (int dy = -1; dy <= 1; ++dy) { + for (int dx = -1; dx <= 1; ++dx) { + auto it = grid.find({gx + dx, gy + dy}); + if (it == grid.end()) + continue; + + for (const auto &neighbor : it->second) { + if (neighbor.cIdx == c) + continue; + const auto &other = contours[neighbor.cIdx]; + int idxB = neighbor.pIdx; + + auto checkSeg = [&](int s, int e) { + Point t = getClosestPoint(currentPos, other[s], other[e]); + if (Point::distSq(currentPos, t) < radiusSq) { + sumTargets += t; + matchCount++; + } + }; + if (idxB < (int)other.size() - 1) + checkSeg(idxB, idxB + 1); + if (idxB > 0) + checkSeg(idxB - 1, idxB); + } + } + } + + if (matchCount > 0) { + float stiffness{1.5f}; + nextContours[c][p] = + (sumTargets + currentPos * stiffness) / (matchCount + stiffness); + } + } + } + + contours = nextContours; + } +} + +} // namespace contours \ No newline at end of file diff --git a/src/wasm/modules/image/src/graph.cpp b/src/wasm/modules/image/src/graph.cpp index 6c6b6dd0a..79c4b5113 100644 --- a/src/wasm/modules/image/src/graph.cpp +++ b/src/wasm/modules/image/src/graph.cpp @@ -1,11 +1,14 @@ #include "graph.h" #include "Pixel.h" +#include "bezier.h" #include +#include #include - +#include +#include /* - *Graph class - manages Node class - */ + Graph class - manages Node class +*/ static inline float colorDistance(const ImageLib::RGBPixel &a, const ImageLib::RGBPixel &b) { @@ -22,11 +25,11 @@ static inline float colorDistance(const ImageLib::RGBPixel &a, } /* - *To quickly search m_nodes (std::vector) for the index of a node id - *create an std::unordered_map of node id - index pairs - *indexing time of std::vector by value is O(N) - *lookup time of std::unordered_map by key is O(log(N)) - */ +To quickly search m_nodes (std::vector) for the index of a node id +create an std::unordered_map of node id - index pairs +indexing time of std::vector by value is O(N) +lookup time of std::unordered_map by key is O(log(N)) +*/ void Graph::hash_node_ids() { for (int32_t i{0}; i < m_nodes->size(); i++) { const int32_t key{m_nodes->at(i)->id()}; @@ -103,8 +106,9 @@ void Graph::clear_unconnected_nodes() { void Graph::discover_edges(const std::vector ®ion_labels, const int32_t width, const int32_t height) { - // Moore 4-connected neighbourhood - constexpr int8_t dirs[4][2]{{1, 0}, {-1, 0}, {0, 1}, {0, -1}}; + // Moore 8-connected neighbourhood + constexpr int8_t dirs[8][2]{{1, 0}, {-1, 0}, {0, 1}, {0, -1}, + {1, 1}, {-1, -1}, {-1, 1}, {1, -1}}; int32_t rneigh[4]; @@ -113,7 +117,7 @@ void Graph::discover_edges(const std::vector ®ion_labels, const int32_t idx{y * width + x}; const int32_t rid{region_labels[idx]}; - for (int32_t k{0}; k < 4; ++k) { + for (int32_t k{0}; k < 8; ++k) { const int32_t nx{x + dirs[k][0]}; const int32_t ny{y + dirs[k][1]}; @@ -133,6 +137,181 @@ void Graph::discover_edges(const std::vector ®ion_labels, } } +void Graph::compute_contours() { + // overlap edge pixels + // then compute contours + + std::set> adjusted_neighbors{}; + + constexpr int8_t dirs[8][2]{{1, 0}, {-1, 0}, {0, 1}, {0, -1}, + {1, 1}, {-1, -1}, {-1, 1}, {1, -1}}; + + for (const Node_ptr &n : get_nodes()) { + if (n->area() == 0) + continue; + + std::vector> full_neighborhood; + std::vector> xywh; + + std::vector node_binary; + std::array xywh0 = n->create_binary_image(node_binary); + + full_neighborhood.push_back(node_binary); + xywh.push_back(xywh0); + + std::vector considered_neigbors; + for (const auto &neighbor : n->edges()) { + if (neighbor->area() == 0) + continue; + + // check if this neighbor pairing has already been addressed + std::pair id1 = std::make_pair(n->id(), neighbor->id()); + std::pair id2 = std::make_pair(neighbor->id(), n->id()); + auto _end{adjusted_neighbors.end()}; + auto _it1{adjusted_neighbors.find(id1)}; + auto _it2{adjusted_neighbors.find(id2)}; + + if (_it1 != _end || _it2 != _end) { + continue; + } + + considered_neigbors.push_back(neighbor); + + std::vector neighbor_binary; + std::array xywh1 = + neighbor->create_binary_image(neighbor_binary); + full_neighborhood.push_back(neighbor_binary); + xywh.push_back(xywh1); + + adjusted_neighbors.insert(id1); + } + + // address overlaps + std::vector neighborhood; + std::array bounds = {std::numeric_limits::max(), + std::numeric_limits::max(), -1, -1}; + for (auto &_xywh : xywh) { + if (_xywh[0] < bounds[0]) { + bounds[0] = _xywh[0]; + } // xmin + if (_xywh[1] < bounds[1]) { + bounds[1] = _xywh[1]; + } // ymin + if (_xywh[2] + _xywh[0] - 1 > bounds[2]) { + bounds[2] = _xywh[2] + _xywh[0] - 1; + } // xmax + if (_xywh[3] + _xywh[1] - 1 > bounds[3]) { + bounds[3] = _xywh[3] + _xywh[1] - 1; + } // ymax + } + + bounds[2] = bounds[2] - bounds[0] + 1; // w + bounds[3] = bounds[3] - bounds[1] + 1; // h + + // build joined neighborhood map + neighborhood.resize(bounds[2] * bounds[3], 0); + + for (int i = 0; i < full_neighborhood.size(); ++i) { + for (int y = 0; y < xywh[i][3]; ++y) { + for (int x = 0; x < xywh[i][2]; ++x) { + int global_y = y + xywh[i][1] - bounds[1]; + int global_x = x + xywh[i][0] - bounds[0]; + uint8_t val = full_neighborhood[i][y * xywh[i][2] + x]; + if (val != 0) { + neighborhood[global_y * bounds[2] + global_x] = (i + 1); + } + } + } + } + + // 0 = background, 1 = this node, 2+ = neighboring nodes + + // find touching edges + for (int y = 0; y < bounds[3]; ++y) { + for (int x = 0; x < bounds[2]; ++x) { + uint8_t val = neighborhood[y * bounds[2] + x]; + // check neighbors + if (val == 1) { + for (int32_t k{0}; k < 8; ++k) { + int32_t nx{x + dirs[k][0]}; + int32_t ny{y + dirs[k][1]}; + nx = std::clamp(nx, 0, bounds[2] - 1); + ny = std::clamp(ny, 0, bounds[3] - 1); + uint8_t n_val = neighborhood[ny * bounds[2] + nx]; + if ((n_val != val) & (n_val != 0)) { + // need a smarter approach to prevent pinching + bool is_too_thin = false; + // check around (nx,ny) if we can stretch into another region, + // then it's too thin + for (int32_t k{0}; k < 8; ++k) { + int32_t mx{nx + dirs[k][0]}; + int32_t my{ny + dirs[k][1]}; + mx = std::clamp(mx, 0, bounds[2] - 1); + my = std::clamp(my, 0, bounds[3] - 1); + uint8_t m_val = neighborhood[my * bounds[2] + mx]; + + if ((m_val != val) & (m_val != n_val)) { + is_too_thin = true; + } + } + + if (is_too_thin) { + considered_neigbors[n_val - 2]->add_edge_pixel( + XY{x + bounds[0], y + bounds[1]}); + } else { + n->add_edge_pixel(XY{nx + bounds[0], ny + bounds[1]}); + } + } + } + } + } + } + } + + // ask each Node to compute contours + for (const Node_ptr &n : get_nodes()) { + if (n->area() == 0) + continue; + n->compute_contour(); + } + + // smoothing + std::vector> all_contours; + for (const Node_ptr &n : get_nodes()) { + if (n->area() == 0) + continue; + + ColoredContours *c0 = &n->m_contours; + for (size_t i = 0; i < c0->contours.size(); ++i) { + all_contours.push_back(c0->contours[i]); + } + } + + contours::coupled_smooth(all_contours, + Rect{0.0f, 0.0f, static_cast(m_width), + static_cast(m_height)}); + + std::vector> all_curves; + fit_curve_reduction(all_contours, all_curves, 0.5f); + + int j = 0; + for (const Node_ptr &n : get_nodes()) { + if (n->area() == 0) + continue; + + ColoredContours *c0 = &n->m_contours; + for (size_t i = 0; i < c0->contours.size(); ++i) { + std::copy(all_contours[j].begin(), all_contours[j].end(), + c0->contours[i].begin()); + + c0->curves[i].resize(all_curves[j].size()); + std::copy(all_curves[j].begin(), all_curves[j].end(), + c0->curves[i].begin()); + j++; + } + } +} + void Graph::merge_small_area_nodes(const int32_t min_area) { int32_t counter{0}; while (!all_areas_bigger_than(min_area)) { @@ -144,19 +323,16 @@ void Graph::merge_small_area_nodes(const int32_t min_area) { std::back_inserter(neighbors)); ImageLib::RGBPixel col = n->color(); - // Sort by size -> a.area < b.area - // std::sort(neighbors.begin(), neighbors.end(), - // [](Node_ptr a, Node_ptr b) { return a->area() < b->area(); - // }); - // sort by size and color similarity - std::sort(neighbors.begin(), neighbors.end(), - [col](Node_ptr a, Node_ptr b) { - float cdista = colorDistance(a->color(), col); - float cdistb = colorDistance(b->color(), col); - return (static_cast(a->area()) + 10.f * cdista) < - (static_cast(b->area()) + 10.f * cdistb); - }); + std::sort( + neighbors.begin(), neighbors.end(), [col](Node_ptr a, Node_ptr b) { + float cdista = + ImageLib::RGBPixel::colorDistance(a->color(), col); + float cdistb = + ImageLib::RGBPixel::colorDistance(b->color(), col); + return (static_cast(a->area()) + 10.f * cdista) < + (static_cast(b->area()) + 10.f * cdistb); + }); int32_t idx{0}; // find first non-zero area neighbor diff --git a/src/wasm/modules/image/src/kmeans.cpp b/src/wasm/modules/image/src/kmeans.cpp index b724fe1a3..5a0bb6123 100644 --- a/src/wasm/modules/image/src/kmeans.cpp +++ b/src/wasm/modules/image/src/kmeans.cpp @@ -16,21 +16,6 @@ #include #include -static inline float colorDistance(const ImageLib::RGBAPixel &a, - const ImageLib::RGBAPixel &b) { - // sqrt un-necessary - return (a.red - b.red) * (a.red - b.red) + - (a.green - b.green) * (a.green - b.green) + - (a.blue - b.blue) * (a.blue - b.blue); -} - -static inline float colorDistance(const ImageLib::LABAPixel &a, - const ImageLib::LABAPixel &b) { - // sqrt un-necessary - return (a.l - b.l) * (a.l - b.l) + (a.a - b.a) * (a.a - b.a) + - (a.b - b.b) * (a.b - b.b); -} - static constexpr uint8_t COLOR_SPACE_OPTION_CIELAB{0}; static constexpr uint8_t COLOR_SPACE_OPTION_RGB{1}; @@ -44,7 +29,7 @@ void _process_dist_per_centroid(const ImageLib::Image &pixels, for (int j{start_centroid}; j < end_centroid; ++j) { std::transform(pixels.begin(), pixels.end(), _res.begin(), [¢roids, j](const PixelT &p) { - return colorDistance(p, centroids[j]); + return PixelT::colorDistance(p, centroids[j]); }); std::copy(_res.begin(), _res.end(), output[j].begin()); } @@ -105,7 +90,7 @@ void kMeansPlusPlusInit(const ImageLib::Image &pixels, // We don't need to recheck previous centroids; min_dist_sq already holds // the best distance to them. for (int j = 0; j < num_pixels; ++j) { - double d = colorDistance(pixels[j], centroids.back()); + double d = PixelT::colorDistance(pixels[j], centroids.back()); // If this new centroid is closer than the previous best, update the min // distance @@ -295,11 +280,13 @@ void kmeans(const uint8_t *data, uint8_t *out_data, int32_t *out_labels, // float dist{colorDistance(pixels[i], centroids[j])}; switch (color_space) { case COLOR_SPACE_OPTION_RGB: { - dist = colorDistance(pixels[i], centroids[j]); + dist = ImageLib::RGBAPixel::colorDistance(pixels[i], + centroids[j]); break; } case COLOR_SPACE_OPTION_CIELAB: { - dist = colorDistance(lab[i], centroids_lab[j]); + dist = ImageLib::LABAPixel::colorDistance(lab[i], + centroids_lab[j]); break; } } diff --git a/src/wasm/modules/image/src/kmeans_graph.cpp b/src/wasm/modules/image/src/kmeans_graph.cpp index b414f7d4e..f758999cf 100644 --- a/src/wasm/modules/image/src/kmeans_graph.cpp +++ b/src/wasm/modules/image/src/kmeans_graph.cpp @@ -2,9 +2,7 @@ // should actually be taken into account. #include "kmeans_graph.h" -#include "Image.h" -#include "PixelConverters.h" -#include "RGBAPixel.h" +#include "bezier.h" #include "contours.h" #include "graph.h" #include @@ -24,12 +22,6 @@ #include #include -struct ColoredContours : ContoursResult { - // inherits: contours, hierarchy, is_hole - - std::vector> colors; -}; - /* Flood fill */ int flood_fill(std::vector &label_array, std::vector ®ion_array, const uint8_t *color_array, @@ -113,7 +105,6 @@ void region_labeling(const uint8_t *data, std::vector &labels, void visualize_contours(const std::vector> &contours, ImageLib::Image> &results, int width, int height, int xmin = 0, int ymin = 0) { - auto index = [width](int x, int y) { return y * width + x; }; // Random generator for colors static std::mt19937 rng(std::random_device{}()); @@ -124,8 +115,8 @@ void visualize_contours(const std::vector> &contours, static_cast(dist(rng)), 255}; for (const auto &p : c) { - int32_t _x{p.x + xmin}; - int32_t _y{p.y + ymin}; + int32_t _x{static_cast(p.x) + xmin}; + int32_t _y{static_cast(p.y) + ymin}; // Ensure within bounds if (_x < 0 || _x >= width || _y < 0 || _y >= height) @@ -141,7 +132,7 @@ std::string contourToSVGPath(const std::vector &contour) { return ""; std::ostringstream path; - path << std::fixed << std::setprecision(0); + path << std::fixed << std::setprecision(2); // Move to the first point path << "M " << contour[0].x << " " << contour[0].y << " "; @@ -156,6 +147,27 @@ std::string contourToSVGPath(const std::vector &contour) { return path.str(); } +std::string contourToSVGCurve(const std::vector &curves) { + + if (curves.empty()) + return ""; + + std::ostringstream path; + path << std::fixed << std::setprecision(2); + + for (size_t i = 0; i < curves.size(); ++i) { + const auto &c = curves[i]; + if (i == 0) + path << "M " << c.p0.x << " " << c.p0.y << " "; + path << "Q " << c.p1.x << " " << c.p1.y << " " << c.p2.x << " " << c.p2.y + << " "; + } + + // Close the path + path << "Z"; + return path.str(); +} + std::string contoursResultToSVG(const ColoredContours &result, const int width, const int height) { std::ostringstream svg; @@ -163,8 +175,8 @@ std::string contoursResultToSVG(const ColoredContours &result, const int width, "width=\"" << width << "\" height=\"" << height << "\">\n"; - for (size_t i = 0; i < result.contours.size(); ++i) { - std::string pathData = contourToSVGPath(result.contours[i]); + for (size_t i = 0; i < result.curves.size(); ++i) { + std::string pathData = contourToSVGCurve(result.curves[i]); const auto &px = result.colors[i]; std::ostringstream oss; @@ -205,7 +217,7 @@ char *kmeans_clustering_graph(uint8_t *data, int32_t *labels, const int width, // 2. initialize Graph from all Nodes std::unique_ptr> node_ptr = std::make_unique>(std::move(nodes)); - Graph G(node_ptr); + Graph G(node_ptr, width, height); // 3. Discover node adjacencies - add edges to Graph G.discover_edges(region_labels, width, height); @@ -226,46 +238,29 @@ char *kmeans_clustering_graph(uint8_t *data, int32_t *labels, const int width, } // 6. Contours - ColoredContours all_contours; + // graph will manage computing contours + G.compute_contours(); + // accumulate all contours for svg export + ColoredContours all_contours; for (auto &n : G.get_nodes()) { if (n->area() == 0) continue; - - std::array xywh; - std::vector binary; - - xywh = n->create_binary_image(binary); - - int xmin = xywh[0]; - int ymin = xywh[1]; - int bw = xywh[2]; - int bh = xywh[3]; - - ContoursResult contour_res = contours::find_contours(binary, bw, bh); - - if (draw_contour_borders) { - visualize_contours(contour_res.contours, results, width, height, xmin, - ymin); - } else { - // shift contour coordinates to image space - for (size_t cidx = 0; cidx < contour_res.contours.size(); ++cidx) { - auto &contour = contour_res.contours[cidx]; - for (auto &p : contour) { - p.x += xmin; - p.y += ymin; - } - - // pick the color from the first pixel of the contour in the recolored - // image - const auto &first_px = contour[0]; - ImageLib::RGBAPixel col = results(first_px.x, first_px.y); - - all_contours.contours.push_back(contour); - all_contours.hierarchy.push_back(contour_res.hierarchy[cidx]); - all_contours.is_hole.push_back(contour_res.is_hole[cidx]); - all_contours.colors.push_back(col); - } + ColoredContours node_contours = n->get_contours(); + for (auto &c : node_contours.contours) { + all_contours.contours.push_back(c); + } + for (auto &c : node_contours.hierarchy) { + all_contours.hierarchy.push_back(c); + } + for (bool b : node_contours.is_hole) { + all_contours.is_hole.push_back(b); + } + for (auto &c : node_contours.colors) { + all_contours.colors.push_back(c); + } + for (auto &c : node_contours.curves) { + all_contours.curves.push_back(c); } } diff --git a/src/wasm/modules/image/src/node.cpp b/src/wasm/modules/image/src/node.cpp index 608aa69e1..9266b5ad1 100644 --- a/src/wasm/modules/image/src/node.cpp +++ b/src/wasm/modules/image/src/node.cpp @@ -75,6 +75,21 @@ std::array Node::bounding_box_xywh() const { } } + for (auto &p : m_edge_pixels) { + if (p.x < x_min) { + x_min = p.x; + } + if (p.x > x_max) { + x_max = p.x; + } + if (p.y < y_min) { + y_min = p.y; + } + if (p.y > y_max) { + y_max = p.y; + } + } + const int32_t w{x_max - x_min + 1}; const int32_t h{y_max - y_min + 1}; @@ -93,16 +108,74 @@ Node::create_binary_image(std::vector &binary) const { binary[_y * xywh[2] + _x] = 1; } + // include the edge pixels to ensure contour overlap with neighbor + for (auto &p : m_edge_pixels) { + int32_t _x = p.x - xywh[0]; + int32_t _y = p.y - xywh[1]; + binary[_y * xywh[2] + _x] = 1; + } + return xywh; } +void Node::clear_contour(void) { + m_contours.contours.clear(); + m_contours.hierarchy.clear(); + m_contours.is_hole.clear(); + m_contours.colors.clear(); + m_contours.curves.clear(); +} + +void Node::compute_contour(void) { + // return list of all contours present in Node. + // usually 1 sometimes more if holes are present + + clear_contour(); + + std::vector binary; + std::array xywh{create_binary_image(binary)}; + + int xmin = xywh[0]; + int ymin = xywh[1]; + int bw = xywh[2]; + int bh = xywh[3]; + + ContoursResult contour_res = contours::find_contours(binary, bw, bh); + + for (size_t cidx = 0; cidx < contour_res.contours.size(); ++cidx) { + auto &contour = contour_res.contours[cidx]; + for (auto &p : contour) { + p.x += xmin; + p.y += ymin; + } + + // if (contour_res.is_hole[cidx]) { continue; } + + ImageLib::RGBPixel _col = color(); + ImageLib::RGBAPixel col{_col.red, _col.green, _col.blue, 255}; + m_contours.contours.push_back(contour); + m_contours.hierarchy.push_back(contour_res.hierarchy[cidx]); + m_contours.is_hole.push_back(contour_res.is_hole[cidx]); + m_contours.colors.push_back(col); + } + + m_contours.curves.resize(m_contours.contours.size()); +} + void Node::add_pixels(const std::vector &new_pixels) { for (auto &c : new_pixels) { m_pixels->push_back(c); } } +void Node::add_edge_pixel(const XY edge_pixel) { + m_edge_pixels.insert(edge_pixel); +} + +void Node::clear_edge_pixels() { m_edge_pixels.clear(); } + void Node::clear_all() { m_edges.clear(); m_pixels->clear(); + m_edge_pixels.clear(); }