flatcv (empty) → 0.3.0

raw patch · 10 files changed

+7279/−0 lines, 10 filesdep +basedep +flatcvdep +hspecsetup-changed

Dependencies added: base, flatcv, hspec, vector

Files

• LICENSE +15/−0
• README.md +120/−0
• Setup.hs +4/−0
• changelog.md +23/−0
• flatcv.cabal +63/−0
• flatcv.h +546/−0
• flatcv_impl.c +5133/−0
• src/FlatCV.hs +610/−0
• src/FlatCV/Raw.chs +515/−0
• test/Spec.hs +250/−0

@@ -0,0 +1,15 @@+ISC License++Copyright (c) 2025 Adrian Sieber++Permission to use, copy, modify, and/or distribute this software for any+purpose with or without fee is hereby granted, provided that the above+copyright notice and this permission notice appear in all copies.++THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH+REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY+AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,+INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM+LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR+OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR+PERFORMANCE OF THIS SOFTWARE.

@@ -0,0 +1,120 @@+# FlatCV Haskell Bindings++Haskell bindings for the [FlatCV](https://github.com/ad-si/FlatCV)+image processing library.+++## Usage++```hs+import FlatCV+import qualified Data.Vector.Storable as V++main :: IO ()+main = do+  -- Assuming you have RGBA image data as a Vector Word8+  let+    width = 640+    height = 480+    -- imageData :: V.Vector Word8++  -- Convert to grayscale+  grayData <- grayscale width height imageData++  -- Apply Gaussian blur+  blurredData <- gaussianBlur width height 2.0 imageData++  -- Detect edges+  edges <- sobelEdgeDetection width height 4 imageData++  -- Extract a document from a photo+  (outW, outH, docData) <- extractDocumentAuto width height imageData+```+++## Available Functions++### Color Conversion++- `grayscale` - Convert RGBA to grayscale (4-channel output)+- `grayscaleStretch` - Grayscale with contrast stretching+- `rgbaToGrayscale` - Convert to single-channel grayscale+- `singleToMultichannel` - Convert grayscale to RGBA+++### Filtering++- `gaussianBlur` - Apply Gaussian blur+++### Thresholding++- `otsuThreshold` - Otsu's automatic thresholding+- `bwSmart` - Smart black & white conversion+- `applyGlobalThreshold` - Apply fixed threshold+++### Edge Detection++- `sobelEdgeDetection` - Sobel edge detection+- `foerstnerCorner` - Foerstner corner detection+++### Morphological Operations++- `binaryDilationDisk` - Binary dilation+- `binaryErosionDisk` - Binary erosion+- `binaryClosingDisk` - Binary closing+- `binaryOpeningDisk` - Binary opening+++### Geometric Transforms++- `flipX`, `flipY` - Flip horizontally/vertically+- `transpose`, `transverse` - Transpose operations+- `rotate90CW`, `rotate180`, `rotate270CW` - Rotations+- `resize` - Resize with bilinear interpolation+++### Cropping & Trimming++- `crop` - Crop rectangular region+- `trim` - Auto-trim borders+- `trimThreshold` - Trim with threshold tolerance+++### Perspective Transform++- `calculatePerspectiveTransform` - Calculate transform matrix+- `applyMatrix3x3` - Apply 3x3 transformation+- `detectCorners` - Detect document corners+++### Document Extraction++- `extractDocument` - Extract with fixed output size+- `extractDocumentAuto` - Extract with automatic sizing+++### Drawing++- `drawCircle` - Draw circle outline+- `drawDisk` - Draw filled circle+- `addBorder` - Add border around image+++### Analysis++- `generateHistogram` - Generate histogram visualization+++### Segmentation++- `watershedSegmentation` - Watershed segmentation+++### Utility++- `getExifOrientation` - Read JPEG EXIF orientation+- `parseHexColor` - Parse hex color string+- `convertToBinary` - Convert to binary with custom colors

@@ -0,0 +1,4 @@+import Distribution.Simple++main :: IO ()+main = defaultMain

@@ -0,0 +1,23 @@+# Changelog++All notable changes to this project will be documented in this file.+++## 2025-01-15 - 0.1++- Initial release of Haskell bindings for FlatCV+- Bindings for all public FlatCV functions:+  - Color conversion (grayscale, contrast stretching)+  - Filtering (Gaussian blur)+  - Thresholding (Otsu, global threshold)+  - Edge detection (Sobel, Foerstner)+  - Morphological operations (dilation, erosion, opening, closing)+  - Geometric transforms (flip, rotate, transpose, resize)+  - Cropping and trimming+  - Perspective transforms and document extraction+  - Drawing primitives (circles, disks, borders)+  - Histogram generation+  - Watershed segmentation+  - EXIF orientation reading+- High-level Haskell API with `Data.Vector.Storable` integration+- Type-safe wrappers for all C structures (Corners, Matrix3x3, Point2D)

@@ -0,0 +1,63 @@+cabal-version:      2.4+name:               flatcv+version:            0.3.0+synopsis:           Haskell bindings for FlatCV image processing library+description:+    FlatCV is a lightweight, dependency-free image processing library+    written in C. This package provides Haskell bindings to the FlatCV+    library, enabling fast image operations like grayscale conversion,+    Gaussian blur, edge detection, perspective transforms, and more.+    .+    The library is particularly useful for document scanning and+    image preprocessing tasks.++license:            ISC+license-file:       LICENSE+author:             Adrian Sieber+maintainer:         mail@adriansieber.com+copyright:          2025 Adrian Sieber+category:           Graphics, Image+build-type:         Simple+extra-source-files:+    flatcv.h+    flatcv_impl.c+    changelog.md+    README.md++homepage:           https://github.com/ad-si/FlatCV+bug-reports:        https://github.com/ad-si/FlatCV/issues++source-repository head+    type:     git+    location: https://github.com/ad-si/FlatCV.git++library+    exposed-modules:    FlatCV+    other-modules:      FlatCV.Raw+    build-tool-depends: c2hs:c2hs >= 0.28+    hs-source-dirs:     src+    default-language:   Haskell2010+    default-extensions:+        RecordWildCards++    build-depends:+        base >= 4.12 && < 5,+        vector >= 0.12++    include-dirs:       .+    c-sources:          flatcv_impl.c++    cc-options:         -O2+    ghc-options:        -Wall++test-suite flatcv-test+    type:               exitcode-stdio-1.0+    main-is:            Spec.hs+    hs-source-dirs:     test+    default-language:   Haskell2010+    ghc-options:        -Wall -threaded -rtsopts -with-rtsopts=-N+    build-depends:+        base >= 4.12 && < 5,+        flatcv,+        hspec >= 2.7,+        vector >= 0.12

@@ -0,0 +1,546 @@+/* FlatCV - Amalgamated public header (auto-generated) */+/*+ * ISC License+ * + * Copyright (c) 2025 Adrian Sieber+ * + * Permission to use, copy, modify, and/or distribute this software for any+ * purpose with or without fee is hereby granted, provided that the above+ * copyright notice and this permission notice appear in all copies.+ * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH+ * REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY+ * AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,+ * INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM+ * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR+ * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR+ * PERFORMANCE OF THIS SOFTWARE.+ */++#ifndef FLATCV_H+#define FLATCV_H+#define FLATCV_AMALGAMATION+// File: include/1_types.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++typedef struct {+  double x;+  double y;+} Point2D;++typedef struct {+  double tl_x, tl_y;+  double tr_x, tr_y;+  double br_x, br_y;+  double bl_x, bl_y;+} Corners;++typedef struct {+  double m00, m01, m02;+  double m10, m11, m12;+  double m20, m21, m22;+} Matrix3x3;+// File: include/binary_closing_disk.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_binary_dilation_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+);++uint8_t *fcv_binary_erosion_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+);++uint8_t *fcv_binary_closing_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+);++uint8_t *fcv_binary_opening_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+);+// File: include/conversion.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdbool.h>+#include <stdint.h>++uint8_t *fcv_apply_gaussian_blur(+  uint32_t width,+  uint32_t height,+  double radius,+  uint8_t const * const data+);++uint8_t *fcv_grayscale(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++uint8_t *fcv_grayscale_stretch(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++void fcv_apply_global_threshold(+  uint32_t img_length,+  uint8_t *data,+  uint8_t threshold+);++uint8_t *fcv_otsu_threshold_rgba(+  uint32_t width,+  uint32_t height,+  bool use_double_threshold,+  uint8_t const * const data+);++uint8_t *fcv_bw_smart(+  uint32_t width,+  uint32_t height,+  bool use_double_threshold,+  uint8_t const * const data+);++uint8_t *fcv_resize(+  uint32_t width,+  uint32_t height,+  double scale_x,+  double scale_y,+  uint32_t* out_width,+  uint32_t* out_height,+  uint8_t const * const data+);+// File: include/convert_to_binary.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_convert_to_binary(+  const uint8_t *image_data,+  int32_t width,+  int32_t height,+  const char *foreground_hex,+  const char *background_hex+);+// File: include/corner_peaks.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#include "1_types.h"+#endif++typedef struct {+  Point2D *points;+  uint32_t count;+} CornerPeaks;++CornerPeaks *fcv_corner_peaks(+  uint32_t width,+  uint32_t height,+  uint8_t const *data,+  uint32_t min_distance,+  double accuracy_thresh,+  double roundness_thresh+);+// File: include/crop.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_crop(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  uint8_t const * const data,+  uint32_t x,+  uint32_t y,+  uint32_t new_width,+  uint32_t new_height+);+// File: include/draw.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++void fcv_set_circle_pixel(+  uint8_t *data,+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  int32_t px,+  int32_t py,+  uint8_t r,+  uint8_t g,+  uint8_t b+);++void fcv_draw_circle(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  const char *hex_color,+  double radius,+  double center_x,+  double center_y,+  uint8_t *data+);++void fcv_draw_disk(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  const char *hex_color,+  double radius,+  double center_x,+  double center_y,+  uint8_t *data+);++void fcv_draw_circle_points(+  uint8_t *data,+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  int32_t cx,+  int32_t cy,+  int32_t x,+  int32_t y,+  uint8_t r,+  uint8_t g,+  uint8_t b+);++void fcv_fill_disk_lines(+  uint8_t *data,+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  int32_t cx,+  int32_t cy,+  int32_t x,+  int32_t y,+  uint8_t r,+  uint8_t g,+  uint8_t b+);++uint8_t *fcv_add_border(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  const char *hex_color,+  uint32_t border_width,+  uint8_t *input_data,+  uint32_t *output_width,+  uint32_t *output_height+);+// File: include/exif.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++/**+ * Get EXIF orientation from a JPEG file.+ *+ * @param filename Path to the JPEG file.+ * @return Orientation value (1-8), or 1 if not found or error.+ */+int32_t fcv_get_exif_orientation(char const *filename);+// File: include/extract_document.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_extract_document(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data,+  uint32_t output_width,+  uint32_t output_height+);++uint8_t *fcv_extract_document_auto(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data,+  uint32_t *output_width,+  uint32_t *output_height+);+// File: include/flip.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_flip_x(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++uint8_t *fcv_flip_y(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++uint8_t *fcv_transpose(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++uint8_t *fcv_transverse(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);+// File: include/foerstner_corner.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_foerstner_corner(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data,+  double sigma+);+// File: include/histogram.h+#ifndef HISTOGRAM_H+#define HISTOGRAM_H++#include <stdint.h>++/**+ * Generate a histogram visualization image from input image data.+ * For grayscale images, creates a single histogram.+ * For RGB(A) images, creates overlapping histograms for each channel.+ *+ * @param width Width of the input image.+ * @param height Height of the input image.+ * @param channels Number of channels in the input image (1, 3, or 4).+ * @param data Pointer to the input pixel data.+ * @param out_width Pointer to store the output histogram width.+ * @param out_height Pointer to store the output histogram height.+ * @return Pointer to the histogram image data (RGBA format).+ */+uint8_t *fcv_generate_histogram(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  uint8_t const *const data,+  uint32_t *out_width,+  uint32_t *out_height+);++#endif // HISTOGRAM_H+// File: include/parse_hex_color.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++void fcv_parse_hex_color(+  const char *hex_color,+  uint8_t *r,+  uint8_t *g,+  uint8_t *b+);+// File: include/perspectivetransform.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#include "1_types.h"+#endif++Matrix3x3 *fcv_calculate_perspective_transform(+  Corners *src_corners,+  Corners *dst_corners+);++uint8_t *fcv_apply_matrix_3x3(+  int32_t in_width,+  int32_t in_height,+  uint8_t *in_data,+  int32_t out_width,+  int32_t out_height,+  Matrix3x3 *tmat+);++Corners fcv_detect_corners(const uint8_t *image, int32_t width, int32_t height);+Corners* fcv_detect_corners_ptr(const uint8_t *image, int32_t width, int32_t height);+// File: include/rgba_to_grayscale.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++// Avoid using floating point arithmetic by pre-multiplying the weights+#define R_WEIGHT 76  // 0.299 * 256+#define G_WEIGHT 150 // 0.587 * 256+#define B_WEIGHT 30  // 0.114 * 256++uint8_t *fcv_rgba_to_grayscale(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);+// File: include/rotate.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_rotate_90_cw(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++uint8_t *fcv_rotate_180(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);++uint8_t *fcv_rotate_270_cw(+  uint32_t width,+  uint32_t height,+  uint8_t const * const data+);+// File: include/single_to_multichannel.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_single_to_multichannel(+  uint32_t width,+  uint32_t height,+  uint8_t const *const data+);+// File: include/sobel_edge_detection.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++uint8_t *fcv_sobel_edge_detection(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  uint8_t const * const data+);+// File: include/sort_corners.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#include "1_types.h"+#endif++Corners sort_corners(+  uint32_t width, uint32_t height,+  uint32_t out_width, uint32_t out_height,+  Point2D *corners, uint32_t num_corners,+  Point2D *result+);+// File: include/trim.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#endif++#include <stdint.h>++/**+ * Trim border pixels that have the same color.+ *+ * @param width Pointer to width of the image (will be updated).+ * @param height Pointer to height of the image (will be updated).+ * @param channels Number of channels in the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the new trimmed image data.+ */+uint8_t *fcv_trim(+  int32_t *width,+  int32_t *height,+  uint32_t channels,+  uint8_t const *const data+);++/**+ * Trim border pixels that are within a threshold of the reference color.+ * This is useful for images with JPEG artifacts or slight vignetting.+ *+ * @param width Pointer to width of the image (will be updated).+ * @param height Pointer to height of the image (will be updated).+ * @param channels Number of channels in the image.+ * @param data Pointer to the pixel data.+ * @param threshold_percent Tolerance percentage (0-100). A value of 2 means+ *        pixels within 2% (≈5 units) of the reference color will be trimmed.+ * @return Pointer to the new trimmed image data.+ */+uint8_t *fcv_trim_threshold(+  int32_t *width,+  int32_t *height,+  uint32_t channels,+  uint8_t const *const data,+  double threshold_percent+);+// File: include/watershed_segmentation.h+#ifndef FLATCV_AMALGAMATION+#pragma once+#include "1_types.h"+#endif++#include <stdint.h>+#include <stdbool.h>++uint8_t *fcv_watershed_segmentation(+  uint32_t width,+  uint32_t height,+  uint8_t const * const grayscale_data,+  Point2D* markers,+  uint32_t num_markers,+  bool create_boundaries+);+#endif /* FLATCV_H */

@@ -0,0 +1,5133 @@+/* FlatCV - Amalgamated implementation (auto-generated) */+/*+ * ISC License+ * + * Copyright (c) 2025 Adrian Sieber+ * + * Permission to use, copy, modify, and/or distribute this software for any+ * purpose with or without fee is hereby granted, provided that the above+ * copyright notice and this permission notice appear in all copies.+ * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH+ * REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY+ * AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,+ * INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM+ * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR+ * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR+ * PERFORMANCE OF THIS SOFTWARE.+ */++#define FLATCV_AMALGAMATION+#include "flatcv.h"+// File: src/binary_closing_disk.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "conversion.h"+#include "perspectivetransform.h"+#else+#include "flatcv.h"+#endif++uint8_t *fcv_binary_dilation_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+) {+  if (!image_data || width <= 0 || height <= 0 || radius < 0) {+    return NULL;+  }++  // Check for overflow: width * height (width and height already validated > 0)+  if ((size_t)width > SIZE_MAX / (size_t)height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * (size_t)height;++  uint8_t *result = malloc(num_pixels);+  if (!result) {+    return NULL;+  }++  // Initialize result to 0 (black)+  memset(result, 0, num_pixels);++  // For each pixel in the image+  for (int32_t y = 0; y < height; y++) {+    for (int32_t x = 0; x < width; x++) {+      int32_t idx = y * width + x;++      // If current pixel is white (255), dilate it+      if (image_data[idx] == 255) {+        // Apply disk-shaped structuring element+        for (int32_t dy = -radius; dy <= radius; dy++) {+          for (int32_t dx = -radius; dx <= radius; dx++) {+            // Check if point32_t is within disk (use slightly larger radius for+            // better connectivity)+            double r_eff = radius + 0.5;+            if (dx * dx + dy * dy <= r_eff * r_eff) {+              int32_t ny = y + dy;+              int32_t nx = x + dx;++              // Check bounds+              if (ny >= 0 && ny < height && nx >= 0 && nx < width) {+                int32_t nidx = ny * width + nx;+                result[nidx] = 255;+              }+            }+          }+        }+      }+    }+  }++  return result;+}++static uint8_t *binary_erosion_disk_internal(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius,+  bool replicate_border+) {+  if (!image_data || width <= 0 || height <= 0 || radius < 0) {+    return NULL;+  }++  // Check for overflow: width * height (width and height already validated > 0)+  if ((size_t)width > SIZE_MAX / (size_t)height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * (size_t)height;++  uint8_t *result = malloc(num_pixels);+  if (!result) {+    return NULL;+  }++  // Initialize result to 0 (black)+  memset(result, 0, num_pixels);++  // For each pixel in the image+  for (int32_t y = 0; y < height; y++) {+    for (int32_t x = 0; x < width; x++) {+      int32_t idx = y * width + x;++      // Check if structuring element fits entirely within white pixels+      bool can_erode = true;++      for (int32_t dy = -radius; dy <= radius && can_erode; dy++) {+        for (int32_t dx = -radius; dx <= radius && can_erode; dx++) {+          // Check if point32_t is within disk (use slightly larger radius for+          // better connectivity)+          double r_eff = radius + 0.5;+          if (dx * dx + dy * dy <= r_eff * r_eff) {+            int32_t ny = y + dy;+            int32_t nx = x + dx;++            uint8_t neighbor_value;++            // Check bounds+            if (ny < 0 || ny >= height || nx < 0 || nx >= width) {+              if (replicate_border) {+                // Replicate border: clamp coordinates to valid range+                int32_t clamped_y =+                  ny < 0 ? 0 : (ny >= height ? height - 1 : ny);+                int32_t clamped_x = nx < 0 ? 0 : (nx >= width ? width - 1 : nx);+                neighbor_value = image_data[clamped_y * width + clamped_x];+              }+              else {+                // Default: treat out-of-bounds as black+                neighbor_value = 0;+              }+            }+            else {+              neighbor_value = image_data[ny * width + nx];+            }++            if (neighbor_value != 255) {+              can_erode = false;+            }+          }+        }+      }++      if (can_erode) {+        result[idx] = 255;+      }+    }+  }++  return result;+}++uint8_t *fcv_binary_erosion_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+) {+  // Public API uses default behavior: treat out-of-bounds as black+  return binary_erosion_disk_internal(image_data, width, height, radius, false);+}++uint8_t *fcv_binary_closing_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+) {+  if (!image_data || width <= 0 || height <= 0 || radius < 0) {+    return NULL;+  }++  // Step 1: Dilation+  uint8_t *dilated =+    fcv_binary_dilation_disk(image_data, width, height, radius);+  if (!dilated) {+    return NULL;+  }++  // Step 2: Erosion of the dilated image+  // Use replicate_border=true to preserve original white pixels at borders+  // In a closing operation, out-of-bounds pixels use "replicate" border mode+  // which clamps coordinates to valid range, preventing border pixels from+  // being erroneously eroded due to artificial black boundary+  uint8_t *result =+    binary_erosion_disk_internal(dilated, width, height, radius, true);++  // Free intermediate result+  free(dilated);++  return result;+}++uint8_t *fcv_binary_opening_disk(+  uint8_t const *image_data,+  int32_t width,+  int32_t height,+  int32_t radius+) {+  if (!image_data || width <= 0 || height <= 0 || radius < 0) {+    return NULL;+  }++  // Step 1: Erosion+  uint8_t *eroded = fcv_binary_erosion_disk(image_data, width, height, radius);+  if (!eroded) {+    return NULL;+  }++  // Step 2: Dilation of the eroded image+  uint8_t *result = fcv_binary_dilation_disk(eroded, width, height, radius);++  // Free intermediate result+  free(eroded);++  return result;+}+// File: src/conversion.c+#include <assert.h>+#include <math.h>+#ifndef M_PI+#define M_PI 3.14159265358979323846+#endif+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "conversion.h"+#include "draw.h"+#include "parse_hex_color.h"+#include "perspectivetransform.h"+#include "rgba_to_grayscale.h"+#include "single_to_multichannel.h"+#include "sobel_edge_detection.h"+#else+#include "flatcv.h"+#endif++/**+ * Convert raw RGBA row-major top-to-bottom image data+ * to RGBA row-major top-to-bottom grayscale image data.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the grayscale image data.+ */+uint8_t *+fcv_grayscale(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *grayscale_data = malloc(img_length_byte);++  if (!grayscale_data) { // Memory allocation failed+    return NULL;+  }++  // Process each pixel row by row+  for (uint32_t i = 0; i < width * height; i++) {+    uint32_t rgba_index = i * 4;++    uint8_t r = data[rgba_index];+    uint8_t g = data[rgba_index + 1];+    uint8_t b = data[rgba_index + 2];++    uint8_t gray = (r * R_WEIGHT + g * G_WEIGHT + b * B_WEIGHT) >> 8;++    grayscale_data[rgba_index] = gray;+    grayscale_data[rgba_index + 1] = gray;+    grayscale_data[rgba_index + 2] = gray;+    grayscale_data[rgba_index + 3] = 255;+  }++  return grayscale_data;+}++/**+ * Convert raw RGBA row-major top-to-bottom image data+ * to RGBA row-major top-to-bottom grayscale image data+ * with a stretched contrast range.+ * Set the 1.5625 % darkest pixels to 0 and the 1.5625 % brightest to 255.+ * Uses this specific value for speed: x * 1.5625 % = x >> 6+ * The rest of the pixel values are linearly scaled to the range [0, 255].++ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the grayscale image data.+ */+uint8_t *fcv_grayscale_stretch(+  uint32_t width,+  uint32_t height,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;+  if (img_length_px > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = img_length_px * 4;++  uint8_t *grayscale_data = malloc(img_length_byte);++  if (!grayscale_data) { // Memory allocation failed+    return NULL;+  }+  // Ignore 1.5625 % of the pixels+  uint32_t num_pixels_to_ignore = img_length_px >> 6;++  uint8_t *gray_values = malloc(img_length_px);+  if (!gray_values) { // Memory allocation failed+    free(grayscale_data);+    return NULL;+  }++  // Process each pixel row by row to get grayscale values+  for (uint32_t i = 0; i < img_length_px; i++) {+    uint32_t rgba_index = i * 4;++    uint8_t r = data[rgba_index];+    uint8_t g = data[rgba_index + 1];+    uint8_t b = data[rgba_index + 2];++    gray_values[i] = (r * R_WEIGHT + g * G_WEIGHT + b * B_WEIGHT) >> 8;+  }++  // Use counting sort to find the 1.5625% darkest and brightest pixels+  uint32_t histogram[256] = {0};+  for (uint32_t i = 0; i < img_length_px; i++) {+    histogram[gray_values[i]]++;+  }++  uint32_t cumulative_count = 0;+  uint8_t min_val = 0;+  for (uint32_t i = 0; i < 256; i++) {+    cumulative_count += histogram[i];+    if (cumulative_count > num_pixels_to_ignore) {+      min_val = i;+      break;+    }+  }++  cumulative_count = 0;+  uint8_t max_val = 255;+  for (int32_t i = 255; i >= 0; i--) {+    cumulative_count += histogram[i];+    if (cumulative_count > num_pixels_to_ignore) {+      max_val = i;+      break;+    }+  }++  free(gray_values);++  uint8_t range = max_val - min_val;++  // Process each pixel row by row+  for (size_t i = 0; i < img_length_px; i++) {+    uint32_t rgba_index = i * 4;++    uint8_t r = data[rgba_index];+    uint8_t g = data[rgba_index + 1];+    uint8_t b = data[rgba_index + 2];++    uint8_t gray = (r * R_WEIGHT + g * G_WEIGHT + b * B_WEIGHT) >> 8;++    if (range == 0) {+      // All pixels have same value, preserve as-is+      // (no stretching possible)+    }+    else if (gray < min_val) {+      gray = 0;+    }+    else if (gray > max_val) {+      gray = 255;+    }+    else {+      gray = (gray - min_val) * 255 / range;+    }++    grayscale_data[rgba_index] = gray;+    grayscale_data[rgba_index + 1] = gray;+    grayscale_data[rgba_index + 2] = gray;+    grayscale_data[rgba_index + 3] = 255;+  }++  return grayscale_data;+}++/**+ * Apply a global threshold to the image data.+ *+ * @param img_length_px Length of the image data in pixels.+ * @param data Pointer to the image data.+ * @param threshold Threshold value.+ *+ */+void fcv_apply_global_threshold(+  uint32_t img_length_px,+  uint8_t *data,+  uint8_t threshold+) {+  for (uint32_t i = 0; i < img_length_px; i++) {+    data[i] = data[i] > threshold ? 255 : 0;+  }+}++/**+ * Applies two thresholds to the image data by blackening pixels+ * below the lower threshold and whitening pixels above the upper threshold.+ * Pixels between the two thresholds are scaled to the range [0, 255].+ *+ * @param img_length_px Length of the image data in pixels.+ * @param data Pointer to the image data.+ * @param lower_threshold Every pixel below this value will be blackened.+ * @param upper_threshold Every pixel above this value will be whitened.+ *+ */+void apply_double_threshold(+  uint32_t img_length_px,+  uint8_t *data,+  uint8_t lower_threshold,+  uint8_t upper_threshold+) {+  for (uint32_t i = 0; i < img_length_px; i++) {+    if (data[i] < lower_threshold) {+      data[i] = 0;+    }+    else if (data[i] > upper_threshold) {+      data[i] = 255;+    }+    else {+      data[i] =+        (data[i] - lower_threshold) * 255 / (upper_threshold - lower_threshold);+    }+  }+}++/**+ * Apply Otsu's thresholding algorithm to the image data.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param use_double_threshold Whether to use double thresholding.+ * @param data Pointer to the pixel data.+ * @return Pointer to the monochrome image data.+ */+uint8_t *fcv_otsu_threshold_rgba(+  uint32_t width,+  uint32_t height,+  bool use_double_threshold,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height+  if (width > SIZE_MAX / height) {+    return NULL;+  }++  uint8_t *grayscale_img = fcv_rgba_to_grayscale(width, height, data);+  if (!grayscale_img) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;++  uint32_t histogram[256] = {0};+  for (uint32_t i = 0; i < img_length_px; i++) {+    histogram[grayscale_img[i]]++;+  }++  float histogram_norm[256] = {0};+  for (uint32_t i = 0; i < 256; i++) {+    histogram_norm[i] = (float)histogram[i] / img_length_px;+  }++  float global_mean = 0.0;++  for (uint32_t i = 0; i < 256; i++) {+    global_mean += i * histogram_norm[i];+  }++  float cumulative_sum = 0.0;+  float cumulative_mean = 0.0;+  float max_variance = 0.0;+  int32_t optimal_threshold = 0;++  for (uint32_t i = 0; i < 256; i++) {+    cumulative_sum += histogram_norm[i];+    cumulative_mean += i * histogram_norm[i];++    if (cumulative_sum == 0 || cumulative_sum == 1) {+      continue;+    }++    float mean1 = cumulative_mean / cumulative_sum;+    float mean2 = (global_mean - cumulative_mean) / (1 - cumulative_sum);++    float class_variance =+      cumulative_sum * (1 - cumulative_sum) * (mean1 - mean2) * (mean1 - mean2);++    if (class_variance > max_variance) {+      max_variance = class_variance;+      optimal_threshold = i;+    }+  }++  const int32_t threshold_range_offset = 16;++  if (use_double_threshold) {+    apply_double_threshold(+      img_length_px,+      grayscale_img,+      optimal_threshold - threshold_range_offset,+      optimal_threshold + threshold_range_offset+    );+  }+  else {+    fcv_apply_global_threshold(img_length_px, grayscale_img, optimal_threshold);+  }++  uint8_t *monochrome_data =+    fcv_single_to_multichannel(width, height, grayscale_img);++  free(grayscale_img);++  return monochrome_data;+}++/**+ * Apply gaussian blur to the image data.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the blurred image data.+ */+uint8_t *fcv_apply_gaussian_blur(+  uint32_t width,+  uint32_t height,+  double radius,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Validate radius+  if (radius < 0 || !isfinite(radius)) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;+  if (img_length_px > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = img_length_px * 4;++  if (radius == 0) {+    uint8_t *copy = malloc(img_length_byte);+    if (!copy) {+      return NULL;+    }+    return memcpy(copy, data, img_length_byte);+  }++  // Reject excessive radius to prevent excessive memory allocation+  if (radius > 1000) {+    return NULL;+  }++  uint8_t *blurred_data = malloc(img_length_byte);++  if (!blurred_data) { // Memory allocation failed+    return NULL;+  }++  uint32_t kernel_size = (uint32_t)(2 * radius + 1);+  float *kernel = malloc(kernel_size * sizeof(float));++  if (!kernel) { // Memory allocation failed+    free(blurred_data);+    return NULL;+  }++  float sigma = radius / 3.0;+  float sigma_sq = sigma * sigma;+  float two_sigma_sq = 2 * sigma_sq;+  float sqrt_two_pi_sigma = sqrt(2 * M_PI) * sigma;++  for (uint32_t i = 0; i < kernel_size; i++) {+    int32_t x = i - radius;+    kernel[i] = exp(-(x * x) / two_sigma_sq) / sqrt_two_pi_sigma;+  }++  // Apply the kernel in the horizontal direction+  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      float r_sum = 0.0;+      float g_sum = 0.0;+      float b_sum = 0.0;+      float weight_sum = 0.0;++      for (int32_t k = -radius; k <= radius; k++) {+        int32_t x_offset = x + k;+        if (x_offset < 0 || (uint32_t)x_offset >= width) {+          continue;+        }++        uint32_t img_index = y * width + x_offset;+        uint32_t img_rgba_index = img_index * 4;++        // Bounds check for array access+        if (img_rgba_index + 2 >= img_length_px * 4) {+          continue;+        }++        float weight = kernel[k + (int32_t)radius];+        weight_sum += weight;++        r_sum += data[img_rgba_index] * weight;+        g_sum += data[img_rgba_index + 1] * weight;+        b_sum += data[img_rgba_index + 2] * weight;+      }++      uint32_t rgba_index = (y * width + x) * 4;+      blurred_data[rgba_index] = r_sum / weight_sum;+      blurred_data[rgba_index + 1] = g_sum / weight_sum;+      blurred_data[rgba_index + 2] = b_sum / weight_sum;+      blurred_data[rgba_index + 3] = 255;+    }+  }++  // Create temporary buffer for vertical pass to avoid reading from buffer+  // being written to+  uint8_t *temp_data = malloc(img_length_px * 4);+  if (!temp_data) {+    free(blurred_data);+    free(kernel);+    return NULL;+  }++  // Copy horizontal blur result to temp buffer+  memcpy(temp_data, blurred_data, img_length_px * 4);++  // Apply the kernel in the vertical direction+  for (uint32_t x = 0; x < width; x++) {+    for (uint32_t y = 0; y < height; y++) {+      float r_sum = 0.0;+      float g_sum = 0.0;+      float b_sum = 0.0;+      float weight_sum = 0.0;++      for (int32_t k = -radius; k <= radius; k++) {+        int32_t y_offset = y + k;+        if (y_offset < 0 || (uint32_t)y_offset >= height) {+          continue;+        }++        uint32_t img_index = y_offset * width + x;+        uint32_t img_rgba_index = img_index * 4;++        float weight = kernel[k + (int32_t)radius];+        weight_sum += weight;++        r_sum += temp_data[img_rgba_index] * weight;+        g_sum += temp_data[img_rgba_index + 1] * weight;+        b_sum += temp_data[img_rgba_index + 2] * weight;+      }++      uint32_t rgba_index = (y * width + x) * 4;+      blurred_data[rgba_index] = r_sum / weight_sum;+      blurred_data[rgba_index + 1] = g_sum / weight_sum;+      blurred_data[rgba_index + 2] = b_sum / weight_sum;+      blurred_data[rgba_index + 3] = 255;+    }+  }++  free(temp_data);++  free(kernel);++  return blurred_data;+}++#include <time.h>++/**+ * Convert image to anti-aliased black and white.+ * 1. Convert the image to grayscale.+ * 2. Subtract blurred image from the original image to get the high+ * frequencies.+ * 3. Apply OTSU's threshold to get the optimal threshold.+ * 4. Apply the threshold + offset to get the anti-aliased image.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param use_double_threshold Whether to use double thresholding.+ * @param data Pointer to the pixel data.+ * @return Pointer to the blurred image data.+ */+uint8_t *fcv_bw_smart(+  uint32_t width,+  uint32_t height,+  bool use_double_threshold,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  uint8_t *grayscale_data = fcv_grayscale(width, height, data);+  if (!grayscale_data) {+    return NULL;+  }++  // Calculate blur radius dependent on image size+  // (Empirical formula after testing)+  double blurRadius = (sqrt((double)width * (double)height)) * 0.1;++  uint8_t *blurred_data =+    fcv_apply_gaussian_blur(width, height, blurRadius, grayscale_data);+  if (!blurred_data) {+    free(grayscale_data);+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    free(grayscale_data);+    free(blurred_data);+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;+  if (img_length_px > SIZE_MAX / 4) {+    free(grayscale_data);+    free(blurred_data);+    return NULL;+  }++  uint8_t *high_freq_data = malloc(img_length_px * 4);++  if (!high_freq_data) { // Memory allocation failed+    free((void *)grayscale_data);+    free((void *)blurred_data);+    return NULL;+  }++  // Subtract blurred image from the original image to get the high frequencies+  // and invert the high frequencies to get a white background.+  for (uint32_t i = 0; i < img_length_px; i++) {+    uint32_t rgba_idx = i * 4;+    int32_t high_freq_val =+      127 + grayscale_data[rgba_idx] - blurred_data[rgba_idx];++    // Clamp the value to [0, 255] to prevent overflow+    if (high_freq_val < 0) {+      high_freq_val = 0;+    }+    else if (high_freq_val > 255) {+      high_freq_val = 255;+    }++    high_freq_data[rgba_idx] = high_freq_val;     // R+    high_freq_data[rgba_idx + 1] = high_freq_val; // G+    high_freq_data[rgba_idx + 2] = high_freq_val; // B+    high_freq_data[rgba_idx + 3] = 255;           // A+  }++  free((void *)grayscale_data);+  free((void *)blurred_data);++  uint8_t *final_data = fcv_otsu_threshold_rgba(+    width,+    height,+    use_double_threshold,+    high_freq_data+  );++  free(high_freq_data);++  return final_data;+}++/**+ * Resize an image by given resize factors using bilinear interpolation.+ *+ * @param width Width of the input image.+ * @param height Height of the input image.+ * @param resize_x Horizontal resize factor (e.g., 2.0 for 2x, 0.5 for half).+ * @param resize_y Vertical resize factor (e.g., 2.0 for 2x, 0.5 for half).+ * @param out_width Pointer to store the output image width.+ * @param out_height Pointer to store the output image height.+ * @param data Pointer to the input pixel data.+ * @return Pointer to the resized image data.+ */+uint8_t *fcv_resize(+  uint32_t width,+  uint32_t height,+  double resize_x,+  double resize_y,+  uint32_t *out_width,+  uint32_t *out_height,+  uint8_t const *const data+) {+  if (!data || !out_width || !out_height) {+    return NULL;+  }++  if (width == 0 || height == 0) {+    return NULL;+  }++  if (resize_x <= 0.0 || resize_y <= 0.0 || !isfinite(resize_x) ||+      !isfinite(resize_y)) {+    return NULL;+  }++  // Check for overflow in output dimensions+  double new_width_d = (double)width * resize_x;+  double new_height_d = (double)height * resize_y;++  if (new_width_d > UINT32_MAX || new_height_d > UINT32_MAX) {+    return NULL;+  }++  *out_width = (uint32_t)new_width_d;+  *out_height = (uint32_t)new_height_d;++  if (*out_width == 0 || *out_height == 0) {+    return NULL;+  }++  // Check for overflow: out_width * out_height * 4+  if (*out_width > SIZE_MAX / *out_height) {+    return NULL;+  }+  size_t out_img_pixels = (size_t)(*out_width) * (*out_height);+  if (out_img_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t out_img_length = out_img_pixels * 4;++  uint8_t *resized_data = malloc(out_img_length);++  if (!resized_data) {+    return NULL;+  }++  for (uint32_t out_y = 0; out_y < *out_height; out_y++) {+    for (uint32_t out_x = 0; out_x < *out_width; out_x++) {+      if (resize_x < 1.0 || resize_y < 1.0) {+        double src_x = (out_x + 0.5) / resize_x - 0.5;+        double src_y = (out_y + 0.5) / resize_y - 0.5;++        double filter_size_x = 1.0 / resize_x;+        double filter_size_y = 1.0 / resize_y;++        double x_start = src_x - filter_size_x * 0.5;+        double y_start = src_y - filter_size_y * 0.5;+        double x_end = src_x + filter_size_x * 0.5;+        double y_end = src_y + filter_size_y * 0.5;++        int32_t ix_start = (int32_t)floor(x_start);+        int32_t iy_start = (int32_t)floor(y_start);+        int32_t ix_end = (int32_t)ceil(x_end);+        int32_t iy_end = (int32_t)ceil(y_end);++        if (ix_start < 0) {+          ix_start = 0;+        }+        if (iy_start < 0) {+          iy_start = 0;+        }+        if (ix_end > (int32_t)width) {+          ix_end = width;+        }+        if (iy_end > (int32_t)height) {+          iy_end = height;+        }++        double r_sum = 0.0, g_sum = 0.0, b_sum = 0.0;+        double total_weight = 0.0;++        for (int32_t sy = iy_start; sy < iy_end; sy++) {+          for (int32_t sx = ix_start; sx < ix_end; sx++) {+            double left = sx;+            double right = sx + 1;+            double top = sy;+            double bottom = sy + 1;++            double overlap_left = left > x_start ? left : x_start;+            double overlap_right = right < x_end ? right : x_end;+            double overlap_top = top > y_start ? top : y_start;+            double overlap_bottom = bottom < y_end ? bottom : y_end;++            if (overlap_right > overlap_left && overlap_bottom > overlap_top) {+              double weight =+                (overlap_right - overlap_left) * (overlap_bottom - overlap_top);+              total_weight += weight;++              uint32_t src_idx = (sy * width + sx) * 4;+              r_sum += data[src_idx] * weight;+              g_sum += data[src_idx + 1] * weight;+              b_sum += data[src_idx + 2] * weight;+            }+          }+        }++        if (total_weight > 0.0) {+          resized_data[(out_y * *out_width + out_x) * 4] =+            (uint8_t)(r_sum / total_weight + 0.5);+          resized_data[(out_y * *out_width + out_x) * 4 + 1] =+            (uint8_t)(g_sum / total_weight + 0.5);+          resized_data[(out_y * *out_width + out_x) * 4 + 2] =+            (uint8_t)(b_sum / total_weight + 0.5);+        }+        else {+          resized_data[(out_y * *out_width + out_x) * 4] = 0;+          resized_data[(out_y * *out_width + out_x) * 4 + 1] = 0;+          resized_data[(out_y * *out_width + out_x) * 4 + 2] = 0;+        }+        resized_data[(out_y * *out_width + out_x) * 4 + 3] = 255;+      }+      else {+        double src_x = (out_x + 0.5) / resize_x - 0.5;+        double src_y = (out_y + 0.5) / resize_y - 0.5;++        int32_t x0 = (int32_t)floor(src_x);+        int32_t y0 = (int32_t)floor(src_y);+        int32_t x1 = x0 + 1;+        int32_t y1 = y0 + 1;++        if (x0 < 0) {+          x0 = 0;+        }+        if (y0 < 0) {+          y0 = 0;+        }+        if (x1 >= (int32_t)width) {+          x1 = width - 1;+        }+        if (y1 >= (int32_t)height) {+          y1 = height - 1;+        }++        double dx = src_x - x0;+        double dy = src_y - y0;++        if (dx < 0) {+          dx = 0;+        }+        if (dy < 0) {+          dy = 0;+        }+        if (dx > 1) {+          dx = 1;+        }+        if (dy > 1) {+          dy = 1;+        }++        for (int32_t c = 0; c < 3; c++) {+          uint8_t p00 = data[(y0 * width + x0) * 4 + c];+          uint8_t p01 = data[(y0 * width + x1) * 4 + c];+          uint8_t p10 = data[(y1 * width + x0) * 4 + c];+          uint8_t p11 = data[(y1 * width + x1) * 4 + c];++          double interpolated = p00 * (1 - dx) * (1 - dy) ++                                p01 * dx * (1 - dy) + p10 * (1 - dx) * dy ++                                p11 * dx * dy;++          resized_data[(out_y * *out_width + out_x) * 4 + c] =+            (uint8_t)(interpolated + 0.5);+        }++        resized_data[(out_y * *out_width + out_x) * 4 + 3] = 255;+      }+    }+  }++  return resized_data;+}+// File: src/convert_to_binary.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "convert_to_binary.h"+#include "parse_hex_color.h"+#else+#include "flatcv.h"+#endif++/**+ * Convert an image to binary format based on foreground and background colors.+ * Return a grayscale image where the foreground color is white and the+ * background color is black.+ *+ * @param image_data Pointer to the input image data (RGBA format).+ * @param width Width of the image.+ * @param height Height of the image.+ * @param foreground_hex Hex color code for the foreground color (e.g.,+ * "FF0000").+ * @param background_hex Hex color code for the background color (e.g.,+ * "00FF00").+ */+uint8_t *fcv_convert_to_binary(+  const uint8_t *image_data,+  int32_t width,+  int32_t height,+  const char *foreground_hex,+  const char *background_hex+) {+  if (!image_data || !foreground_hex || !background_hex) {+    return NULL;+  }++  if (width <= 0 || height <= 0) {+    return NULL;+  }++  // Check for overflow: width * height (width and height already validated > 0)+  if ((size_t)width > SIZE_MAX / (size_t)height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * (size_t)height;++  uint8_t *result = malloc(num_pixels);+  if (!result) {+    return NULL;+  }++  uint8_t r_fg, g_fg, b_fg;+  fcv_parse_hex_color(foreground_hex, &r_fg, &g_fg, &b_fg);++  uint8_t r_bg, g_bg, b_bg;+  fcv_parse_hex_color(background_hex, &r_bg, &g_bg, &b_bg);++  for (size_t i = 0; i < num_pixels; i++) {+    size_t rgba_index = i * 4;+    uint8_t r = image_data[rgba_index];+    uint8_t g = image_data[rgba_index + 1];+    uint8_t b = image_data[rgba_index + 2];++    if (r == r_fg && g == g_fg && b == b_fg) {+      // Convert foreground color to white+      result[i] = 255;+    }+    else if (r == r_bg && g == g_bg && b == b_bg) {+      // Convert background color to black+      result[i] = 0;+    }+    else {+      // Default to black for unmatched colors+      result[i] = 0;+    }+  }++  return result;+}+// File: src/corner_detection.c+#include <assert.h>+#include <math.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "binary_closing_disk.h"+#include "conversion.h"+#include "convert_to_binary.h"+#include "corner_peaks.h"+#include "draw.h"+#include "foerstner_corner.h"+#include "parse_hex_color.h"+#include "perspectivetransform.h"+#include "sobel_edge_detection.h"+#include "sort_corners.h"+#include "watershed_segmentation.h"++#ifdef DEBUG_LOGGING+#define STB_IMAGE_WRITE_IMPLEMENTATION+#include "stb_image_write.h"+#endif+#else+#include "flatcv.h"+#endif++/** Count number of distinct colors in RGBA image+ * 1. Convert to grayscale+ * 2. Scale to 256x256 (save scale ratio for x and y)+ * 3. Blur image+ * 4. Create elevation map with Sobel+ * 5. Flatten elevation map at center seed to avoid being trapped in a local+ * minimum+ * 6. Add black border to flood from all directions at once+ * 7. `watershed(image=elevation_map, markers=markers)`+ *     Set center as marker for foreground basin and border for background basin+ * 8. Check region count equals 2+ * 9. Smooth result with binary closing+ * 10. Use Foerstner corner detector+ *     (Harris detector corners are shifted inwards)+ * 11. Sort corners+ * 12. Select 4 corners with the largest angle while maintaining their order+ * 13. Normalize corners based on scale ratio+ */+int32_t count_colors(const uint8_t *image, int32_t width, int32_t height) {+  assert(image != NULL);+  assert(width > 0);+  assert(height > 0);++  int32_t color_count = 0;+  int32_t *color_map = (int32_t *)calloc(256 * 256 * 256, sizeof(int32_t));+  if (!color_map) {+    fprintf(stderr, "Error: Failed to allocate memory for color map\n");+    return -1;+  }++  for (int32_t i = 0; i < width * height * 4; i += 4) {+    uint8_t r = image[i];+    uint8_t g = image[i + 1];+    uint8_t b = image[i + 2];+    uint8_t a = image[i + 3];++    // Skip fully transparent pixels+    if (a == 0) {+      continue;+    }++    // Create a unique color key+    int32_t color_key = (r << 16) | (g << 8) | b;++    // Increment the count for this color+    if (color_map[color_key] == 0) {+      color_count++;+    }+    color_map[color_key]++;+  }++  free(color_map);+  return color_count;+}++typedef struct {+  int32_t width;+  int32_t height;+  int32_t channels;+  const uint8_t *data;+} Image;++#ifdef DEBUG_LOGGING+void write_debug_img(Image img, const char *filename) {+  int32_t result = stbi_write_png(+    filename,+    img.width,+    img.height,+    img.channels,+    img.data,+    img.width * img.channels+  );+  if (result == 0) {+    fprintf(stderr, "Error: Failed to write debug image %s\n", filename);+  }+}+#endif++/**+ * Detect corners in the input image.+ */+Corners+fcv_detect_corners(const uint8_t *image, int32_t width, int32_t height) {+  assert(image != NULL);+  assert(width > 0);+  assert(height > 0);++  Corners default_corners = {+    .tl_x = 0,+    .tl_y = 0,+    .tr_x = width - 1,+    .tr_y = 0,+    .br_x = width - 1,+    .br_y = height - 1,+    .bl_x = 0,+    .bl_y = height - 1+  };++  // 1. Convert to grayscale+  uint8_t const *grayscale_image = fcv_grayscale(width, height, image);+  if (!grayscale_image) {+    fprintf(stderr, "Error: Failed to convert image to grayscale\n");+    return default_corners;+  }++  // 2. Resize image to 256x256+  uint32_t out_width = 256;+  uint32_t out_height = 256;+  uint8_t const *resized_image = fcv_resize(+    width,+    height,+    (double)out_width / width,+    (double)out_height / height,+    &out_width,+    &out_height,+    grayscale_image+  );+  free((void *)grayscale_image);+  if (!resized_image) {+    fprintf(stderr, "Error: Failed to resize image\n");+    exit(EXIT_FAILURE);+  }++#ifdef DEBUG_LOGGING+  Image out_img = {+    .width = out_width,+    .height = out_height,+    .channels = 4,+    .data = resized_image+  };+  write_debug_img(out_img, "temp_1_resized_image.png");+#endif++  // 3. Apply Gaussian blur+  uint8_t const *blurred_image =+    fcv_apply_gaussian_blur(out_width, out_height, 3.0, resized_image);+#ifndef DEBUG_LOGGING+  free((void *)resized_image);+#endif+  if (!blurred_image) {+    fprintf(stderr, "Error: Failed to apply Gaussian blur\n");+    exit(EXIT_FAILURE);+  }++  // 4. Create elevation map with Sobel edge detection+  uint8_t *elevation_map = (uint8_t *)+    fcv_sobel_edge_detection(out_width, out_height, 4, blurred_image);+  free((void *)blurred_image);+  if (!elevation_map) {+    fprintf(stderr, "Error: Failed to create elevation map with Sobel\n");+    exit(EXIT_FAILURE);+  }++  // 5. Flatten elevation map at center to not get trapped in a local minimum+  fcv_draw_disk(+    out_width,+    out_height,+    1, // Single channel grayscale+    "000000",+    24,+    out_width / 2.0,+    out_height / 2.0,+    elevation_map+  );++  // 6. Add black border to flood from all directions at once+  uint32_t bordered_width, bordered_height;+  uint8_t *bordered_elevation_map = fcv_add_border(+    out_width,+    out_height,+    1,        // Single channel grayscale+    "000000", // Black border+    1,        // Border width of 1 pixel+    elevation_map,+    &bordered_width,+    &bordered_height+  );++  free(elevation_map);+  if (!bordered_elevation_map) {+    fprintf(stderr, "Error: Failed to add black border to elevation map\n");+    exit(EXIT_FAILURE);+  }++#ifdef DEBUG_LOGGING+  out_img.data = bordered_elevation_map;+  out_img.channels = 1; // Grayscale+  out_img.width = bordered_width;+  out_img.height = bordered_height;+  write_debug_img(out_img, "temp_2_elevation_map.png");+#endif++  // 7. Perform watershed segmentation+  int32_t num_markers = 2;+  Point2D *markers = malloc(num_markers * sizeof(Point2D));+  if (!markers) {+    fprintf(stderr, "Error: Failed to allocate memory for markers\n");+    free((void *)bordered_elevation_map);+    exit(EXIT_FAILURE);+  }+  // Set center as foreground marker and upper left corner as background marker+  // Adjust for the added border (border width = 1)+  markers[0] = (Point2D){.x = bordered_width / 2.0, .y = bordered_height / 2.0};+  markers[1] = (Point2D){.x = 0, .y = 0};++  uint8_t *segmented_image_wide = fcv_watershed_segmentation(+    bordered_width,+    bordered_height,+    bordered_elevation_map,+    markers,+    num_markers,+    false // No boundaries+  );+  free((void *)bordered_elevation_map);+  free((void *)markers);++  // Remove 1 pixel border from segmented image+  uint8_t *segmented_image = malloc(out_width * out_height * 4);+  if (!segmented_image) {+    fprintf(stderr, "Error: Failed to allocate memory for segmented image\n");+    free((void *)segmented_image_wide);+    exit(EXIT_FAILURE);+  }+  for (uint32_t y = 1; y < bordered_height - 1; y++) {+    memcpy(+      &segmented_image[(y - 1) * out_width * 4],+      &segmented_image_wide[(y * bordered_width + 1) * 4],+      out_width * 4+    );+  }++#ifdef DEBUG_LOGGING+  out_img.width = out_width;+  out_img.height = out_height;+  out_img.channels = 4;+  out_img.data = segmented_image;+  write_debug_img(out_img, "temp_3_watershed.png");+#endif++  // 8. Check if the image has exactly 2 regions (foreground and background)+  int32_t region_count = count_colors(segmented_image, out_width, out_height);+  if (region_count != 2) {+    fprintf(stderr, "Error: Expected 2 regions, found %d\n", region_count);+    free((void *)segmented_image);+    exit(EXIT_FAILURE);+  }++  // Convert red pixel (marker 1 - foreground) to white+  // and green pixel (marker 2 - background) to black+  uint8_t *segmented_binary = fcv_convert_to_binary(+    segmented_image,+    out_width,+    out_height,+    "FF0000", // red -> white+    "00FF00"  // green -> black+  );+  free((void *)segmented_image);++#ifdef DEBUG_LOGGING+  out_img.data = segmented_binary;+  out_img.channels = 1;+  write_debug_img(out_img, "temp_4_watershed_binary.png");+#endif++  // 9. Smooth the result+  uint8_t *segmented_closed = fcv_binary_closing_disk(+    segmented_binary,+    out_width,+    out_height,+    12 // Closing radius+  );+  free((void *)segmented_binary);+  if (!segmented_closed) {+    fprintf(stderr, "Error: Failed to perform binary closing\n");+    exit(EXIT_FAILURE);+  }++#ifdef DEBUG_LOGGING+  out_img.data = segmented_closed;+  write_debug_img(out_img, "temp_5_segmented_closed.png");+#endif++  // 10. Find corners in the closed image+  uint8_t *corner_response = fcv_foerstner_corner(+    out_width,+    out_height,+    segmented_closed,+    1.5 // Sigma for Gaussian smoothing+  );+  free((void *)segmented_closed);+  if (!corner_response) {+    fprintf(stderr, "Error: Failed to compute corner response\n");+    exit(EXIT_FAILURE);+  }++  // Extract `w` channel (error ellipse size) for visualization+  uint8_t *w_channel = malloc(out_width * out_height);+  if (!w_channel) {+    fprintf(stderr, "Error: Failed to allocate memory for w channel\n");+    free((void *)corner_response);+    exit(EXIT_FAILURE);+  }++  for (uint32_t i = 0; i < out_width * out_height; i++) {+    w_channel[i] = corner_response[i * 2]; // Extract `w` channel+  }++#ifdef DEBUG_LOGGING+  out_img.data = w_channel;+  write_debug_img(out_img, "temp_6_corner_response.png");+#endif+  free(w_channel);++  // 10. Find corner peaks using thresholds+  // Gradually decrease thresholds until it finds at least 4 corners+  CornerPeaks *peaks = NULL;+  double accuracy_thresh = 0.4;+  double roundness_thresh = 0.2;+  const double thresh_decrement = 0.05;+  const double min_thresh = 0.01;++  do {+    if (peaks) {+      free(peaks->points);+      free(peaks);+    }++    peaks = fcv_corner_peaks(+      out_width,+      out_height,+      corner_response,+      16, // Minimum distance between peaks+      accuracy_thresh,+      roundness_thresh+    );++    if (!peaks || peaks->count < 4) {+      accuracy_thresh -= thresh_decrement;+      roundness_thresh -= thresh_decrement;++      // Prevent thresholds from going negative+      if (accuracy_thresh < min_thresh) {+        accuracy_thresh = min_thresh;+      }+      if (roundness_thresh < min_thresh) {+        roundness_thresh = min_thresh;+      }+    }++  } while (peaks && peaks->count < 4 &&+           (accuracy_thresh > min_thresh || roundness_thresh > min_thresh));+  free((void *)corner_response);+  if (!peaks) {+    fprintf(stderr, "Error: Failed to find corner peaks\n");+    exit(EXIT_FAILURE);+  }++  // First, sort all corners to get them in clockwise order+  Point2D *sorted_result = malloc(peaks->count * sizeof(Point2D));+  sort_corners(+    width,+    height,+    out_width,+    out_height,+    peaks->points,+    peaks->count,+    sorted_result+  );++  // Scale corners back to original image dimensions+  double scale_x = (double)width / out_width;+  double scale_y = (double)height / out_height;++  Corners sorted_corners;++  if (peaks->count > 4) {+    // Calculate angles for each corner using cross product method+    double *angles = malloc(peaks->count * sizeof(double));+    typedef struct {+      uint32_t index;+      double angle_abs;+    } AngleIndex;+    AngleIndex *angle_indices = malloc(peaks->count * sizeof(AngleIndex));++    for (uint32_t i = 0; i < peaks->count; i++) {+      // Get the three consecutive points for angle calculation+      Point2D prev = sorted_result[(i - 1 + peaks->count) % peaks->count];+      Point2D curr = sorted_result[i];+      Point2D next = sorted_result[(i + 1) % peaks->count];++      // Calculate vectors from current point+      double a_x = prev.x - curr.x;+      double a_y = prev.y - curr.y;+      double b_x = next.x - curr.x;+      double b_y = next.y - curr.y;++      // Calculate lengths of vectors+      double a_len = sqrt(a_x * a_x + a_y * a_y);+      double b_len = sqrt(b_x * b_x + b_y * b_y);++      if (a_len == 0 || b_len == 0) {+        angles[i] = 0;+      }+      else {+        // Calculate cross product (z-component in 2D)+        double cross_product = (a_x * b_y - a_y * b_x) / (a_len * b_len);+        // Clamp to [-1, 1] to avoid numerical errors+        if (cross_product > 1.0) {+          cross_product = 1.0;+        }+        if (cross_product < -1.0) {+          cross_product = -1.0;+        }++        // Calculate angle in degrees+        angles[i] = asin(cross_product) * 180.0 / M_PI;+      }++      angle_indices[i].index = i;+      angle_indices[i].angle_abs = fabs(angles[i]);+    }++    // Sort indices by descending angle magnitude+    for (uint32_t i = 0; i < peaks->count - 1; i++) {+      for (uint32_t j = 0; j < peaks->count - 1 - i; j++) {+        if (angle_indices[j].angle_abs < angle_indices[j + 1].angle_abs) {+          AngleIndex temp = angle_indices[j];+          angle_indices[j] = angle_indices[j + 1];+          angle_indices[j + 1] = temp;+        }+      }+    }++    // Take top 4 indices and sort them to maintain original clockwise order+    uint32_t top_4_indices[4];+    for (uint32_t i = 0; i < 4; i++) {+      top_4_indices[i] = angle_indices[i].index;+    }++    // Sort the top 4 indices to maintain original clockwise order+    for (uint32_t i = 0; i < 3; i++) {+      for (uint32_t j = 0; j < 3 - i; j++) {+        if (top_4_indices[j] > top_4_indices[j + 1]) {+          uint32_t temp = top_4_indices[j];+          top_4_indices[j] = top_4_indices[j + 1];+          top_4_indices[j + 1] = temp;+        }+      }+    }++    // Create final corners struct directly from selected corners in clockwise+    // order The first sorted corner becomes top-left, and we maintain clockwise+    // sequence+    sorted_corners =+      (Corners){.tl_x = sorted_result[top_4_indices[0]].x * scale_x,+                .tl_y = sorted_result[top_4_indices[0]].y * scale_y,+                .tr_x = sorted_result[top_4_indices[1]].x * scale_x,+                .tr_y = sorted_result[top_4_indices[1]].y * scale_y,+                .br_x = sorted_result[top_4_indices[2]].x * scale_x,+                .br_y = sorted_result[top_4_indices[2]].y * scale_y,+                .bl_x = sorted_result[top_4_indices[3]].x * scale_x,+                .bl_y = sorted_result[top_4_indices[3]].y * scale_y};++    free(angles);+    free(angle_indices);+  }+  else {+    // Use all available corners if we have 4 or fewer, already in clockwise+    // order+    sorted_corners =+      (Corners){.tl_x = sorted_result[0].x * scale_x,+                .tl_y = sorted_result[0].y * scale_y,+                .tr_x = sorted_result[1 % peaks->count].x * scale_x,+                .tr_y = sorted_result[1 % peaks->count].y * scale_y,+                .br_x = sorted_result[2 % peaks->count].x * scale_x,+                .br_y = sorted_result[2 % peaks->count].y * scale_y,+                .bl_x = sorted_result[3 % peaks->count].x * scale_x,+                .bl_y = sorted_result[3 % peaks->count].y * scale_y};+  }++  free(sorted_result);++#ifdef DEBUG_LOGGING+  // Print peaks for debugging+  fprintf(stderr, "Peaks:\n");+  fprintf(stderr, "%.0f,%.0f\n", sorted_corners.tl_x, sorted_corners.tl_y);+  fprintf(stderr, "%.0f,%.0f\n", sorted_corners.tr_x, sorted_corners.tr_y);+  fprintf(stderr, "%.0f,%.0f\n", sorted_corners.br_x, sorted_corners.br_y);+  fprintf(stderr, "%.0f,%.0f\n", sorted_corners.bl_x, sorted_corners.bl_y);+  // Draw corner peaks on the grayscale and resized image+  for (uint32_t i = 0; i < peaks->count; i++) {+    // Print peak coordinates for debugging+    fprintf(+      stderr,+      "Peak %u: (%.1f, %.1f)\n",+      i,+      peaks->points[i].x,+      peaks->points[i].y+    );+    fcv_draw_circle(+      out_width,+      out_height,+      4,        // GRAY as RGBA (each value is the same)+      "FFFF00", // Yellow color for corners+      2,        // Radius+      peaks->points[i].x,+      peaks->points[i].y,+      (uint8_t *)resized_image+    );+  }++  // Draw sorted corners with different colors to distinguish from peaks.+  // Reuse the scale variables from above.+  // (Inverted for drawing on resized image.)+  double draw_scale_x = (double)out_width / width;+  double draw_scale_y = (double)out_height / height;++  // Draw top-left corner in red+  fcv_draw_circle(+    out_width,+    out_height,+    4,+    "FF0000", // Red color for top-left+    3,        // Slightly larger radius to distinguish from peaks+    sorted_corners.tl_x * draw_scale_x,+    sorted_corners.tl_y * draw_scale_y,+    (uint8_t *)resized_image+  );++  // Draw top-right corner in green+  fcv_draw_circle(+    out_width,+    out_height,+    4,+    "00FF00", // Green color for top-right+    3,+    sorted_corners.tr_x * draw_scale_x,+    sorted_corners.tr_y * draw_scale_y,+    (uint8_t *)resized_image+  );++  // Draw bottom-right corner in blue+  fcv_draw_circle(+    out_width,+    out_height,+    4,+    "0000FF", // Blue color for bottom-right+    3,+    sorted_corners.br_x * draw_scale_x,+    sorted_corners.br_y * draw_scale_y,+    (uint8_t *)resized_image+  );++  // Draw bottom-left corner in magenta+  fcv_draw_circle(+    out_width,+    out_height,+    4,+    "FF00FF", // Magenta color for bottom-left+    3,+    sorted_corners.bl_x * draw_scale_x,+    sorted_corners.bl_y * draw_scale_y,+    (uint8_t *)resized_image+  );++  Image corners_img = {+    .width = out_width,+    .height = out_height,+    .channels = 4,+    .data = resized_image+  };+  write_debug_img(corners_img, "temp_7_corners.png");+  free((void *)resized_image);+#endif++  free(peaks);+  return sorted_corners;+}++/**+ * Wrapper function that allocates corners on heap for FFI.+ */+Corners *+fcv_detect_corners_ptr(const uint8_t *image, int32_t width, int32_t height) {+  Corners corners = fcv_detect_corners(image, width, height);+  Corners *result = malloc(sizeof(Corners));+  if (result == NULL) {+    return NULL;+  }+  *result = corners;+  return result;+}+// File: src/corner_peaks.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "1_types.h"+#include "conversion.h"+#include "corner_peaks.h"+#else+#include "flatcv.h"+#endif++static bool is_local_maximum(+  uint8_t const *data,+  uint32_t width,+  uint32_t height,+  uint32_t x,+  uint32_t y,+  uint32_t channel+) {+  if (x == 0 || y == 0 || x >= width - 1 || y >= height - 1) {+    return false;+  }++  uint32_t center_idx = (y * width + x) * 2 + channel;+  uint8_t center_val = data[center_idx];++  if (center_val == 0) {+    return false;+  }++  for (int32_t dy = -1; dy <= 1; dy++) {+    for (int32_t dx = -1; dx <= 1; dx++) {+      if (dx == 0 && dy == 0) {+        continue;+      }++      uint32_t neighbor_idx = ((y + dy) * width + (x + dx)) * 2 + channel;+      if (data[neighbor_idx] > center_val) {+        return false;+      }+    }+  }++  return true;+}++static double euclidean_distance(Point2D a, Point2D b) {+  double dx = a.x - b.x;+  double dy = a.y - b.y;+  return sqrt(dx * dx + dy * dy);+}++/** Detect corner peaks in the image.+ * This function finds local maxima in the corner response image+ * and filters them based on specified thresholds.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the 2 channel (w and q) corner response data.+ * @param min_distance Minimum distance between peaks.+ * @param accuracy_thresh Threshold for accuracy measure (w).+ * @param roundness_thresh Threshold for roundness measure (q).+ * @return Pointer to CornerPeaks structure containing detected peaks.+ */+CornerPeaks *fcv_corner_peaks(+  uint32_t width,+  uint32_t height,+  uint8_t const *data,+  uint32_t min_distance,+  double accuracy_thresh,+  double roundness_thresh+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Validate thresholds+  if (!isfinite(accuracy_thresh) || !isfinite(roundness_thresh)) {+    return NULL;+  }++  // Check for overflow: width * height+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;++  // Check for overflow in candidates allocation+  if (num_pixels > SIZE_MAX / sizeof(Point2D)) {+    return NULL;+  }++  Point2D *candidates = malloc(sizeof(Point2D) * num_pixels);+  if (!candidates) {+    return NULL;+  }++  uint32_t candidate_count = 0;++  for (uint32_t y = 1; y < height - 1; y++) {+    for (uint32_t x = 1; x < width - 1; x++) {+      uint32_t idx = (y * width + x) * 2;++      // Convert normalized values (0-255) back to 0-1 range+      double w = data[idx] / 255.0;     // accuracy measure+      double q = data[idx + 1] / 255.0; // roundness measure++      if (q > roundness_thresh && w > accuracy_thresh &&+          is_local_maximum(data, width, height, x, y, 0)) {+        candidates[candidate_count].x = (double)x;+        candidates[candidate_count].y = (double)y;+        candidate_count++;+      }+    }+  }++  if (candidate_count == 0) {+    free(candidates);+    CornerPeaks *result = malloc(sizeof(CornerPeaks));+    if (result) {+      result->points = NULL;+      result->count = 0;+    }+    return result;+  }++  bool *rejected = calloc(candidate_count, sizeof(bool));+  if (!rejected) {+    free(candidates);+    return NULL;+  }++  for (uint32_t i = 0; i < candidate_count; i++) {+    if (rejected[i]) {+      continue;+    }++    for (uint32_t j = i + 1; j < candidate_count; j++) {+      if (rejected[j]) {+        continue;+      }++      if (euclidean_distance(candidates[i], candidates[j]) < min_distance) {+        uint32_t i_idx =+          ((uint32_t)candidates[i].y * width + (uint32_t)candidates[i].x) * 2;+        uint32_t j_idx =+          ((uint32_t)candidates[j].y * width + (uint32_t)candidates[j].x) * 2;++        if (data[i_idx] >= data[j_idx]) {+          rejected[j] = true;+        }+        else {+          rejected[i] = true;+          break;+        }+      }+    }+  }++  uint32_t final_count = 0;+  for (uint32_t i = 0; i < candidate_count; i++) {+    if (!rejected[i]) {+      final_count++;+    }+  }++  Point2D *final_points = malloc(sizeof(Point2D) * final_count);+  if (!final_points) {+    free(candidates);+    free(rejected);+    return NULL;+  }++  uint32_t idx = 0;+  for (uint32_t i = 0; i < candidate_count; i++) {+    if (!rejected[i]) {+      final_points[idx] = candidates[i];+      idx++;+    }+  }++  free(candidates);+  free(rejected);++  CornerPeaks *result = malloc(sizeof(CornerPeaks));+  if (!result) {+    free(final_points);+    return NULL;+  }++  result->points = final_points;+  result->count = final_count;++  return result;+}+// File: src/crop.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "crop.h"+#else+#include "flatcv.h"+#endif++/**+ * Crop an image.+ *+ * @param width Width of the original image.+ * @param height Height of the original image.+ * @param channels Number of channels in the image.+ * @param data Pointer to the pixel data.+ * @param x The x-coordinate of the top-left corner of the crop area.+ * @param y The y-coordinate of the top-left corner of the crop area.+ * @param new_width The width of the crop area.+ * @param new_height The height of the crop area.+ * @return Pointer to the new cropped image data.+ */+uint8_t *fcv_crop(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  uint8_t const *const data,+  uint32_t x,+  uint32_t y,+  uint32_t new_width,+  uint32_t new_height+) {+  if (!data) {+    return NULL;+  }++  if (width == 0 || height == 0 || channels == 0 || new_width == 0 ||+      new_height == 0) {+    return NULL;+  }++  // Check for overflow in x + new_width and y + new_height+  if (x > UINT32_MAX - new_width || y > UINT32_MAX - new_height) {+    return NULL;+  }++  if (x + new_width > width || y + new_height > height) {+    fprintf(stderr, "Crop area is outside the original image bounds.\n");+    return NULL;+  }++  // Check for overflow in memory allocation: new_width * new_height * channels+  if (new_width > SIZE_MAX / new_height) {+    return NULL;+  }+  size_t num_pixels = (size_t)new_width * new_height;+  if (num_pixels > SIZE_MAX / channels) {+    return NULL;+  }+  size_t alloc_size = num_pixels * channels;++  uint8_t *cropped_data = malloc(alloc_size);+  if (!cropped_data) {+    fprintf(stderr, "Memory allocation failed for cropped image.\n");+    return NULL;+  }++  size_t row_bytes = (size_t)new_width * channels;+  for (uint32_t i = 0; i < new_height; ++i) {+    size_t src_index = ((size_t)(y + i) * width + x) * channels;+    size_t dst_index = (size_t)i * row_bytes;+    memcpy(cropped_data + dst_index, data + src_index, row_bytes);+  }++  return cropped_data;+}+// File: src/draw.c+#include <assert.h>+#include <limits.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "draw.h"+#include "parse_hex_color.h"+#else+#include "flatcv.h"+#endif++/**+ * Helper function to set a pixel to specified color (for circle drawing)+ */+void fcv_set_circle_pixel(+  uint8_t *data,+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  int32_t px,+  int32_t py,+  uint8_t r,+  uint8_t g,+  uint8_t b+) {+  if (!data || width == 0 || height == 0 || channels == 0) {+    return;+  }+  if (px >= 0 && px < (int32_t)width && py >= 0 && py < (int32_t)height) {+    size_t pixel_index = ((size_t)py * width + px) * channels;++    if (channels == 1) {+      // Grayscale: use luminance formula+      data[pixel_index] = (uint8_t)(0.299 * r + 0.587 * g + 0.114 * b);+    }+    else if (channels >= 3) {+      data[pixel_index] = r;     // R+      data[pixel_index + 1] = g; // G+      data[pixel_index + 2] = b; // B+      if (channels == 4) {+        data[pixel_index + 3] = 255; // A+      }+    }+  }+}++/**+ * Helper function to draw 8 symmetric points of a circle+ */+void fcv_draw_circle_points(+  uint8_t *data,+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  int32_t cx,+  int32_t cy,+  int32_t x,+  int32_t y,+  uint8_t r,+  uint8_t g,+  uint8_t b+) {+  fcv_set_circle_pixel(data, width, height, channels, cx + x, cy + y, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx - x, cy + y, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx + x, cy - y, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx - x, cy - y, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx + y, cy + x, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx - y, cy + x, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx + y, cy - x, r, g, b);+  fcv_set_circle_pixel(data, width, height, channels, cx - y, cy - x, r, g, b);+}++/**+ * Helper function to fill horizontal lines for a filled circle (disk)+ */+void fcv_fill_disk_lines(+  uint8_t *data,+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  int32_t cx,+  int32_t cy,+  int32_t x,+  int32_t y,+  uint8_t r,+  uint8_t g,+  uint8_t b+) {+  // Fill horizontal lines at y offsets+  for (int32_t i = cx - x; i <= cx + x; i++) {+    fcv_set_circle_pixel(data, width, height, channels, i, cy + y, r, g, b);+    fcv_set_circle_pixel(data, width, height, channels, i, cy - y, r, g, b);+  }++  // Fill horizontal lines at x offsets+  // (avoid duplicating center line when x == y)+  if (x != y) {+    for (int32_t i = cx - y; i <= cx + y; i++) {+      fcv_set_circle_pixel(data, width, height, channels, i, cy + x, r, g, b);+      fcv_set_circle_pixel(data, width, height, channels, i, cy - x, r, g, b);+    }+  }+}++/**+ * Draw a circle on an image using Bresenham's circle algorithm+ * (modifies in-place).+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param channels Number of channels in the image (1, 3, or 4).+ * @param hex_color Hex color code (e.g., "FF0000" for red, "#00FF00" for+ * green).+ * @param radius Radius of the circle.+ * @param center_x X coordinate of the circle center.+ * @param center_y Y coordinate of the circle center.+ * @param data Pointer to the pixel data (modified in-place).+ */+void fcv_draw_circle(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  const char *hex_color,+  double radius,+  double center_x,+  double center_y,+  uint8_t *data+) {+  if (!data || !hex_color || width == 0 || height == 0 || channels == 0) {+    return;+  }++  // Validate radius and center coordinates+  if (radius < 0 || !isfinite(radius) || !isfinite(center_x) ||+      !isfinite(center_y)) {+    return;+  }++  // Reject excessive radius to prevent integer overflow in Bresenham algorithm+  if (radius > INT32_MAX / 2) {+    return;+  }++  // Parse color+  uint8_t r, g, b;+  fcv_parse_hex_color(hex_color, &r, &g, &b);++  int32_t radius_int32_t = (int32_t)radius;+  int32_t cx = (int32_t)center_x;+  int32_t cy = (int32_t)center_y;++  // Bresenham's circle algorithm+  int32_t x = 0;+  int32_t y = radius_int32_t;+  int32_t d = 3 - 2 * radius_int32_t;++  // Initial points+  fcv_draw_circle_points(data, width, height, channels, cx, cy, x, y, r, g, b);++  while (y >= x) {+    x++;+    if (d > 0) {+      y--;+      d = d + 4 * (x - y) + 10;+    }+    else {+      d = d + 4 * x + 6;+    }+    fcv_draw_circle_points(+      data,+      width,+      height,+      channels,+      cx,+      cy,+      x,+      y,+      r,+      g,+      b+    );+  }+}++/**+ * Draw a filled circle (disk) on an image using Bresenham's circle+ * algorithm (modifies in-place).+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param channels Number of channels in the image (1, 3, or 4).+ * @param hex_color Hex color code (e.g., "FF0000" for red, "#00FF00" for+ * green).+ * @param radius Radius of the disk.+ * @param center_x X coordinate of the disk center.+ * @param center_y Y coordinate of the disk center.+ * @param data Pointer to the pixel data (modified in-place).+ */+void fcv_draw_disk(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  const char *hex_color,+  double radius,+  double center_x,+  double center_y,+  uint8_t *data+) {+  if (!data || !hex_color || width == 0 || height == 0 || channels == 0) {+    return;+  }++  // Validate radius and center coordinates+  if (radius < 0 || !isfinite(radius) || !isfinite(center_x) ||+      !isfinite(center_y)) {+    return;+  }++  // Reject excessive radius to prevent integer overflow in Bresenham algorithm+  if (radius > INT32_MAX / 2) {+    return;+  }++  // Parse color+  uint8_t r, g, b;+  fcv_parse_hex_color(hex_color, &r, &g, &b);++  int32_t radius_int32_t = (int32_t)radius;+  int32_t cx = (int32_t)center_x;+  int32_t cy = (int32_t)center_y;++  // Bresenham's circle algorithm for filled disk+  int32_t x = 0;+  int32_t y = radius_int32_t;+  int32_t d = 3 - 2 * radius_int32_t;++  // Initial filled lines+  fcv_fill_disk_lines(data, width, height, channels, cx, cy, x, y, r, g, b);++  while (y >= x) {+    x++;+    if (d > 0) {+      y--;+      d = d + 4 * (x - y) + 10;+    }+    else {+      d = d + 4 * x + 6;+    }+    fcv_fill_disk_lines(data, width, height, channels, cx, cy, x, y, r, g, b);+  }+}++/**+ * Add a colored border around an image.+ * Returns a new image with the border added.+ *+ * @param width Width of the input image.+ * @param height Height of the input image.+ * @param channels Number of channels in the image (1, 3, or 4).+ * @param hex_color Hex color code for the border (e.g., "FF0000" for red).+ * @param border_width Width of the border in pixels.+ * @param input_data Pointer to the input pixel data.+ * @param output_width Pointer to store the output image width.+ * @param output_height Pointer to store the output image height.+ * @return Pointer to the new image data with border, or NULL on failure.+ */+uint8_t *fcv_add_border(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  const char *hex_color,+  uint32_t border_width,+  uint8_t *input_data,+  uint32_t *output_width,+  uint32_t *output_height+) {+  if (!input_data || !output_width || !output_height || !hex_color ||+      border_width == 0 || channels == 0) {+    return NULL;+  }++  if (width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow in output dimensions+  if (border_width > (UINT32_MAX - width) / 2 ||+      border_width > (UINT32_MAX - height) / 2) {+    return NULL;+  }++  // Parse border color+  uint8_t r, g, b;+  fcv_parse_hex_color(hex_color, &r, &g, &b);++  // Calculate output dimensions+  *output_width = width + 2 * border_width;+  *output_height = height + 2 * border_width;++  // Check for overflow in memory allocation: output_width * output_height *+  // channels+  if (*output_width > SIZE_MAX / *output_height) {+    return NULL;+  }+  size_t num_pixels = (size_t)(*output_width) * (*output_height);+  if (num_pixels > SIZE_MAX / channels) {+    return NULL;+  }+  size_t output_size = num_pixels * channels;++  uint8_t *output_data = malloc(output_size);+  if (!output_data) {+    return NULL;+  }++  // Fill the entire output image with border color+  for (uint32_t y = 0; y < *output_height; y++) {+    for (uint32_t x = 0; x < *output_width; x++) {+      uint32_t pixel_index = (y * *output_width + x) * channels;++      if (channels == 1) {+        // Grayscale: use luminance formula+        output_data[pixel_index] = (uint8_t)(0.299 * r + 0.587 * g + 0.114 * b);+      }+      else if (channels >= 3) {+        output_data[pixel_index] = r;     // R+        output_data[pixel_index + 1] = g; // G+        output_data[pixel_index + 2] = b; // B+        if (channels == 4) {+          output_data[pixel_index + 3] = 255; // A+        }+      }+    }+  }++  // Copy the original image into the center+  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      uint32_t src_index = (y * width + x) * channels;+      uint32_t dst_index =+        ((y + border_width) * *output_width + (x + border_width)) * channels;++      for (uint32_t c = 0; c < channels; c++) {+        output_data[dst_index + c] = input_data[src_index + c];+      }+    }+  }++  return output_data;+}+// File: src/exif.c+#include <stdint.h>+#include <stdio.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "exif.h"+#endif++static uint16_t read_u16(FILE *file, int is_little_endian) {+  uint8_t buf[2];+  if (fread(buf, 1, 2, file) != 2) {+    return 0;+  }+  if (is_little_endian) {+    return (uint16_t)(buf[0] | (buf[1] << 8));+  }+  return (uint16_t)((buf[0] << 8) | buf[1]);+}++static uint32_t read_u32(FILE *file, int is_little_endian) {+  uint8_t buf[4];+  if (fread(buf, 1, 4, file) != 4) {+    return 0;+  }+  if (is_little_endian) {+    return (uint32_t)(buf[0] | (buf[1] << 8) | (buf[2] << 16) | (buf[3] << 24));+  }+  return (uint32_t)((buf[0] << 24) | (buf[1] << 16) | (buf[2] << 8) | buf[3]);+}++int32_t fcv_get_exif_orientation(char const *filename) {+  FILE *file = fopen(filename, "rb");+  if (!file) {+    return 1;+  }++  uint8_t header[2];+  if (fread(header, 1, 2, file) != 2 || header[0] != 0xFF ||+      header[1] != 0xD8) {+    fclose(file);+    return 1;+  }++  uint8_t marker[2];+  while (fread(marker, 1, 2, file) == 2) {+    if (marker[0] != 0xFF) {+      break;+    }++    if (marker[1] == 0xD9 || marker[1] == 0xDA) { // EOI or SOS+      break;+    }++    uint16_t size = read_u16(file, 0); // Markers are always big-endian++    if (marker[1] == 0xE1) { // APP1 - EXIF+      uint8_t exif_header[6];+      if (fread(exif_header, 1, 6, file) != 6) {+        break;+      }++      if (memcmp(exif_header, "Exif\0\0", 6) != 0) {+        fseek(file, size - 2 - 6, SEEK_CUR);+        continue;+      }++      long tiff_start = ftell(file);+      uint8_t tiff_header[2];+      if (fread(tiff_header, 1, 2, file) != 2) {+        break;+      }++      int is_little_endian = (tiff_header[0] == 'I' && tiff_header[1] == 'I');+      int is_big_endian = (tiff_header[0] == 'M' && tiff_header[1] == 'M');++      if (!is_little_endian && !is_big_endian) {+        break;+      }++      uint16_t magic = read_u16(file, is_little_endian);+      if (magic != 42) {+        break;+      }++      uint32_t ifd_offset = read_u32(file, is_little_endian);++      fseek(file, tiff_start + ifd_offset, SEEK_SET);++      uint16_t num_entries = read_u16(file, is_little_endian);++      for (uint16_t i = 0; i < num_entries; i++) {+        uint16_t tag = read_u16(file, is_little_endian);+        (void)read_u16(file, is_little_endian); // Skip type+        (void)read_u32(file, is_little_endian); // Skip count++        if (tag == 0x0112) { // Orientation+          uint16_t orientation = read_u16(file, is_little_endian);+          fclose(file);+          return (int32_t)orientation;+        }+        else {+          fseek(file, 4, SEEK_CUR); // Skip value/offset+        }+      }+      break; // IFD0 done+    }+    else {+      fseek(file, size - 2, SEEK_CUR);+    }+  }++  fclose(file);+  return 1;+}+// File: src/extract_document.c+#include <math.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>++#ifndef FLATCV_AMALGAMATION+#include "perspectivetransform.h"+#else+#include "flatcv.h"+#endif++uint8_t *fcv_extract_document(+  uint32_t width,+  uint32_t height,+  uint8_t const *const data,+  uint32_t output_width,+  uint32_t output_height+) {+  if (!data) {+    fprintf(stderr, "Error: NULL image data\n");+    return NULL;+  }++  // Step 1: Detect corners of the document+  Corners detected_corners = fcv_detect_corners(data, width, height);++  // Step 2: Define destination corners (rectangular output)+  Corners dst_corners = {+    .tl_x = 0.0,+    .tl_y = 0.0,+    .tr_x = (double)(output_width - 1),+    .tr_y = 0.0,+    .br_x = (double)(output_width - 1),+    .br_y = (double)(output_height - 1),+    .bl_x = 0.0,+    .bl_y = (double)(output_height - 1)+  };++  // Step 3: Calculate perspective transformation matrix+  // Note: We need to map FROM destination TO source for the inverse transform+  Matrix3x3 *transform_matrix =+    fcv_calculate_perspective_transform(&dst_corners, &detected_corners);++  if (!transform_matrix) {+    fprintf(stderr, "Error: Failed to calculate perspective transform\n");+    return NULL;+  }++  // Step 4: Apply the transformation+  uint8_t *result = fcv_apply_matrix_3x3(+    width,+    height,+    (uint8_t *)data,+    output_width,+    output_height,+    transform_matrix+  );++  // Free the transformation matrix if it's not the identity matrix+  if (transform_matrix &&+      !(transform_matrix->m00 == 1.0 && transform_matrix->m01 == 0.0 &&+        transform_matrix->m02 == 0.0 && transform_matrix->m10 == 0.0 &&+        transform_matrix->m11 == 1.0 && transform_matrix->m12 == 0.0 &&+        transform_matrix->m20 == 0.0 && transform_matrix->m21 == 0.0 &&+        transform_matrix->m22 == 1.0)) {+    free(transform_matrix);+  }++  return result;+}++uint8_t *fcv_extract_document_auto(+  uint32_t width,+  uint32_t height,+  uint8_t const *const data,+  uint32_t *output_width,+  uint32_t *output_height+) {+  if (!data || !output_width || !output_height) {+    fprintf(stderr, "Error: NULL parameters\n");+    return NULL;+  }++  // Step 1: Detect corners of the document+  Corners detected_corners = fcv_detect_corners(data, width, height);++  // Step 2: Calculate optimal output dimensions based on detected corners+  // Calculate distances between corners to determine aspect ratio+  double top_width = sqrt(+    pow(detected_corners.tr_x - detected_corners.tl_x, 2) ++    pow(detected_corners.tr_y - detected_corners.tl_y, 2)+  );+  double bottom_width = sqrt(+    pow(detected_corners.br_x - detected_corners.bl_x, 2) ++    pow(detected_corners.br_y - detected_corners.bl_y, 2)+  );+  double left_height = sqrt(+    pow(detected_corners.bl_x - detected_corners.tl_x, 2) ++    pow(detected_corners.bl_y - detected_corners.tl_y, 2)+  );+  double right_height = sqrt(+    pow(detected_corners.br_x - detected_corners.tr_x, 2) ++    pow(detected_corners.br_y - detected_corners.tr_y, 2)+  );++  // Use maximum dimensions to preserve detail+  double max_width = fmax(top_width, bottom_width);+  double max_height = fmax(left_height, right_height);++  // Round to reasonable dimensions+  *output_width = (uint32_t)(max_width + 0.5);+  *output_height = (uint32_t)(max_height + 0.5);++  // Ensure minimum reasonable size+  if (*output_width < 100) {+    *output_width = 100;+  }+  if (*output_height < 100) {+    *output_height = 100;+  }++  // Step 3: Define destination corners (rectangular output)+  Corners dst_corners = {+    .tl_x = 0.0,+    .tl_y = 0.0,+    .tr_x = (double)(*output_width - 1),+    .tr_y = 0.0,+    .br_x = (double)(*output_width - 1),+    .br_y = (double)(*output_height - 1),+    .bl_x = 0.0,+    .bl_y = (double)(*output_height - 1)+  };++  // Step 4: Calculate perspective transformation matrix+  // Note: We need to map FROM destination TO source for the inverse transform+  Matrix3x3 *transform_matrix =+    fcv_calculate_perspective_transform(&dst_corners, &detected_corners);++  if (!transform_matrix) {+    fprintf(stderr, "Error: Failed to calculate perspective transform\n");+    return NULL;+  }++  // Step 5: Apply the transformation+  uint8_t *result = fcv_apply_matrix_3x3(+    width,+    height,+    (uint8_t *)data,+    *output_width,+    *output_height,+    transform_matrix+  );++  // Free the transformation matrix if it's not the identity matrix+  if (transform_matrix &&+      !(transform_matrix->m00 == 1.0 && transform_matrix->m01 == 0.0 &&+        transform_matrix->m02 == 0.0 && transform_matrix->m10 == 0.0 &&+        transform_matrix->m11 == 1.0 && transform_matrix->m12 == 0.0 &&+        transform_matrix->m20 == 0.0 && transform_matrix->m21 == 0.0 &&+        transform_matrix->m22 == 1.0)) {+    free(transform_matrix);+  }++  return result;+}+// File: src/flip.c+#include <stdint.h>+#include <stdlib.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "flip.h"+#else+#include "flatcv.h"+#endif++/**+ * Flip an image horizontally (mirror along vertical axis).+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the flipped image data.+ */+uint8_t *+fcv_flip_x(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *flipped_data = malloc(img_length_byte);++  if (!flipped_data) { // Memory allocation failed+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index = ((size_t)y * width + (width - 1 - x)) * 4;++      flipped_data[dst_index] = data[src_index];         // R+      flipped_data[dst_index + 1] = data[src_index + 1]; // G+      flipped_data[dst_index + 2] = data[src_index + 2]; // B+      flipped_data[dst_index + 3] = data[src_index + 3]; // A+    }+  }++  return flipped_data;+}++/**+ * Flip an image vertically (mirror along horizontal axis).+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the flipped image data.+ */+uint8_t *+fcv_flip_y(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *flipped_data = malloc(img_length_byte);++  if (!flipped_data) { // Memory allocation failed+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index = ((size_t)(height - 1 - y) * width + x) * 4;++      flipped_data[dst_index] = data[src_index];         // R+      flipped_data[dst_index + 1] = data[src_index + 1]; // G+      flipped_data[dst_index + 2] = data[src_index + 2]; // B+      flipped_data[dst_index + 3] = data[src_index + 3]; // A+    }+  }++  return flipped_data;+}++/**+ * Transpose an image (flip along main diagonal).+ */+uint8_t *+fcv_transpose(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *transposed_data = malloc(img_length_byte);++  if (!transposed_data) {+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index = ((size_t)x * height + y) * 4;++      transposed_data[dst_index] = data[src_index];+      transposed_data[dst_index + 1] = data[src_index + 1];+      transposed_data[dst_index + 2] = data[src_index + 2];+      transposed_data[dst_index + 3] = data[src_index + 3];+    }+  }++  return transposed_data;+}++/**+ * Transverse an image (flip along anti-diagonal).+ */+uint8_t *+fcv_transverse(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *transposed_data = malloc(img_length_byte);++  if (!transposed_data) {+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index =+        ((size_t)(width - 1 - x) * height + (height - 1 - y)) * 4;++      transposed_data[dst_index] = data[src_index];+      transposed_data[dst_index + 1] = data[src_index + 1];+      transposed_data[dst_index + 2] = data[src_index + 2];+      transposed_data[dst_index + 3] = data[src_index + 3];+    }+  }++  return transposed_data;+}+// File: src/foerstner_corner.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "conversion.h"+#include "perspectivetransform.h"+#else+#include "flatcv.h"+#endif++/** Implementation of the Foerstner corner measure response image.+ * Expected input is a grayscale image.+ */+uint8_t *fcv_foerstner_corner(+  uint32_t width,+  uint32_t height,+  uint8_t const *const gray_data,+  double sigma+) {+  if (!gray_data || width == 0 || height == 0) {+    return NULL;+  }++  // Allocate memory for gradient images+  double *grad_x = calloc(width * height, sizeof(double));+  double *grad_y = calloc(width * height, sizeof(double));++  if (!grad_x || !grad_y) {+    free(grad_x);+    free(grad_y);+    return NULL;+  }++  // Compute image gradients using Sobel-like operators+  for (uint32_t y = 1; y < height - 1; y++) {+    for (uint32_t x = 1; x < width - 1; x++) {+      uint32_t idx = y * width + x;++      // Sobel X kernel: [-1 0 1; -2 0 2; -1 0 1]+      double gx = 0.0;+      gx += -1.0 * gray_data[(y - 1) * width + (x - 1)];+      gx += 1.0 * gray_data[(y - 1) * width + (x + 1)];+      gx += -2.0 * gray_data[y * width + (x - 1)];+      gx += 2.0 * gray_data[y * width + (x + 1)];+      gx += -1.0 * gray_data[(y + 1) * width + (x - 1)];+      gx += 1.0 * gray_data[(y + 1) * width + (x + 1)];++      // Sobel Y kernel: [-1 -2 -1; 0 0 0; 1 2 1]+      double gy = 0.0;+      gy += -1.0 * gray_data[(y - 1) * width + (x - 1)];+      gy += -2.0 * gray_data[(y - 1) * width + x];+      gy += -1.0 * gray_data[(y - 1) * width + (x + 1)];+      gy += 1.0 * gray_data[(y + 1) * width + (x - 1)];+      gy += 2.0 * gray_data[(y + 1) * width + x];+      gy += 1.0 * gray_data[(y + 1) * width + (x + 1)];++      grad_x[idx] = gx / 8.0; // Normalize+      grad_y[idx] = gy / 8.0;+    }+  }++  // Compute gradient products for structure tensor+  double *Axx = malloc(width * height * sizeof(double));+  double *Axy = malloc(width * height * sizeof(double));+  double *Ayy = malloc(width * height * sizeof(double));++  if (!Axx || !Axy || !Ayy) {+    free(grad_x);+    free(grad_y);+    free(Axx);+    free(Axy);+    free(Ayy);+    return NULL;+  }++  for (uint32_t i = 0; i < width * height; i++) {+    Axx[i] = grad_x[i] * grad_x[i];+    Axy[i] = grad_x[i] * grad_y[i];+    Ayy[i] = grad_y[i] * grad_y[i];+  }++  // Apply Gaussian smoothing to structure tensor elements+  // For simplicity, use a basic box filter approximation when sigma is small+  int32_t kernel_size = (int32_t)(3 * sigma) | 1; // Ensure odd size+  if (kernel_size < 3) {+    kernel_size = 3;+  }+  int32_t half_kernel = kernel_size / 2;++  double *Axx_smooth = calloc(width * height, sizeof(double));+  double *Axy_smooth = calloc(width * height, sizeof(double));+  double *Ayy_smooth = calloc(width * height, sizeof(double));++  if (!Axx_smooth || !Axy_smooth || !Ayy_smooth) {+    free(grad_x);+    free(grad_y);+    free(Axx);+    free(Axy);+    free(Ayy);+    free(Axx_smooth);+    free(Axy_smooth);+    free(Ayy_smooth);+    return NULL;+  }++  // Simple box filter smoothing+  for (uint32_t y = half_kernel; y < height - half_kernel; y++) {+    for (uint32_t x = half_kernel; x < width - half_kernel; x++) {+      uint32_t idx = y * width + x;+      double sum_xx = 0.0, sum_xy = 0.0, sum_yy = 0.0;+      int32_t count = 0;++      for (int32_t ky = -half_kernel; ky <= half_kernel; ky++) {+        for (int32_t kx = -half_kernel; kx <= half_kernel; kx++) {+          int32_t sample_idx = (y + ky) * width + (x + kx);+          sum_xx += Axx[sample_idx];+          sum_xy += Axy[sample_idx];+          sum_yy += Ayy[sample_idx];+          count++;+        }+      }++      Axx_smooth[idx] = sum_xx / count;+      Axy_smooth[idx] = sum_xy / count;+      Ayy_smooth[idx] = sum_yy / count;+    }+  }++  // Compute Foerstner measures w and q+  uint8_t *result = malloc(width * height * 2); // 2 channels: w, q+  if (!result) {+    free(grad_x);+    free(grad_y);+    free(Axx);+    free(Axy);+    free(Ayy);+    free(Axx_smooth);+    free(Axy_smooth);+    free(Ayy_smooth);+    return NULL;+  }++  double max_w = 0.0, max_q = 0.0;+  double *w_values = malloc(width * height * sizeof(double));+  double *q_values = malloc(width * height * sizeof(double));+  if (!w_values || !q_values) {+    free(result);+    free(w_values);+    free(q_values);+    free(grad_x);+    free(grad_y);+    free(Axx);+    free(Axy);+    free(Ayy);+    free(Axx_smooth);+    free(Axy_smooth);+    free(Ayy_smooth);+    return NULL;+  }++  // First pass: compute w and q values and find maximum for normalization+  for (uint32_t i = 0; i < width * height; i++) {+    double det_A =+      Axx_smooth[i] * Ayy_smooth[i] - Axy_smooth[i] * Axy_smooth[i];+    double trace_A = Axx_smooth[i] + Ayy_smooth[i];++    double w = 0.0, q = 0.0;++    if (fabs(trace_A) > 1e-10) { // Avoid division by zero+      w = det_A / trace_A;+      q = 4.0 * det_A / (trace_A * trace_A);+    }++    w_values[i] = fmax(0.0, w); // Ensure non-negative+    q_values[i] = fmax(0.0, q);++    if (w_values[i] > max_w) {+      max_w = w_values[i];+    }+    if (q_values[i] > max_q) {+      max_q = q_values[i];+    }+  }++  // Second pass: normalize and convert to bytes+  for (uint32_t i = 0; i < width * height; i++) {+    uint8_t w_byte = 0, q_byte = 0;++    if (max_w > 0.0) {+      w_byte = (uint8_t)(255.0 * w_values[i] / max_w);+    }+    if (max_q > 0.0) {+      q_byte = (uint8_t)(255.0 * q_values[i] / max_q);+    }++    result[i * 2] = w_byte;     // w measure+    result[i * 2 + 1] = q_byte; // q measure+  }++  free(w_values);+  free(q_values);++  // Cleanup+  free(grad_x);+  free(grad_y);+  free(Axx);+  free(Axy);+  free(Ayy);+  free(Axx_smooth);+  free(Axy_smooth);+  free(Ayy_smooth);++  return result;+}+// File: src/histogram.c+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "histogram.h"+#else+#include "flatcv.h"+#endif++/**+ * Generate a histogram visualization image from input image data.+ * For grayscale images, creates a single histogram.+ * For RGB(A) images, creates overlapping histograms for each channel.+ *+ * @param width Width of the input image.+ * @param height Height of the input image.+ * @param channels Number of channels in the input image (1, 3, or 4).+ * @param data Pointer to the input pixel data.+ * @param out_width Pointer to store the output histogram width.+ * @param out_height Pointer to store the output histogram height.+ * @return Pointer to the histogram image data (RGBA format).+ */+uint8_t *fcv_generate_histogram(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  uint8_t const *const data,+  uint32_t *out_width,+  uint32_t *out_height+) {+  if (!data || !out_width || !out_height) {+    return NULL;+  }++  if (width == 0 || height == 0 || channels == 0) {+    return NULL;+  }++  // Check for overflow in img_size calculation+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_size = (size_t)width * height;++  // Validate that we can access all pixels+  if (img_size > SIZE_MAX / channels) {+    return NULL;+  }++  // Output histogram dimensions+  const uint32_t hist_width = 256;  // One pixel per brightness value+  const uint32_t hist_height = 200; // Fixed height for visualization++  *out_width = hist_width;+  *out_height = hist_height;++  // Output size is fixed: 256 * 200 * 4 = 204800 bytes (no overflow possible)+  size_t output_size = (size_t)hist_width * hist_height * 4;++  // Allocate output image (initialized to black background)+  uint8_t *output = calloc(output_size, sizeof(uint8_t));+  if (!output) {+    return NULL;+  }++  // Initialize with black background and full alpha+  for (uint32_t i = 0; i < (*out_width) * (*out_height); i++) {+    output[i * 4 + 3] = 255; // Alpha channel+  }++  // Count histograms for each channel+  uint32_t hist_r[256] = {0};+  uint32_t hist_g[256] = {0};+  uint32_t hist_b[256] = {0};+  uint32_t hist_gray[256] = {0};++  bool is_grayscale = (channels == 1);+  if (channels == 4) {+    // Check if it's actually grayscale (all RGB channels equal)+    is_grayscale = true;+    for (size_t i = 0; i < img_size && is_grayscale; i++) {+      size_t idx = i * 4;+      if (data[idx] != data[idx + 1] || data[idx] != data[idx + 2]) {+        is_grayscale = false;+      }+    }+  }++  if (is_grayscale) {+    // Single channel or grayscale image+    for (size_t i = 0; i < img_size; i++) {+      uint8_t value;+      if (channels == 1) {+        value = data[i];+      }+      else {+        // Use the red channel for grayscale RGBA+        value = data[i * 4];+      }+      hist_gray[value]++;+    }+  }+  else {+    // RGB or RGBA image - count each channel+    for (size_t i = 0; i < img_size; i++) {+      size_t idx = i * channels;+      hist_r[data[idx]]++;+      if (channels >= 3) {+        hist_g[data[idx + 1]]++;+        hist_b[data[idx + 2]]++;+      }+    }+  }++  // Find maximum count for scaling+  uint32_t max_count = 0;+  if (is_grayscale) {+    for (int i = 0; i < 256; i++) {+      if (hist_gray[i] > max_count) {+        max_count = hist_gray[i];+      }+    }+  }+  else {+    for (int i = 0; i < 256; i++) {+      if (hist_r[i] > max_count) {+        max_count = hist_r[i];+      }+      if (channels >= 3) {+        if (hist_g[i] > max_count) {+          max_count = hist_g[i];+        }+        if (hist_b[i] > max_count) {+          max_count = hist_b[i];+        }+      }+    }+  }++  if (max_count == 0) {+    return output; // Return black image if no data+  }++  // Draw histogram bars+  for (int x = 0; x < 256; x++) {+    if (is_grayscale) {+      // Draw white histogram for grayscale+      uint32_t bar_height = (hist_gray[x] * hist_height) / max_count;+      for (uint32_t y = 0; y < bar_height; y++) {+        uint32_t output_x = x;+        uint32_t output_y = hist_height - 1 - y;+        uint32_t output_idx = (output_y * (*out_width) + output_x) * 4;++        output[output_idx] = 255;     // R+        output[output_idx + 1] = 255; // G+        output[output_idx + 2] = 255; // B+        output[output_idx + 3] = 255; // A+      }+    }+    else {+      // Draw overlapping colored histograms for RGB+      uint32_t bar_height_r = (hist_r[x] * hist_height) / max_count;+      uint32_t bar_height_g =+        (channels >= 3) ? (hist_g[x] * hist_height) / max_count : 0;+      uint32_t bar_height_b =+        (channels >= 3) ? (hist_b[x] * hist_height) / max_count : 0;++      uint32_t max_bar_height = bar_height_r;+      if (bar_height_g > max_bar_height) {+        max_bar_height = bar_height_g;+      }+      if (bar_height_b > max_bar_height) {+        max_bar_height = bar_height_b;+      }++      for (uint32_t y = 0; y < max_bar_height; y++) {+        uint32_t output_x = x;+        uint32_t output_y = hist_height - 1 - y;+        uint32_t output_idx = (output_y * (*out_width) + output_x) * 4;++        // Initialize to black+        uint8_t r = 0, g = 0, b = 0;++        // Add red channel contribution+        if (y < bar_height_r) {+          r = 255;+        }++        // Add green channel contribution+        if (channels >= 3 && y < bar_height_g) {+          g = 255;+        }++        // Add blue channel contribution+        if (channels >= 3 && y < bar_height_b) {+          b = 255;+        }++        output[output_idx] = r;       // R+        output[output_idx + 1] = g;   // G+        output[output_idx + 2] = b;   // B+        output[output_idx + 3] = 255; // A+      }+    }+  }++  return output;+}+// File: src/parse_hex_color.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "parse_hex_color.h"+#else+#include "flatcv.h"+#endif++// Helper function to parse hex color code to RGB values+void fcv_parse_hex_color(+  const char *hex_color,+  uint8_t *r,+  uint8_t *g,+  uint8_t *b+) {+  const char *hex = hex_color;++  // Skip '#' if present+  if (hex[0] == '#') {+    hex++;+  }++  // Validate hex string length (should be 6 characters)+  if (strlen(hex) != 6) {+    // Default to white for invalid hex codes+    *r = 255;+    *g = 255;+    *b = 255;+    return;+  }++  // Validate that all characters are valid hex digits+  for (int32_t i = 0; i < 6; i++) {+    if (!((hex[i] >= '0' && hex[i] <= '9') ||+          (hex[i] >= 'A' && hex[i] <= 'F') ||+          (hex[i] >= 'a' && hex[i] <= 'f'))) {+      // Default to white for invalid hex codes+      *r = 255;+      *g = 255;+      *b = 255;+      return;+    }+  }++  // Parse hex values+  uint32_t rgb_value;+  sscanf(hex, "%x", &rgb_value);++  *r = (rgb_value >> 16) & 0xFF;+  *g = (rgb_value >> 8) & 0xFF;+  *b = rgb_value & 0xFF;+}+// File: src/perspectivetransform.c+#include <assert.h>+#include <math.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "perspectivetransform.h"+#else+#include "flatcv.h"+#endif++// #define DEBUG_LOGGING++#ifdef DEBUG_LOGGING+#define log(msg) printf("DEBUG: %s\n", msg)+#else+#define log(msg) // No operation+#endif++/**+ * Helper function to solve 8x8 linear system using Gaussian elimination+ * Returns 1 on success, 0 on failure+ */+int32_t solve_linear_system(double A[8][8], double b[8], double x[8]) {+  const int32_t n = 8;+  const double epsilon = 1e-10;+  int32_t i, j, k;++  // Create augmented matrix [A|b] with extra safety margin+  double aug[8][10]; // One extra column for safety++  // Initialize augmented matrix+  for (i = 0; i < n; i++) {+    for (j = 0; j < n; j++) {+      aug[i][j] = A[i][j];+    }+    aug[i][n] = b[i];+  }++  // Gaussian elimination with partial pivoting+  for (i = 0; i < n; i++) {+    // Find pivot+    int32_t max_row = i;+    double max_val = fabs(aug[i][i]);++    for (k = i + 1; k < n; k++) {+      if (fabs(aug[k][i]) > max_val) {+        max_val = fabs(aug[k][i]);+        max_row = k;+      }+    }++    // Check for singularity+    if (max_val < epsilon) {+      log("Warning: Matrix is nearly singular\n");+      return 0;+    }++    // Swap maximum row with current row+    if (max_row != i) {+      for (j = 0; j <= n; j++) {+        double temp = aug[i][j];+        aug[i][j] = aug[max_row][j];+        aug[max_row][j] = temp;+      }+    }++    // Eliminate column i+    for (j = i + 1; j < n; j++) {+      double factor = aug[j][i] / aug[i][i];+      for (k = i; k <= n; k++) {+        aug[j][k] -= factor * aug[i][k];+      }+    }+  }++  // Back substitution+  for (i = n - 1; i >= 0; i--) {+    if (fabs(aug[i][i]) < epsilon) {+      log("Warning: Zero pivot encountered\n");+      return 0;+    }++    x[i] = aug[i][n];+    for (j = i + 1; j < n; j++) {+      x[i] -= aug[i][j] * x[j];+    }+    x[i] /= aug[i][i];++    // Check for invalid results+    if (isnan(x[i]) || isinf(x[i])) {+      log("Warning: Invalid result detected\n");+      return 0;+    }+  }++  return 1;+}++/**+ * Calculate the perspective transformation matrix+ * from the source and destination corner coordinates.+ */+Matrix3x3 *fcv_calculate_perspective_transform(+  Corners *src_corners,+  Corners *dst_corners+) {+  // Initialize matrices with zeros+  double A[8][8] = {{0}};+  double b[8] = {0};+  double x[8] = {0};++  // Identity matrix as fallback+  static Matrix3x3 identity = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};++  if (!src_corners || !dst_corners) {+    log("Error: NULL pointer passed to calculate_perspective_transform\n");+    return &identity;+  }++#ifdef DEBUG_LOGGING+  fprintf(stderr, "[C] Calculating perspective transform:\n");+  fprintf(+    stderr,+    "src_corners:\ntl(%f, %f)\ntr(%f, %f)\nbr(%f, %f)\nbl(%f, %f)\n\n",+    src_corners->tl_x,+    src_corners->tl_y,+    src_corners->tr_x,+    src_corners->tr_y,+    src_corners->br_x,+    src_corners->br_y,+    src_corners->bl_x,+    src_corners->bl_y+  );+  fprintf(+    stderr,+    "dst_corners:\ntl(%f, %f)\ntr(%f, %f)\nbr(%f, %f)\nbl(%f, %f)\n\n",+    dst_corners->tl_x,+    dst_corners->tl_y,+    dst_corners->tr_x,+    dst_corners->tr_y,+    dst_corners->br_x,+    dst_corners->br_y,+    dst_corners->bl_x,+    dst_corners->bl_y+  );+#endif++  // Validate input coordinates+  if (isnan(src_corners->tl_x) || isnan(src_corners->tl_y) ||+      isnan(src_corners->tr_x) || isnan(src_corners->tr_y) ||+      isnan(src_corners->br_x) || isnan(src_corners->br_y) ||+      isnan(src_corners->bl_x) || isnan(src_corners->bl_y) ||+      isnan(dst_corners->tl_x) || isnan(dst_corners->tl_y) ||+      isnan(dst_corners->tr_x) || isnan(dst_corners->tr_y) ||+      isnan(dst_corners->br_x) || isnan(dst_corners->br_y) ||+      isnan(dst_corners->bl_x) || isnan(dst_corners->bl_y)) {+    log("Error: Invalid coordinates (NaN) detected\n");+    return &identity;+  }++  // Set up the system of equations+  for (int32_t i = 0; i < 4; i++) {+    double srcX = 0.0, srcY = 0.0, dstX = 0.0, dstY = 0.0;++    // Safely extract coordinates+    switch (i) {+    case 0: // Top-left+      srcX = src_corners->tl_x;+      srcY = src_corners->tl_y;+      dstX = dst_corners->tl_x;+      dstY = dst_corners->tl_y;+      break;+    case 1: // Top-right+      srcX = src_corners->tr_x;+      srcY = src_corners->tr_y;+      dstX = dst_corners->tr_x;+      dstY = dst_corners->tr_y;+      break;+    case 2: // Bottom-right+      srcX = src_corners->br_x;+      srcY = src_corners->br_y;+      dstX = dst_corners->br_x;+      dstY = dst_corners->br_y;+      break;+    case 3: // Bottom-left+      srcX = src_corners->bl_x;+      srcY = src_corners->bl_y;+      dstX = dst_corners->bl_x;+      dstY = dst_corners->bl_y;+      break;+    }++    // Validate extracted coordinates+    if (isinf(srcX) || isinf(srcY) || isinf(dstX) || isinf(dstY)) {+      log("Error: Invalid coordinates (Inf) detected\n");+      return &identity;+    }++    // First four equations for x coordinates+    A[i][0] = srcX;+    A[i][1] = srcY;+    A[i][2] = 1.0;+    A[i][6] = -srcX * dstX;+    A[i][7] = -srcY * dstX;+    b[i] = dstX;++    // Last four equations for y coordinates+    A[i + 4][3] = srcX;+    A[i + 4][4] = srcY;+    A[i + 4][5] = 1.0;+    A[i + 4][6] = -srcX * dstY;+    A[i + 4][7] = -srcY * dstY;+    b[i + 4] = dstY;+  }++  log("Solve the system of equations …\n");+  if (!solve_linear_system(A, b, x)) {+    log("Failed to solve system, returning identity matrix\n");+    return &identity;+  }++  // Validate solution+  for (int32_t i = 0; i < 8; i++) {+    if (isnan(x[i]) || isinf(x[i]) || fabs(x[i]) > 1e6) {+      log("Error: Invalid solution values detected\n");+      return &identity;+    }+  }++  Matrix3x3 *result = malloc(sizeof(Matrix3x3));+  *result = (Matrix3x3){x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], 1.0};++#ifdef DEBUG_LOGGING+  fprintf(stderr, "Result matrix:\n");+  fprintf(stderr, "%f, %f, %f\n", result->m00, result->m01, result->m02);+  fprintf(stderr, "%f, %f, %f\n", result->m10, result->m11, result->m12);+  fprintf(stderr, "%f, %f, %f\n", result->m20, result->m21, result->m22);+#endif++  // Final validation of the result matrix+  if (isnan(result->m00) || isnan(result->m01) || isnan(result->m02) ||+      isnan(result->m10) || isnan(result->m11) || isnan(result->m12) ||+      isnan(result->m20) || isnan(result->m21) || isnan(result->m22) ||+      isinf(result->m00) || isinf(result->m01) || isinf(result->m02) ||+      isinf(result->m10) || isinf(result->m11) || isinf(result->m12) ||+      isinf(result->m20) || isinf(result->m21) || isinf(result->m22)) {+    log("Error: Invalid values in result matrix\n");+    return &identity;+  }++  return result;+}++/**+ * Apply the transformation matrix to the input image+ * and store the result in the output image.+ * Use bilinear interpolation to calculate final pixel values.+ */+uint8_t *fcv_apply_matrix_3x3(+  int32_t in_width,+  int32_t in_height,+  uint8_t *in_data,+  int32_t out_width,+  int32_t out_height,+  Matrix3x3 *tmat+) {+  if (!in_data || !tmat) {+    return NULL;+  }++  if (in_width <= 0 || in_height <= 0 || out_width <= 0 || out_height <= 0) {+    return NULL;+  }++  // Check for overflow: out_width * out_height * 4 (dimensions already+  // validated > 0)+  if ((size_t)out_width > SIZE_MAX / (size_t)out_height) {+    return NULL;+  }+  size_t out_pixels = (size_t)out_width * (size_t)out_height;+  if (out_pixels > SIZE_MAX / 4) {+    return NULL;+  }++  // Patch flip matrix if needed+  if (fabs(tmat->m00 + 1.0) < 1e-9 && fabs(tmat->m11 + 1.0) < 1e-9 &&+      tmat->m02 == 0.0 && tmat->m12 == 0.0) {+    tmat->m02 = in_width - 1;+    tmat->m12 = in_height - 1;+  }++  uint8_t *out_data = calloc(out_pixels * 4, sizeof(uint8_t));++  if (!out_data) { // Memory allocation failed+    return NULL;+  }++  // Iterate through every pixel in the output image+  for (int32_t out_y = 0; out_y < out_height; ++out_y) {+    for (int32_t out_x = 0; out_x < out_width; ++out_x) {+      // Apply the inverse transformation to find the corresponding source pixel+      double w = tmat->m20 * out_x + tmat->m21 * out_y + tmat->m22;+      if (fabs(w) < 1e-10) {+        continue; // Skip if w is too close to zero+      }++      double srcX = (tmat->m00 * out_x + tmat->m01 * out_y + tmat->m02) / w;+      double srcY = (tmat->m10 * out_x + tmat->m11 * out_y + tmat->m12) / w;++      // Convert source coordinates to integers+      int32_t x0 = (int32_t)floor(srcX);+      int32_t y0 = (int32_t)floor(srcY);+      int32_t x1 = x0 + 1;+      int32_t y1 = y0 + 1;++      // Verify that the anchor pixel is inside the source image+      if (x0 >= 0 && x0 < in_width && y0 >= 0 && y0 < in_height) {++        // Clamp the neighbor coordinates so that a (degenerated)+        // bilinear interpolation can be applied at the image borders.+        int32_t x1c = (x1 < in_width) ? x1 : x0;+        int32_t y1c = (y1 < in_height) ? y1 : y0;++        double dx = srcX - x0;+        double dy = srcY - y0;++        // If a neighbour got clamped we force the corresponding weight to 0+        if (x1c == x0) {+          dx = 0.0;+        }+        if (y1c == y0) {+          dy = 0.0;+        }++        uint8_t *p00 = &in_data[(y0 * in_width + x0) * 4];+        uint8_t *p01 = &in_data[(y0 * in_width + x1c) * 4];+        uint8_t *p10 = &in_data[(y1c * in_width + x0) * 4];+        uint8_t *p11 = &in_data[(y1c * in_width + x1c) * 4];++        for (int32_t c = 0; c < 4; ++c) {+          out_data[(out_y * out_width + out_x) * 4 + c] =+            (uint8_t)(p00[c] * (1 - dx) * (1 - dy) + p01[c] * dx * (1 - dy) ++                      p10[c] * (1 - dx) * dy + p11[c] * dx * dy);+        }+      }+    }+  }++  return out_data;+}+// File: src/rgba_to_grayscale.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "rgba_to_grayscale.h"+#else+#include "flatcv.h"+#endif++/**+ * Convert raw RGBA row-major top-to-bottom image data+ * to a single channel grayscale image data.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the single channel grayscale image data.+ */+uint8_t *fcv_rgba_to_grayscale(+  uint32_t width,+  uint32_t height,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;++  uint8_t *grayscale_data = malloc(img_length_px);++  if (!grayscale_data) { // Memory allocation failed+    return NULL;+  }++  // Process each pixel row by row+  for (size_t i = 0; i < img_length_px; i++) {+    size_t rgba_index = i * 4;++    uint8_t r = data[rgba_index];+    uint8_t g = data[rgba_index + 1];+    uint8_t b = data[rgba_index + 2];++    uint8_t gray = (r * R_WEIGHT + g * G_WEIGHT + b * B_WEIGHT) >> 8;++    grayscale_data[i] = gray;+  }++  return grayscale_data;+}+// File: src/rotate.c+#include <stdint.h>+#include <stdlib.h>+#include <string.h>++#ifndef FLATCV_AMALGAMATION+#include "rotate.h"+#else+#include "flatcv.h"+#endif++/**+ * Rotate an image 90 degrees clockwise.+ */+uint8_t *+fcv_rotate_90_cw(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *rotated_data = malloc(img_length_byte);++  if (!rotated_data) {+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index = ((size_t)x * height + (height - 1 - y)) * 4;++      rotated_data[dst_index] = data[src_index];+      rotated_data[dst_index + 1] = data[src_index + 1];+      rotated_data[dst_index + 2] = data[src_index + 2];+      rotated_data[dst_index + 3] = data[src_index + 3];+    }+  }++  return rotated_data;+}++/**+ * Rotate an image 180 degrees.+ */+uint8_t *+fcv_rotate_180(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *rotated_data = malloc(img_length_byte);++  if (!rotated_data) {+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index =+        ((size_t)(height - 1 - y) * width + (width - 1 - x)) * 4;++      rotated_data[dst_index] = data[src_index];+      rotated_data[dst_index + 1] = data[src_index + 1];+      rotated_data[dst_index + 2] = data[src_index + 2];+      rotated_data[dst_index + 3] = data[src_index + 3];+    }+  }++  return rotated_data;+}++/**+ * Rotate an image 270 degrees clockwise (90 degrees counter-clockwise).+ */+uint8_t *+fcv_rotate_270_cw(uint32_t width, uint32_t height, uint8_t const *const data) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t num_pixels = (size_t)width * height;+  if (num_pixels > SIZE_MAX / 4) {+    return NULL;+  }+  size_t img_length_byte = num_pixels * 4;++  uint8_t *rotated_data = malloc(img_length_byte);++  if (!rotated_data) {+    return NULL;+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      size_t src_index = ((size_t)y * width + x) * 4;+      size_t dst_index = ((size_t)(width - 1 - x) * height + y) * 4;++      rotated_data[dst_index] = data[src_index];+      rotated_data[dst_index + 1] = data[src_index + 1];+      rotated_data[dst_index + 2] = data[src_index + 2];+      rotated_data[dst_index + 3] = data[src_index + 3];+    }+  }++  return rotated_data;+}+// File: src/single_to_multichannel.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "single_to_multichannel.h"+#else+#include "flatcv.h"+#endif++/**+ * Convert single channel grayscale image data to+ * RGBA row-major top-to-bottom image data.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param data Pointer to the pixel data.+ */+uint8_t *fcv_single_to_multichannel(+  uint32_t width,+  uint32_t height,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0) {+    return NULL;+  }++  // Check for overflow: width * height * 4+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;+  if (img_length_px > SIZE_MAX / 4) {+    return NULL;+  }++  uint8_t *multichannel_data = malloc(img_length_px * 4);++  if (!multichannel_data) { // Memory allocation failed+    return NULL;+  }++  for (size_t i = 0; i < img_length_px; i++) {+    size_t rgba_index = i * 4;+    multichannel_data[rgba_index] = data[i];+    multichannel_data[rgba_index + 1] = data[i];+    multichannel_data[rgba_index + 2] = data[i];+    multichannel_data[rgba_index + 3] = 255;+  }++  return multichannel_data;+}+// File: src/sobel_edge_detection.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "conversion.h"+#include "perspectivetransform.h"+#include "rgba_to_grayscale.h"+#include "sobel_edge_detection.h"+#else+#include "flatcv.h"+#endif++/**+ * Apply Sobel edge detection to the image data and return single-channel+ * grayscale data. Uses Sobel kernels to detect edges in horizontal and vertical+ * directions, then combines them to get the edge magnitude.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param channels Number of channels in the input image.+ *                 (1: Grayscale, 3: RGB, 4: RGBA)+ * @param data Pointer to the input pixel data.+ * @return Pointer to the single-channel grayscale edge-detected image data.+ */+uint8_t *fcv_sobel_edge_detection(+  uint32_t width,+  uint32_t height,+  uint32_t channels,+  uint8_t const *const data+) {+  if (!data || width == 0 || height == 0 || channels == 0) {+    return NULL;+  }++  // Check for overflow: width * height+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;++  // Check for overflow in magnitudes allocation+  if (img_length_px > SIZE_MAX / sizeof(double)) {+    return NULL;+  }++  uint8_t *grayscale_data;+  bool allocated_grayscale = false;++  if (channels == 1) {+    // Single-channel input, use data directly+    grayscale_data = (uint8_t *)data;+    allocated_grayscale = false;+  }+  else {+    // Multi-channel input, convert to grayscale+    grayscale_data = fcv_rgba_to_grayscale(width, height, data);+    if (!grayscale_data) {+      return NULL;+    }+    allocated_grayscale = true;+  }++  uint8_t *sobel_data = malloc(img_length_px);+  if (!sobel_data) {+    if (allocated_grayscale) {+      free(grayscale_data);+    }+    return NULL;+  }++  // Sobel kernels+  int32_t sobel_x[3][3] = {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}};+  int32_t sobel_y[3][3] = {{-1, -2, -1}, {0, 0, 0}, {1, 2, 1}};++  // Temporary buffer to store magnitudes for normalization+  double *magnitudes = malloc(img_length_px * sizeof(double));+  if (!magnitudes) {+    if (allocated_grayscale) {+      free(grayscale_data);+    }+    free(sobel_data);+    return NULL;+  }++  double min_magnitude = INFINITY;+  double max_magnitude = 0.0;++  // First pass: calculate all magnitudes and find min/max+  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      int32_t gx = 0;+      int32_t gy = 0;++      // Apply Sobel kernels+      for (int32_t ky = -1; ky <= 1; ky++) {+        for (int32_t kx = -1; kx <= 1; kx++) {+          int32_t px = x + kx;+          int32_t py = y + ky;++          // Handle boundaries by using nearest pixel values+          if (px < 0) {+            px = 0;+          }+          if ((uint32_t)px >= width) {+            px = width - 1;+          }+          if (py < 0) {+            py = 0;+          }+          if ((uint32_t)py >= height) {+            py = height - 1;+          }++          uint8_t pixel = grayscale_data[(size_t)py * width + px];+          gx += pixel * sobel_x[ky + 1][kx + 1];+          gy += pixel * sobel_y[ky + 1][kx + 1];+        }+      }++      double magnitude = sqrt((double)gx * gx + (double)gy * gy);+      magnitudes[(size_t)y * width + x] = magnitude;++      if (magnitude < min_magnitude) {+        min_magnitude = magnitude;+      }+      if (magnitude > max_magnitude) {+        max_magnitude = magnitude;+      }+    }+  }++  // Second pass: normalize and assign to output+  double range = max_magnitude - min_magnitude;+  if (range == 0.0) {+    range = 1.0; // Avoid division by zero+  }++  for (uint32_t y = 0; y < height; y++) {+    for (uint32_t x = 0; x < width; x++) {+      double magnitude = magnitudes[(size_t)y * width + x];++      // Normalize to 0-255 range based on actual min/max+      int32_t final_magnitude =+        (int32_t)(((magnitude - min_magnitude) / range) * 255.0);++      // Clamp to valid range+      if (final_magnitude < 0) {+        final_magnitude = 0;+      }+      if (final_magnitude > 255) {+        final_magnitude = 255;+      }++      sobel_data[(size_t)y * width + x] = (uint8_t)final_magnitude;+    }+  }++  free(magnitudes);++  if (allocated_grayscale) {+    free(grayscale_data);+  }+  return sobel_data;+}+// File: src/sort_corners.c+#include <math.h>+#include <stdbool.h>+#include <stdio.h>+#include <stdlib.h>++#ifndef FLATCV_AMALGAMATION+#include "1_types.h"+#else+#include "flatcv.h"+#endif++typedef struct {+  double r;     // radius+  double theta; // angle+} PolarPoint;++// Helper function to convert Cartesian to polar coordinates+PolarPoint cartesian_to_polar(Point2D cart) {+  PolarPoint polar;+  polar.r = sqrt(cart.x * cart.x + cart.y * cart.y);+  double theta = atan2(cart.y, cart.x);++  // Normalize theta: Convert to degrees and apply normalization+  double theta_deg = theta * 180.0 / M_PI;+  if (theta_deg < 0) {+    polar.theta = -theta_deg;+  }+  else {+    polar.theta = -(theta_deg - 360.0);+  }++  return polar;+}++// Comparison function for qsort (sorting by angle in descending order)+int compare_polar_angles(const void *a, const void *b) {+  const PolarPoint *pa = (const PolarPoint *)a;+  const PolarPoint *pb = (const PolarPoint *)b;++  if (pa->theta < pb->theta) {+    return 1; // Descending order+  }+  if (pa->theta > pb->theta) {+    return -1;+  }+  return 0;+}++// Helper structure to keep track of original indices during sorting+typedef struct {+  PolarPoint polar;+  uint32_t original_index;+} IndexedPolarPoint;++int compare_indexed_polar_angles(const void *a, const void *b) {+  const IndexedPolarPoint *pa = (const IndexedPolarPoint *)a;+  const IndexedPolarPoint *pb = (const IndexedPolarPoint *)b;++  if (pa->polar.theta < pb->polar.theta) {+    return 1; // Descending order+  }+  if (pa->polar.theta > pb->polar.theta) {+    return -1;+  }+  return 0;+}++Corners sort_corners(+  uint32_t width,+  uint32_t height,+  uint32_t out_width,+  uint32_t out_height,+  Point2D *corners,+  uint32_t num_corners,+  Point2D *result+) {+  if (num_corners < 3) {+    // Not enough corners, return zeros+    Corners empty_corners = {0, 0, 0, 0, 0, 0, 0, 0};+    return empty_corners;+  }++  // Allocate memory for working corners (max of num_corners or 4 for 3-corner+  // case)+  uint32_t max_corners = (num_corners >= 4) ? num_corners : 4;+  Point2D *working_corners = malloc(max_corners * sizeof(Point2D));+  if (!working_corners) {+    Corners empty_corners = {0, 0, 0, 0, 0, 0, 0, 0};+    return empty_corners;+  }++  uint32_t corners_to_process;++  if (num_corners == 3) {+    // Copy the 3 corners we have+    working_corners[0] = corners[0];+    working_corners[1] = corners[1];+    working_corners[2] = corners[2];++    // Calculate the fourth corner to complete a parallelogram+    // For 3 corners A, B, C, we need to determine which corner is missing+    // and calculate it using vector addition++    // First, sort the 3 corners by position to understand the layout+    Point2D sorted[3] = {corners[0], corners[1], corners[2]};++    // Find bounding box of the 3 corners+    double min_x = sorted[0].x, max_x = sorted[0].x;+    double min_y = sorted[0].y, max_y = sorted[0].y;+    for (int i = 1; i < 3; i++) {+      if (sorted[i].x < min_x) {+        min_x = sorted[i].x;+      }+      if (sorted[i].x > max_x) {+        max_x = sorted[i].x;+      }+      if (sorted[i].y < min_y) {+        min_y = sorted[i].y;+      }+      if (sorted[i].y > max_y) {+        max_y = sorted[i].y;+      }+    }++    // Determine which corner is missing by finding which position in the+    // rectangle is empty+    bool has_tl = false, has_tr = false, has_bl = false, has_br = false;+    for (int i = 0; i < 3; i++) {+      double x = sorted[i].x;+      double y = sorted[i].y;++      // Use small tolerance for floating point comparison+      double tolerance = 1.0;++      if (fabs(x - min_x) < tolerance && fabs(y - min_y) < tolerance) {+        has_tl = true;+      }+      else if (fabs(x - max_x) < tolerance && fabs(y - min_y) < tolerance) {+        has_tr = true;+      }+      else if (fabs(x - min_x) < tolerance && fabs(y - max_y) < tolerance) {+        has_bl = true;+      }+      else if (fabs(x - max_x) < tolerance && fabs(y - max_y) < tolerance) {+        has_br = true;+      }+    }++    // Calculate the missing corner+    Point2D missing_corner;+    if (!has_tl) {+      missing_corner.x = min_x;+      missing_corner.y = min_y;+    }+    else if (!has_tr) {+      missing_corner.x = max_x;+      missing_corner.y = min_y;+    }+    else if (!has_bl) {+      missing_corner.x = min_x;+      missing_corner.y = max_y;+    }+    else if (!has_br) {+      missing_corner.x = max_x;+      missing_corner.y = max_y;+    }+    else {+      // Fallback: use parallelogram completion+      // For a parallelogram ABCD where we have A, B, C and need D:+      // D = A + C - B (assuming B and D are diagonal)+      // Try all combinations and pick the best one+      Point2D candidates[3];+      candidates[0].x = corners[0].x + corners[1].x - corners[2].x;+      candidates[0].y = corners[0].y + corners[1].y - corners[2].y;+      candidates[1].x = corners[0].x + corners[2].x - corners[1].x;+      candidates[1].y = corners[0].y + corners[2].y - corners[1].y;+      candidates[2].x = corners[1].x + corners[2].x - corners[0].x;+      candidates[2].y = corners[1].y + corners[2].y - corners[0].y;++      // Choose the candidate that gives the most reasonable result+      missing_corner = candidates[2]; // Default to the last one+      double best_score = -1;++      for (int i = 0; i < 3; i++) {+        double score = 0;+        if (candidates[i].x >= 0) {+          score += 1;+        }+        if (candidates[i].y >= 0) {+          score += 1;+        }++        // Check if it's not a duplicate+        bool is_duplicate = false;+        for (int j = 0; j < 3; j++) {+          if (fabs(candidates[i].x - corners[j].x) < 0.1 &&+              fabs(candidates[i].y - corners[j].y) < 0.1) {+            is_duplicate = true;+            break;+          }+        }+        if (!is_duplicate) {+          score += 5;+        }++        if (score > best_score) {+          best_score = score;+          missing_corner = candidates[i];+        }+      }+    }++    working_corners[3] = missing_corner;+    corners_to_process = 4;+  }+  else {+    // Copy all existing corners+    for (uint32_t i = 0; i < num_corners; i++) {+      working_corners[i] = corners[i];+    }+    corners_to_process = num_corners;+  }++  // Clockwise corner sorting starting from top-left corner+  // Works with any number of corners >= 3++  // Step 1: Find the center point (centroid)+  double center_x = 0.0, center_y = 0.0;+  for (uint32_t i = 0; i < corners_to_process; i++) {+    center_x += working_corners[i].x;+    center_y += working_corners[i].y;+  }+  center_x /= corners_to_process;+  center_y /= corners_to_process;++  // Step 2: Find the top-left corner (minimum x+y sum)+  uint32_t top_left_idx = 0;+  double min_sum = working_corners[0].x + working_corners[0].y;+  for (uint32_t i = 1; i < corners_to_process; i++) {+    double sum = working_corners[i].x + working_corners[i].y;+    if (sum < min_sum) {+      min_sum = sum;+      top_left_idx = i;+    }+  }++  // Step 3: Create array with indices and calculate angles relative to center+  typedef struct {+    uint32_t index;+    double angle;+    double distance;+  } CornerInfo;++  CornerInfo *corner_info = malloc(corners_to_process * sizeof(CornerInfo));+  if (!corner_info) {+    free(working_corners);+    Corners empty_corners = {0, 0, 0, 0, 0, 0, 0, 0};+    return empty_corners;+  }++  for (uint32_t i = 0; i < corners_to_process; i++) {+    corner_info[i].index = i;+    corner_info[i].distance = sqrt(+      (working_corners[i].x - center_x) * (working_corners[i].x - center_x) ++      (working_corners[i].y - center_y) * (working_corners[i].y - center_y)+    );++    // Calculate angle from center to corner+    double dx = working_corners[i].x - center_x;+    double dy = working_corners[i].y - center_y;+    corner_info[i].angle = atan2(dy, dx);++    // Normalize angle to [0, 2π) range+    if (corner_info[i].angle < 0) {+      corner_info[i].angle += 2.0 * M_PI;+    }+  }++  // Step 4: Calculate reference angle from center to top-left corner+  double top_left_angle = corner_info[top_left_idx].angle;++  // Step 5: Adjust all angles relative to top-left corner and sort clockwise+  for (uint32_t i = 0; i < corners_to_process; i++) {+    // Adjust angle relative to top-left corner+    corner_info[i].angle = corner_info[i].angle - top_left_angle;++    // Normalize to [0, 2π) range+    while (corner_info[i].angle < 0) {+      corner_info[i].angle += 2.0 * M_PI;+    }+    while (corner_info[i].angle >= 2.0 * M_PI) {+      corner_info[i].angle -= 2.0 * M_PI;+    }+  }++  // Step 6: Sort corners by angle (clockwise from top-left)+  for (uint32_t i = 0; i < corners_to_process - 1; i++) {+    for (uint32_t j = 0; j < corners_to_process - 1 - i; j++) {+      bool should_swap = false;++      if (corner_info[j].angle > corner_info[j + 1].angle) {+        should_swap = true;+      }+      else if (fabs(corner_info[j].angle - corner_info[j + 1].angle) < 1e-6) {+        // If angles are very close, sort by distance from center (closer first)+        if (corner_info[j].distance > corner_info[j + 1].distance) {+          should_swap = true;+        }+      }++      if (should_swap) {+        CornerInfo temp = corner_info[j];+        corner_info[j] = corner_info[j + 1];+        corner_info[j + 1] = temp;+      }+    }+  }++  // Step 7: Copy sorted corners to result+  for (uint32_t i = 0; i < corners_to_process; i++) {+    result[i] = working_corners[corner_info[i].index];+  }++  // Scale corners back to original image dimensions+  double scale_x = (double)width / out_width;+  double scale_y = (double)height / out_height;++  // Return corners struct using the sorted result+  Corners sorted_corners_result = {+    .tl_x = result[0].x * scale_x,+    .tl_y = result[0].y * scale_y,+    .tr_x = result[1].x * scale_x,+    .tr_y = result[1].y * scale_y,+    .br_x = result[2].x * scale_x,+    .br_y = result[2].y * scale_y,+    .bl_x = result[3].x * scale_x,+    .bl_y = result[3].y * scale_y+  };++  // Clean up allocated memory+  free(corner_info);+  free(working_corners);++  return sorted_corners_result;+}+// File: src/trim.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "crop.h"+#include "trim.h"+#else+#include "flatcv.h"+#endif++static bool+pixels_match(const uint8_t *pixel1, const uint8_t *pixel2, uint32_t channels) {+  for (uint32_t c = 0; c < channels; c++) {+    if (pixel1[c] != pixel2[c]) {+      return false;+    }+  }+  return true;+}++/**+ * Check if two pixels match within a given tolerance.+ *+ * @param pixel1 First pixel data.+ * @param pixel2 Second pixel data.+ * @param channels Number of channels.+ * @param tolerance Maximum allowed difference per channel (0-255).+ * @return true if all channels are within tolerance, false otherwise.+ */+static bool pixels_match_threshold(+  const uint8_t *pixel1,+  const uint8_t *pixel2,+  uint32_t channels,+  uint8_t tolerance+) {+  for (uint32_t c = 0; c < channels; c++) {+    int32_t diff = (int32_t)pixel1[c] - (int32_t)pixel2[c];+    if (diff < 0) {+      diff = -diff;+    }+    if ((uint8_t)diff > tolerance) {+      return false;+    }+  }+  return true;+}++/**+ * Trim border pixels that have the same color.+ *+ * @param width Pointer to width of the image (will be updated).+ * @param height Pointer to height of the image (will be updated).+ * @param channels Number of channels in the image.+ * @param data Pointer to the pixel data.+ * @return Pointer to the new trimmed image data.+ */+uint8_t *fcv_trim(+  int32_t *width,+  int32_t *height,+  uint32_t channels,+  uint8_t const *const data+) {+  if (!data || !width || !height || *width <= 0 || *height <= 0 ||+      channels == 0) {+    return NULL;+  }++  uint32_t w = (uint32_t)*width;+  uint32_t h = (uint32_t)*height;++  // Check for potential overflow in pixel index calculations+  if (w > SIZE_MAX / h || (size_t)w * h > SIZE_MAX / channels) {+    return NULL;+  }+  uint32_t left = 0;+  uint32_t top = 0;+  uint32_t right = w;+  uint32_t bottom = h;++  // Trim from the left+  while (left < right) {+    const uint8_t *ref_pixel = &data[(top * w + left) * channels];+    bool can_trim = true;++    // Check if the entire left column has the same color as the reference pixel+    for (uint32_t y = top; y < bottom; y++) {+      const uint8_t *pixel = &data[(y * w + left) * channels];+      if (!pixels_match(pixel, ref_pixel, channels)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim && right - left > 1) {+      left++;+    }+    else {+      break;+    }+  }++  // Trim from the right+  while (left < right) {+    const uint8_t *ref_pixel = &data[(top * w + (right - 1)) * channels];+    bool can_trim = true;++    // Check if the entire right column has the same color as the reference+    // pixel+    for (uint32_t y = top; y < bottom; y++) {+      const uint8_t *pixel = &data[(y * w + (right - 1)) * channels];+      if (!pixels_match(pixel, ref_pixel, channels)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim && right - left > 1) {+      right--;+    }+    else {+      break;+    }+  }++  // Trim from the top+  while (top < bottom && bottom - top > 1) {+    const uint8_t *ref_pixel = &data[(top * w + left) * channels];+    bool can_trim = true;++    // Check if the entire top row has the same color as the reference pixel+    for (uint32_t x = left; x < right; x++) {+      const uint8_t *pixel = &data[(top * w + x) * channels];+      if (!pixels_match(pixel, ref_pixel, channels)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim) {+      top++;+    }+    else {+      break;+    }+  }++  // Trim from the bottom+  while (top < bottom && bottom - top > 1) {+    const uint8_t *ref_pixel = &data[((bottom - 1) * w + left) * channels];+    bool can_trim = true;++    // Check if the entire bottom row has the same color as the reference pixel+    for (uint32_t x = left; x < right; x++) {+      const uint8_t *pixel = &data[((bottom - 1) * w + x) * channels];+      if (!pixels_match(pixel, ref_pixel, channels)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim) {+      bottom--;+    }+    else {+      break;+    }+  }++  // Calculate new dimensions+  uint32_t new_width = right - left;+  uint32_t new_height = bottom - top;++  // If no trimming was done, return a copy+  if (left == 0 && top == 0 && right == w && bottom == h) {+    size_t alloc_size = (size_t)w * h * channels;+    uint8_t *result = malloc(alloc_size);+    if (result) {+      memcpy(result, data, alloc_size);+    }+    return result;+  }++  // Use existing crop function to extract the trimmed region+  uint8_t *trimmed_data =+    fcv_crop(w, h, channels, data, left, top, new_width, new_height);++  // Update dimensions+  *width = (int32_t)new_width;+  *height = (int32_t)new_height;++  return trimmed_data;+}++/**+ * Trim border pixels that are within a threshold of the reference color.+ * This is useful for images with JPEG artifacts or slight vignetting.+ *+ * @param width Pointer to width of the image (will be updated).+ * @param height Pointer to height of the image (will be updated).+ * @param channels Number of channels in the image.+ * @param data Pointer to the pixel data.+ * @param threshold_percent Tolerance percentage (0-100). A value of 2 means+ *        pixels within 2% (≈5 units) of the reference color will be trimmed.+ * @return Pointer to the new trimmed image data.+ */+uint8_t *fcv_trim_threshold(+  int32_t *width,+  int32_t *height,+  uint32_t channels,+  uint8_t const *const data,+  double threshold_percent+) {+  if (!data || !width || !height || *width <= 0 || *height <= 0 ||+      channels == 0) {+    return NULL;+  }++  // Validate threshold_percent - reject invalid values+  if (!isfinite(threshold_percent)) {+    return NULL;+  }++  // Clamp threshold to valid range and convert to absolute value (0-255)+  if (threshold_percent < 0.0) {+    threshold_percent = 0.0;+  }+  if (threshold_percent > 100.0) {+    threshold_percent = 100.0;+  }+  uint8_t tolerance = (uint8_t)(threshold_percent * 255.0 / 100.0);++  // If tolerance is 0, use the exact match version+  if (tolerance == 0) {+    return fcv_trim(width, height, channels, data);+  }++  uint32_t w = (uint32_t)*width;+  uint32_t h = (uint32_t)*height;++  // Check for potential overflow in pixel index calculations+  if (w > SIZE_MAX / h || (size_t)w * h > SIZE_MAX / channels) {+    return NULL;+  }+  uint32_t left = 0;+  uint32_t top = 0;+  uint32_t right = w;+  uint32_t bottom = h;++  // Trim from the left+  while (left < right) {+    const uint8_t *ref_pixel = &data[(top * w + left) * channels];+    bool can_trim = true;++    // Check if the entire left column matches within threshold+    for (uint32_t y = top; y < bottom; y++) {+      const uint8_t *pixel = &data[(y * w + left) * channels];+      if (!pixels_match_threshold(pixel, ref_pixel, channels, tolerance)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim && right - left > 1) {+      left++;+    }+    else {+      break;+    }+  }++  // Trim from the right+  while (left < right) {+    const uint8_t *ref_pixel = &data[(top * w + (right - 1)) * channels];+    bool can_trim = true;++    // Check if the entire right column matches within threshold+    for (uint32_t y = top; y < bottom; y++) {+      const uint8_t *pixel = &data[(y * w + (right - 1)) * channels];+      if (!pixels_match_threshold(pixel, ref_pixel, channels, tolerance)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim && right - left > 1) {+      right--;+    }+    else {+      break;+    }+  }++  // Trim from the top+  while (top < bottom && bottom - top > 1) {+    const uint8_t *ref_pixel = &data[(top * w + left) * channels];+    bool can_trim = true;++    // Check if the entire top row matches within threshold+    for (uint32_t x = left; x < right; x++) {+      const uint8_t *pixel = &data[(top * w + x) * channels];+      if (!pixels_match_threshold(pixel, ref_pixel, channels, tolerance)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim) {+      top++;+    }+    else {+      break;+    }+  }++  // Trim from the bottom+  while (top < bottom && bottom - top > 1) {+    const uint8_t *ref_pixel = &data[((bottom - 1) * w + left) * channels];+    bool can_trim = true;++    // Check if the entire bottom row matches within threshold+    for (uint32_t x = left; x < right; x++) {+      const uint8_t *pixel = &data[((bottom - 1) * w + x) * channels];+      if (!pixels_match_threshold(pixel, ref_pixel, channels, tolerance)) {+        can_trim = false;+        break;+      }+    }++    if (can_trim) {+      bottom--;+    }+    else {+      break;+    }+  }++  // Calculate new dimensions+  uint32_t new_width = right - left;+  uint32_t new_height = bottom - top;++  // If no trimming was done, return a copy+  if (left == 0 && top == 0 && right == w && bottom == h) {+    size_t alloc_size = (size_t)w * h * channels;+    uint8_t *result = malloc(alloc_size);+    if (result) {+      memcpy(result, data, alloc_size);+    }+    return result;+  }++  // Use existing crop function to extract the trimmed region+  uint8_t *trimmed_data =+    fcv_crop(w, h, channels, data, left, top, new_width, new_height);++  // Update dimensions+  *width = (int32_t)new_width;+  *height = (int32_t)new_height;++  return trimmed_data;+}+// File: src/watershed_segmentation.c+#include <assert.h>+#include <math.h>+#include <stdbool.h>+#include <stdint.h>+#include <stdio.h>+#include <stdlib.h>+#include <string.h>+#include <time.h>++#ifndef FLATCV_AMALGAMATION+#include "1_types.h"+#include "watershed_segmentation.h"+#else+#include "flatcv.h"+#endif++typedef struct {+  int32_t x, y;+  int32_t label;+} QueueItem;++typedef struct {+  QueueItem *items;+  int32_t front;+  int32_t rear;+  int32_t size;+  int32_t capacity;+} Queue;++static Queue *create_queue(size_t capacity) {+  if (capacity == 0 || capacity > SIZE_MAX / sizeof(QueueItem) ||+      capacity > INT32_MAX) {+    return NULL;+  }++  Queue *q = malloc(sizeof(Queue));+  if (!q) {+    return NULL;+  }++  q->items = malloc(capacity * sizeof(QueueItem));+  if (!q->items) {+    free(q);+    return NULL;+  }++  q->front = 0;+  q->rear = -1;+  q->size = 0;+  q->capacity = (int32_t)capacity;+  return q;+}++static void destroy_queue(Queue *q) {+  if (q) {+    free(q->items);+    free(q);+  }+}++static int32_t is_queue_empty(Queue *q) { return q->size == 0; }++static int32_t enqueue(Queue *q, int32_t x, int32_t y, int32_t label) {+  if (q->size >= q->capacity) {+    return 0;+  }++  q->rear = (q->rear + 1) % q->capacity;+  q->items[q->rear].x = x;+  q->items[q->rear].y = y;+  q->items[q->rear].label = label;+  q->size++;+  return 1;+}++static QueueItem dequeue(Queue *q) {+  QueueItem item = {-1, -1, -1};+  if (q->size == 0) {+    return item;+  }++  item = q->items[q->front];+  q->front = (q->front + 1) % q->capacity;+  q->size--;+  return item;+}++/**+ * Watershed segmentation+ * using (x, y) coordinate markers with elevation-based flooding+ *+ * Implements the watershed transform for image segmentation by treating the+ * grayscale image as an elevation map. Water floods from the marker points,+ * and watershed lines form where different regions would meet. Lower+ * intensity values represent valleys where water accumulates, and higher+ * values represent hills/ridges.+ *+ * @param width Width of the image.+ * @param height Height of the image.+ * @param grayscale_data Pointer to the input grayscale image data.+ * @param markers Array of Point2D markers with x,y coordinates.+ * @param num_markers Number of markers in the array.+ * @return Pointer to the segmented image data.+ */+uint8_t *fcv_watershed_segmentation(+  uint32_t width,+  uint32_t height,+  uint8_t const *const grayscale_data,+  Point2D *markers,+  uint32_t num_markers,+  bool create_boundaries+) {+  if (!grayscale_data || !markers || num_markers == 0) {+    return NULL;+  }++  if (width == 0 || height == 0) {+    return NULL;+  }++  // Validate all marker positions and coordinates+  for (uint32_t m = 0; m < num_markers; m++) {+    // Check for NaN/Inf+    if (!isfinite(markers[m].x) || !isfinite(markers[m].y)) {+      return NULL;+    }++    int32_t marker_x = (int32_t)markers[m].x;+    int32_t marker_y = (int32_t)markers[m].y;++    // Check bounds - if any marker is invalid, fail the entire operation+    if (marker_x < 0 || marker_x >= (int32_t)width || marker_y < 0 ||+        marker_y >= (int32_t)height) {+      return NULL;+    }+  }++  // Check for overflow: width * height+  if (width > SIZE_MAX / height) {+    return NULL;+  }+  size_t img_length_px = (size_t)width * height;++  // Check for overflow in output_data allocation (pixels * 4)+  if (img_length_px > SIZE_MAX / 4) {+    return NULL;+  }++  // Check for overflow in labels allocation+  if (img_length_px > SIZE_MAX / sizeof(int32_t)) {+    return NULL;+  }++  // Create output data+  uint8_t *output_data = malloc(img_length_px * 4);+  if (!output_data) {+    return NULL;+  }++  // Initialize labels array (-1 = unvisited, 0+ = region labels)+  int32_t *labels = malloc(img_length_px * sizeof(int32_t));+  if (!labels) {+    free(output_data);+    return NULL;+  }++  // Initialize all pixels as unvisited+  for (uint32_t i = 0; i < img_length_px; i++) {+    labels[i] = -1;+  }++  // Create queue for flood fill+  Queue *queue = create_queue(img_length_px);+  if (!queue) {+    free(output_data);+    free(labels);+    return NULL;+  }++  // Initialize markers and add to queue+  for (uint32_t m = 0; m < num_markers; m++) {+    int32_t marker_x = (int32_t)markers[m].x;+    int32_t marker_y = (int32_t)markers[m].y;+    int32_t idx = marker_y * width + marker_x;+    int32_t label = (int32_t)m; // Label starts from 0+    labels[idx] = label;+    enqueue(queue, marker_x, marker_y, label);+  }++  // Level-wise watershed algorithm - process pixels in elevation order+  int32_t dx[] = {-1, 1, 0, 0};+  int32_t dy[] = {0, 0, -1, 1};++  // Process each elevation level from 0 to 255+  for (int32_t current_level = 0; current_level <= 255; current_level++) {+    // Single-step expansion: each region grows by exactly one pixel layer per+    // iteration+    bool changed = true;+    while (changed) {+      changed = false;++      // Find all pixels at current level that should be added in this step+      for (uint32_t y = 0; y < height; y++) {+        for (uint32_t x = 0; x < width; x++) {+          int32_t idx = y * width + x;++          // Skip if already labeled+          if (labels[idx] != -1) {+            continue;+          }++          // Allow water to flow downhill:+          // Process pixels at or below current level.+          // This enables proper watershed behavior+          // where water flows from high to low elevation.+          bool can_flood = (grayscale_data[idx] <= current_level);++          if (!can_flood) {+            continue;+          }++          // Check if this pixel neighbors any labeled pixel+          int32_t neighbor_label = -1;+          bool multiple_labels = false;++          for (int32_t d = 0; d < 4; d++) {+            int32_t nx = x + dx[d];+            int32_t ny = y + dy[d];++            if (nx >= 0 && nx < (int32_t)width && ny >= 0 &&+                ny < (int32_t)height) {+              int32_t neighbor_idx = ny * width + nx;+              if (labels[neighbor_idx] != -1) {+                if (neighbor_label == -1) {+                  neighbor_label = labels[neighbor_idx];+                }+                else if (neighbor_label != labels[neighbor_idx]) {+                  multiple_labels = true;+                  break;+                }+              }+            }+          }++          // If found exactly one labeled neighbor, mark for expansion+          // If multiple labels and boundaries disabled, use first found label+          if (neighbor_label != -1 &&+              (!multiple_labels || !create_boundaries)) {+            enqueue(queue, x, y, neighbor_label);+            changed = true;+          }+        }+      }++      // Apply all expansions from this step simultaneously+      while (!is_queue_empty(queue)) {+        QueueItem item = dequeue(queue);+        int32_t idx = item.y * width + item.x;+        labels[idx] = item.label;+      }+    }+  }++  destroy_queue(queue);++  // Generate distinct colors for each region+  uint8_t colors[10][3] = {+    {255, 0, 0},     // red+    {0, 255, 0},     // green+    {0, 0, 255},     // blue+    {255, 255, 0},   // yellow+    {255, 0, 255},   // magenta+    {0, 255, 255},   // cyan+    {255, 128, 0},   // orange+    {128, 0, 255},   // purple+    {255, 192, 203}, // pink+    {128, 128, 128}  // gray+  };++  // Create output image+  for (uint32_t i = 0; i < img_length_px; i++) {+    uint32_t rgba_idx = i * 4;+    int32_t label = labels[i];++    if (label == -1) {+      // Unassigned pixels -> background (black)+      output_data[rgba_idx] = 0;     // R+      output_data[rgba_idx + 1] = 0; // G+      output_data[rgba_idx + 2] = 0; // B+    }+    else if (label >= 0 && label < 10) {+      // Assigned to a region+      output_data[rgba_idx] = colors[label][0];     // R+      output_data[rgba_idx + 1] = colors[label][1]; // G+      output_data[rgba_idx + 2] = colors[label][2]; // B+    }+    else {+      // Fallback for more than 10 regions+      output_data[rgba_idx] = 128;     // R+      output_data[rgba_idx + 1] = 128; // G+      output_data[rgba_idx + 2] = 128; // B+    }+    output_data[rgba_idx + 3] = 255; // A+  }++  free(labels);++  return output_data;+}+/* End of FlatCV amalgamation */

@@ -0,0 +1,610 @@+{-# LANGUAGE RecordWildCards #-}++-- |+-- Haskell bindings for the FlatCV image processing library.+--+-- FlatCV is a lightweight, dependency-free image processing library+-- useful for document scanning, image preprocessing, and basic+-- computer vision tasks.+--+-- = Usage Example+--+-- @+-- import FlatCV+-- import qualified Data.Vector.Storable as V+--+-- main :: IO ()+-- main = do+--     -- Load your image data as a Vector of Word8 (RGBA format)+--     let width = 640+--         height = 480+--         imageData = ... -- your RGBA pixel data+--+--     -- Convert to grayscale+--     grayData <- grayscale width height imageData+--+--     -- Apply Gaussian blur+--     blurredData <- gaussianBlur width height 2.0 imageData+-- @+--+module FlatCV+  ( -- * Types+    Image(..)+  , Corners(..)+  , Matrix3x3(..)+  , Point2D(..)++    -- * Pretty printing+  , prettyShowCorners+  , prettyShowMatrix3x3++    -- * Color conversion+  , grayscale+  , grayscaleStretch+  , rgbaToGrayscale+  , singleToMultichannel++    -- * Filtering+  , gaussianBlur++    -- * Thresholding+  , otsuThreshold+  , bwSmart+  , applyGlobalThreshold++    -- * Edge detection+  , sobelEdgeDetection+  , foerstnerCorner++    -- * Morphological operations+  , binaryDilationDisk+  , binaryErosionDisk+  , binaryClosingDisk+  , binaryOpeningDisk++    -- * Geometric transforms+  , flipX+  , flipY+  , transpose+  , transverse+  , rotate90CW+  , rotate180+  , rotate270CW+  , resize++    -- * Cropping and trimming+  , crop+  , trim+  , trimThreshold++    -- * Perspective transform+  , calculatePerspectiveTransform+  , applyMatrix3x3+  , detectCorners++    -- * Document extraction+  , extractDocument+  , extractDocumentAuto++    -- * Drawing+  , drawCircle+  , drawDisk+  , addBorder++    -- * Analysis+  , generateHistogram++    -- * Segmentation+  , watershedSegmentation++    -- * EXIF+  , getExifOrientation++    -- * Utility+  , parseHexColor+  , convertToBinary+  ) where++import Control.Monad (when)+import Data.Int (Int32)+import Data.Word (Word8, Word32)+import qualified Data.Vector.Storable as V+import qualified Data.Vector.Storable.Mutable as MV+import Foreign.C.String (withCString)+import Foreign.C.Types (CUChar)+import Foreign.ForeignPtr (newForeignPtr, FinalizerPtr, castForeignPtr)+import Foreign.Marshal.Alloc (alloca, free, mallocBytes)+import Foreign.Marshal.Array (pokeArray)+import Foreign.Ptr (Ptr, castPtr, nullPtr)+import Foreign.Storable (Storable(..), peek, poke, sizeOf, alignment)+import Text.Printf (printf)++import qualified FlatCV.Raw as Raw++-- | An image with dimensions and pixel data+data Image = Image+  { imageWidth    :: {-# UNPACK #-} !Word32+  , imageHeight   :: {-# UNPACK #-} !Word32+  , imageChannels :: {-# UNPACK #-} !Word32+  , imageData     :: !(V.Vector Word8)+  } deriving (Show, Eq)++-- | 2D point+data Point2D = Point2D+  { pointX :: {-# UNPACK #-} !Double+  , pointY :: {-# UNPACK #-} !Double+  } deriving (Show, Eq)++instance Storable Point2D where+  sizeOf _ = 2 * sizeOf (0.0 :: Double)+  alignment _ = alignment (0.0 :: Double)+  peek ptr = do+    x <- peekByteOff ptr 0+    y <- peekByteOff ptr 8+    return $ Point2D x y+  poke ptr Point2D{..} = do+    pokeByteOff ptr 0 pointX+    pokeByteOff ptr 8 pointY++-- | Quadrilateral corners (top-left, top-right, bottom-right, bottom-left)+data Corners = Corners+  { tlX :: {-# UNPACK #-} !Double+  , tlY :: {-# UNPACK #-} !Double+  , trX :: {-# UNPACK #-} !Double+  , trY :: {-# UNPACK #-} !Double+  , brX :: {-# UNPACK #-} !Double+  , brY :: {-# UNPACK #-} !Double+  , blX :: {-# UNPACK #-} !Double+  , blY :: {-# UNPACK #-} !Double+  } deriving (Show, Eq)++instance Storable Corners where+  sizeOf _ = 8 * sizeOf (0.0 :: Double)+  alignment _ = alignment (0.0 :: Double)+  peek ptr = do+    tlX <- peekByteOff ptr 0+    tlY <- peekByteOff ptr 8+    trX <- peekByteOff ptr 16+    trY <- peekByteOff ptr 24+    brX <- peekByteOff ptr 32+    brY <- peekByteOff ptr 40+    blX <- peekByteOff ptr 48+    blY <- peekByteOff ptr 56+    return Corners{..}+  poke ptr Corners{..} = do+    pokeByteOff ptr 0 tlX+    pokeByteOff ptr 8 tlY+    pokeByteOff ptr 16 trX+    pokeByteOff ptr 24 trY+    pokeByteOff ptr 32 brX+    pokeByteOff ptr 40 brY+    pokeByteOff ptr 48 blX+    pokeByteOff ptr 56 blY++-- | Pretty print corners for debugging+prettyShowCorners :: Corners -> String+prettyShowCorners Corners{..} =+  show (tlX, tlY) ++ " " ++ show (trX, trY) ++ "\n" +++  show (blX, blY) ++ " " ++ show (brX, brY)++-- | 3x3 transformation matrix+data Matrix3x3 = Matrix3x3+  { m00 :: {-# UNPACK #-} !Double+  , m01 :: {-# UNPACK #-} !Double+  , m02 :: {-# UNPACK #-} !Double+  , m10 :: {-# UNPACK #-} !Double+  , m11 :: {-# UNPACK #-} !Double+  , m12 :: {-# UNPACK #-} !Double+  , m20 :: {-# UNPACK #-} !Double+  , m21 :: {-# UNPACK #-} !Double+  , m22 :: {-# UNPACK #-} !Double+  } deriving (Show, Eq)++instance Storable Matrix3x3 where+  sizeOf _ = 9 * sizeOf (0.0 :: Double)+  alignment _ = alignment (0.0 :: Double)+  peek ptr = do+    m00 <- peekByteOff ptr 0+    m01 <- peekByteOff ptr 8+    m02 <- peekByteOff ptr 16+    m10 <- peekByteOff ptr 24+    m11 <- peekByteOff ptr 32+    m12 <- peekByteOff ptr 40+    m20 <- peekByteOff ptr 48+    m21 <- peekByteOff ptr 56+    m22 <- peekByteOff ptr 64+    return Matrix3x3{..}+  poke ptr Matrix3x3{..} = do+    pokeByteOff ptr 0 m00+    pokeByteOff ptr 8 m01+    pokeByteOff ptr 16 m02+    pokeByteOff ptr 24 m10+    pokeByteOff ptr 32 m11+    pokeByteOff ptr 40 m12+    pokeByteOff ptr 48 m20+    pokeByteOff ptr 56 m21+    pokeByteOff ptr 64 m22++-- | Pretty print 3x3 matrix for debugging+prettyShowMatrix3x3 :: Matrix3x3 -> String+prettyShowMatrix3x3 Matrix3x3{..} =+  let fNum n = printf "% .5f" (n :: Double)+  in fNum m00 ++ " " ++ fNum m01 ++ " " ++ fNum m02 ++ "\n" +++     fNum m10 ++ " " ++ fNum m11 ++ " " ++ fNum m12 ++ "\n" +++     fNum m20 ++ " " ++ fNum m21 ++ " " ++ fNum m22++-- Foreign finalizer for C-allocated memory+foreign import ccall "stdlib.h &free" finalizerFree :: FinalizerPtr CUChar++-- Helper to create a Vector from a C-allocated buffer+makeVector :: Ptr Word8 -> Int -> IO (V.Vector Word8)+makeVector outPtr len = do+  fp <- newForeignPtr finalizerFree (castPtr outPtr)+  return $ V.unsafeFromForeignPtr0 (castForeignPtr fp) len++-- Helper to run a function that returns a new image buffer+withImageData :: V.Vector Word8 -> (Ptr Word8 -> IO (Ptr Word8)) -> Word32 -> IO (V.Vector Word8)+withImageData input action outputLen = do+  V.unsafeWith input $ \inPtr -> do+    outPtr <- action inPtr+    when (outPtr == nullPtr) $ error "FlatCV: operation failed (null pointer returned)"+    makeVector outPtr (fromIntegral outputLen)++--------------------------------------------------------------------------------+-- Color conversion+--------------------------------------------------------------------------------++-- | Convert RGBA image to grayscale (output is still 4-channel RGBA)+grayscale :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+grayscale width height input =+  withImageData input (\p -> Raw.fcv_grayscale width height p) (width * height * 4)++-- | Convert to grayscale with contrast stretching+grayscaleStretch :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+grayscaleStretch width height input =+  withImageData input (\p -> Raw.fcv_grayscale_stretch width height p) (width * height * 4)++-- | Convert RGBA to single-channel grayscale+rgbaToGrayscale :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+rgbaToGrayscale width height input =+  withImageData input (\p -> Raw.fcv_rgba_to_grayscale width height p) (width * height)++-- | Convert single-channel grayscale to RGBA+singleToMultichannel :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+singleToMultichannel width height input =+  withImageData input (\p -> Raw.fcv_single_to_multichannel width height p) (width * height * 4)++--------------------------------------------------------------------------------+-- Filtering+--------------------------------------------------------------------------------++-- | Apply Gaussian blur with given radius+gaussianBlur :: Word32 -> Word32 -> Double -> V.Vector Word8 -> IO (V.Vector Word8)+gaussianBlur width height radius input =+  withImageData input (\p -> Raw.fcv_apply_gaussian_blur width height radius p) (width * height * 4)++--------------------------------------------------------------------------------+-- Thresholding+--------------------------------------------------------------------------------++-- | Apply Otsu threshold to convert image to binary+-- If useDoubleThreshold is True, uses double thresholding for better results+otsuThreshold :: Word32 -> Word32 -> Bool -> V.Vector Word8 -> IO (V.Vector Word8)+otsuThreshold width height useDouble input =+  withImageData input (\p -> Raw.fcv_otsu_threshold_rgba width height useDouble p) (width * height * 4)++-- | Smart black and white conversion with anti-aliasing+bwSmart :: Word32 -> Word32 -> Bool -> V.Vector Word8 -> IO (V.Vector Word8)+bwSmart width height useDouble input =+  withImageData input (\p -> Raw.fcv_bw_smart width height useDouble p) (width * height * 4)++-- | Apply global threshold in-place+applyGlobalThreshold :: Word32 -> Word8 -> V.Vector Word8 -> IO (V.Vector Word8)+applyGlobalThreshold imgLen threshold input = do+  output <- V.thaw input+  MV.unsafeWith output $ \ptr ->+    Raw.fcv_apply_global_threshold imgLen ptr threshold+  V.unsafeFreeze output++--------------------------------------------------------------------------------+-- Edge detection+--------------------------------------------------------------------------------++-- | Apply Sobel edge detection+sobelEdgeDetection :: Word32 -> Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+sobelEdgeDetection width height channels input =+  withImageData input (\p -> Raw.fcv_sobel_edge_detection width height channels p) (width * height * channels)++-- | Foerstner corner detection with given sigma+foerstnerCorner :: Word32 -> Word32 -> Double -> V.Vector Word8 -> IO (V.Vector Word8)+foerstnerCorner width height sigma input =+  withImageData input (\p -> Raw.fcv_foerstner_corner width height p sigma) (width * height * 4)++--------------------------------------------------------------------------------+-- Morphological operations+--------------------------------------------------------------------------------++-- | Binary dilation with disk structuring element+binaryDilationDisk :: Int32 -> Int32 -> Int32 -> V.Vector Word8 -> IO (V.Vector Word8)+binaryDilationDisk width height radius input =+  withImageData input (\p -> Raw.fcv_binary_dilation_disk p width height radius)+    (fromIntegral $ width * height * 4)++-- | Binary erosion with disk structuring element+binaryErosionDisk :: Int32 -> Int32 -> Int32 -> V.Vector Word8 -> IO (V.Vector Word8)+binaryErosionDisk width height radius input =+  withImageData input (\p -> Raw.fcv_binary_erosion_disk p width height radius)+    (fromIntegral $ width * height * 4)++-- | Binary closing (dilation followed by erosion)+binaryClosingDisk :: Int32 -> Int32 -> Int32 -> V.Vector Word8 -> IO (V.Vector Word8)+binaryClosingDisk width height radius input =+  withImageData input (\p -> Raw.fcv_binary_closing_disk p width height radius)+    (fromIntegral $ width * height * 4)++-- | Binary opening (erosion followed by dilation)+binaryOpeningDisk :: Int32 -> Int32 -> Int32 -> V.Vector Word8 -> IO (V.Vector Word8)+binaryOpeningDisk width height radius input =+  withImageData input (\p -> Raw.fcv_binary_opening_disk p width height radius)+    (fromIntegral $ width * height * 4)++--------------------------------------------------------------------------------+-- Geometric transforms+--------------------------------------------------------------------------------++-- | Flip image horizontally+flipX :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+flipX width height input =+  withImageData input (\p -> Raw.fcv_flip_x width height p) (width * height * 4)++-- | Flip image vertically+flipY :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+flipY width height input =+  withImageData input (\p -> Raw.fcv_flip_y width height p) (width * height * 4)++-- | Transpose image (swap rows and columns)+transpose :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+transpose width height input =+  withImageData input (\p -> Raw.fcv_transpose width height p) (width * height * 4)++-- | Transverse (transpose + flip both axes)+transverse :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+transverse width height input =+  withImageData input (\p -> Raw.fcv_transverse width height p) (width * height * 4)++-- | Rotate image 90 degrees clockwise+rotate90CW :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+rotate90CW width height input =+  withImageData input (\p -> Raw.fcv_rotate_90_cw width height p) (width * height * 4)++-- | Rotate image 180 degrees+rotate180 :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+rotate180 width height input =+  withImageData input (\p -> Raw.fcv_rotate_180 width height p) (width * height * 4)++-- | Rotate image 270 degrees clockwise (90 counter-clockwise)+rotate270CW :: Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+rotate270CW width height input =+  withImageData input (\p -> Raw.fcv_rotate_270_cw width height p) (width * height * 4)++-- | Resize image with bilinear interpolation+-- Returns (new width, new height, output data)+resize :: Word32 -> Word32 -> Double -> Double -> V.Vector Word8 -> IO (Word32, Word32, V.Vector Word8)+resize width height scaleX scaleY input = do+  alloca $ \outWPtr ->+    alloca $ \outHPtr -> do+      V.unsafeWith input $ \inPtr -> do+        outPtr <- Raw.fcv_resize width height scaleX scaleY outWPtr outHPtr inPtr+        when (outPtr == nullPtr) $ error "FlatCV: resize failed"+        outW <- peek outWPtr+        outH <- peek outHPtr+        outData <- makeVector outPtr (fromIntegral $ outW * outH * 4)+        return (outW, outH, outData)++--------------------------------------------------------------------------------+-- Cropping and trimming+--------------------------------------------------------------------------------++-- | Crop a rectangular region from the image+crop :: Word32 -> Word32 -> Word32 -> Word32 -> Word32 -> Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+crop width height channels x y newWidth newHeight input =+  withImageData input (\p -> Raw.fcv_crop width height channels p x y newWidth newHeight)+    (newWidth * newHeight * channels)++-- | Trim border pixels with same color+-- Returns (new width, new height, output data)+trim :: Int32 -> Int32 -> Word32 -> V.Vector Word8 -> IO (Int32, Int32, V.Vector Word8)+trim width height channels input = do+  alloca $ \wPtr ->+    alloca $ \hPtr -> do+      poke wPtr width+      poke hPtr height+      V.unsafeWith input $ \inPtr -> do+        outPtr <- Raw.fcv_trim wPtr hPtr channels inPtr+        when (outPtr == nullPtr) $ error "FlatCV: trim failed"+        outW <- peek wPtr+        outH <- peek hPtr+        outData <- makeVector outPtr (fromIntegral $ outW * outH * fromIntegral channels)+        return (outW, outH, outData)++-- | Trim border pixels within threshold percentage+-- Returns (new width, new height, output data)+trimThreshold :: Int32 -> Int32 -> Word32 -> Double -> V.Vector Word8 -> IO (Int32, Int32, V.Vector Word8)+trimThreshold width height channels thresholdPercent input = do+  alloca $ \wPtr ->+    alloca $ \hPtr -> do+      poke wPtr width+      poke hPtr height+      V.unsafeWith input $ \inPtr -> do+        outPtr <- Raw.fcv_trim_threshold wPtr hPtr channels inPtr thresholdPercent+        when (outPtr == nullPtr) $ error "FlatCV: trim_threshold failed"+        outW <- peek wPtr+        outH <- peek hPtr+        outData <- makeVector outPtr (fromIntegral $ outW * outH * fromIntegral channels)+        return (outW, outH, outData)++--------------------------------------------------------------------------------+-- Perspective transform+--------------------------------------------------------------------------------++-- | Calculate perspective transformation matrix from source to destination corners+calculatePerspectiveTransform :: Corners -> Corners -> IO Matrix3x3+calculatePerspectiveTransform src dst = do+  alloca $ \srcPtr ->+    alloca $ \dstPtr -> do+      poke srcPtr src+      poke dstPtr dst+      matPtr <- Raw.fcv_calculate_perspective_transform (castPtr srcPtr) (castPtr dstPtr)+      when (matPtr == nullPtr) $ error "FlatCV: perspective transform calculation failed"+      mat <- peek (castPtr matPtr)+      free matPtr+      return mat++-- | Apply 3x3 transformation matrix to image+applyMatrix3x3 :: Int32 -> Int32 -> Int32 -> Int32 -> Matrix3x3 -> V.Vector Word8 -> IO (V.Vector Word8)+applyMatrix3x3 inWidth inHeight outWidth outHeight matrix input = do+  alloca $ \matPtr -> do+    poke matPtr matrix+    withImageData input (\p -> Raw.fcv_apply_matrix_3x3 inWidth inHeight p outWidth outHeight (castPtr matPtr))+      (fromIntegral $ outWidth * outHeight * 4)++-- | Detect document corners in an image+detectCorners :: Int32 -> Int32 -> V.Vector Word8 -> IO Corners+detectCorners width height input = do+  V.unsafeWith input $ \inPtr -> do+    cornersPtr <- Raw.fcv_detect_corners_ptr inPtr width height+    when (cornersPtr == nullPtr) $ error "FlatCV: corner detection failed"+    corners <- peek (castPtr cornersPtr)+    free cornersPtr+    return corners++--------------------------------------------------------------------------------+-- Document extraction+--------------------------------------------------------------------------------++-- | Extract document with specified output dimensions+extractDocument :: Word32 -> Word32 -> Word32 -> Word32 -> V.Vector Word8 -> IO (V.Vector Word8)+extractDocument width height outWidth outHeight input =+  withImageData input (\p -> Raw.fcv_extract_document width height p outWidth outHeight)+    (outWidth * outHeight * 4)++-- | Extract document with automatic size detection+-- Returns (output width, output height, output data)+extractDocumentAuto :: Word32 -> Word32 -> V.Vector Word8 -> IO (Word32, Word32, V.Vector Word8)+extractDocumentAuto width height input = do+  alloca $ \outWPtr ->+    alloca $ \outHPtr -> do+      V.unsafeWith input $ \inPtr -> do+        outPtr <- Raw.fcv_extract_document_auto width height inPtr outWPtr outHPtr+        when (outPtr == nullPtr) $ error "FlatCV: document extraction failed"+        outW <- peek outWPtr+        outH <- peek outHPtr+        outData <- makeVector outPtr (fromIntegral $ outW * outH * 4)+        return (outW, outH, outData)++--------------------------------------------------------------------------------+-- Drawing+--------------------------------------------------------------------------------++-- | Draw a circle outline on an image (modifies input in place)+drawCircle :: Word32 -> Word32 -> Word32 -> String -> Double -> Double -> Double -> V.Vector Word8 -> IO (V.Vector Word8)+drawCircle width height channels hexColor radius centerX centerY input = do+  output <- V.thaw input+  MV.unsafeWith output $ \ptr ->+    withCString hexColor $ \colorPtr ->+      Raw.fcv_draw_circle width height channels colorPtr radius centerX centerY ptr+  V.unsafeFreeze output++-- | Draw a filled disk on an image (modifies input in place)+drawDisk :: Word32 -> Word32 -> Word32 -> String -> Double -> Double -> Double -> V.Vector Word8 -> IO (V.Vector Word8)+drawDisk width height channels hexColor radius centerX centerY input = do+  output <- V.thaw input+  MV.unsafeWith output $ \ptr ->+    withCString hexColor $ \colorPtr ->+      Raw.fcv_draw_disk width height channels colorPtr radius centerX centerY ptr+  V.unsafeFreeze output++-- | Add a border around an image+-- Returns (output width, output height, output data)+addBorder :: Word32 -> Word32 -> Word32 -> String -> Word32 -> V.Vector Word8 -> IO (Word32, Word32, V.Vector Word8)+addBorder width height channels hexColor borderWidth input = do+  alloca $ \outWPtr ->+    alloca $ \outHPtr -> do+      V.unsafeWith input $ \inPtr ->+        withCString hexColor $ \colorPtr -> do+          outPtr <- Raw.fcv_add_border width height channels colorPtr borderWidth inPtr outWPtr outHPtr+          when (outPtr == nullPtr) $ error "FlatCV: add_border failed"+          outW <- peek outWPtr+          outH <- peek outHPtr+          outData <- makeVector outPtr (fromIntegral $ outW * outH * channels)+          return (outW, outH, outData)++--------------------------------------------------------------------------------+-- Analysis+--------------------------------------------------------------------------------++-- | Generate histogram visualization+-- Returns (histogram width, histogram height, RGBA output data)+generateHistogram :: Word32 -> Word32 -> Word32 -> V.Vector Word8 -> IO (Word32, Word32, V.Vector Word8)+generateHistogram width height channels input = do+  alloca $ \outWPtr ->+    alloca $ \outHPtr -> do+      V.unsafeWith input $ \inPtr -> do+        outPtr <- Raw.fcv_generate_histogram width height channels inPtr outWPtr outHPtr+        when (outPtr == nullPtr) $ error "FlatCV: histogram generation failed"+        outW <- peek outWPtr+        outH <- peek outHPtr+        outData <- makeVector outPtr (fromIntegral $ outW * outH * 4)+        return (outW, outH, outData)++--------------------------------------------------------------------------------+-- Segmentation+--------------------------------------------------------------------------------++-- | Watershed segmentation with markers+watershedSegmentation :: Word32 -> Word32 -> [Point2D] -> Bool -> V.Vector Word8 -> IO (V.Vector Word8)+watershedSegmentation width height markers createBoundaries input = do+  let numMarkers = fromIntegral $ length markers+      markerSize = numMarkers * sizeOf (undefined :: Point2D)+  markersPtr <- mallocBytes markerSize+  pokeArray markersPtr markers+  result <- withImageData input+    (\p -> Raw.fcv_watershed_segmentation width height p (castPtr markersPtr) (fromIntegral numMarkers) createBoundaries)+    (width * height * 4)+  free markersPtr+  return result++--------------------------------------------------------------------------------+-- EXIF+--------------------------------------------------------------------------------++-- | Get EXIF orientation from a JPEG file+-- Returns orientation value 1-8, or 1 if not found+getExifOrientation :: FilePath -> IO Int32+getExifOrientation path =+  withCString path Raw.fcv_get_exif_orientation++--------------------------------------------------------------------------------+-- Utility+--------------------------------------------------------------------------------++-- | Parse hex color string (e.g., "#FF0000" or "FF0000") to RGB+parseHexColor :: String -> IO (Word8, Word8, Word8)+parseHexColor hexColor =+  alloca $ \rPtr ->+    alloca $ \gPtr ->+      alloca $ \bPtr ->+        withCString hexColor $ \colorPtr -> do+          Raw.fcv_parse_hex_color colorPtr rPtr gPtr bPtr+          r <- peek rPtr+          g <- peek gPtr+          b <- peek bPtr+          return (r, g, b)++-- | Convert image to binary with custom foreground and background colors+convertToBinary :: Int32 -> Int32 -> String -> String -> V.Vector Word8 -> IO (V.Vector Word8)+convertToBinary width height fgHex bgHex input =+  withCString fgHex $ \fgPtr ->+    withCString bgHex $ \bgPtr ->+      withImageData input (\p -> Raw.fcv_convert_to_binary p width height fgPtr bgPtr)+        (fromIntegral $ width * height * 4)

@@ -0,0 +1,515 @@+{-# LANGUAGE ForeignFunctionInterface #-}++module FlatCV.Raw where++import Data.Int (Int32)+import Data.Word (Word8, Word32)+import Foreign.C.String (CString)+import Foreign.C.Types (CInt, CUChar, CUInt)+import Foreign.Ptr (Ptr, castPtr)++#include "flatcv.h"++-- | Point2D structure+data Point2D = Point2D+  { point2D_x :: {-# UNPACK #-} !Double+  , point2D_y :: {-# UNPACK #-} !Double+  } deriving (Show, Eq)++{#pointer *Point2D as Point2DPtr -> Point2D#}++-- | Corners structure for quadrilateral coordinates+data Corners = Corners+  { corners_tl_x :: {-# UNPACK #-} !Double+  , corners_tl_y :: {-# UNPACK #-} !Double+  , corners_tr_x :: {-# UNPACK #-} !Double+  , corners_tr_y :: {-# UNPACK #-} !Double+  , corners_br_x :: {-# UNPACK #-} !Double+  , corners_br_y :: {-# UNPACK #-} !Double+  , corners_bl_x :: {-# UNPACK #-} !Double+  , corners_bl_y :: {-# UNPACK #-} !Double+  } deriving (Show, Eq)++{#pointer *Corners as CornersPtr -> Corners#}++-- | 3x3 transformation matrix+data Matrix3x3 = Matrix3x3+  { matrix_m00 :: {-# UNPACK #-} !Double+  , matrix_m01 :: {-# UNPACK #-} !Double+  , matrix_m02 :: {-# UNPACK #-} !Double+  , matrix_m10 :: {-# UNPACK #-} !Double+  , matrix_m11 :: {-# UNPACK #-} !Double+  , matrix_m12 :: {-# UNPACK #-} !Double+  , matrix_m20 :: {-# UNPACK #-} !Double+  , matrix_m21 :: {-# UNPACK #-} !Double+  , matrix_m22 :: {-# UNPACK #-} !Double+  } deriving (Show, Eq)++{#pointer *Matrix3x3 as Matrix3x3Ptr -> Matrix3x3#}++-- | Corner peaks result structure+data CornerPeaks = CornerPeaks+  { cornerPeaks_points :: Ptr Point2D+  , cornerPeaks_count  :: {-# UNPACK #-} !Word32+  } deriving (Show, Eq)++{#pointer *CornerPeaks as CornerPeaksPtr -> CornerPeaks#}++-- Marshalling helpers+toWord8Ptr :: Ptr CUChar -> Ptr Word8+toWord8Ptr = castPtr++fromWord8Ptr :: Ptr Word8 -> Ptr CUChar+fromWord8Ptr = castPtr++toInt32Ptr :: Ptr CInt -> Ptr Int32+toInt32Ptr = castPtr++fromInt32Ptr :: Ptr Int32 -> Ptr CInt+fromInt32Ptr = castPtr++toWord32Ptr :: Ptr CUInt -> Ptr Word32+toWord32Ptr = castPtr++fromWord32Ptr :: Ptr Word32 -> Ptr CUInt+fromWord32Ptr = castPtr++--------------------------------------------------------------------------------+-- Conversion functions+--------------------------------------------------------------------------------++-- | Apply Gaussian blur to an RGBA image+{#fun fcv_apply_gaussian_blur+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Double'                      -- ^ radius+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Convert RGBA image to grayscale (still 4-channel output)+{#fun fcv_grayscale+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Convert to grayscale with contrast stretching+{#fun fcv_grayscale_stretch+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Apply global threshold to image data (in-place)+{#fun fcv_apply_global_threshold+  { `Word32'                      -- ^ image length (width * height * channels)+  , fromWord8Ptr `Ptr Word8'      -- ^ data (modified in place)+  , `Word8'                       -- ^ threshold value+  } -> `()'+#}++-- | Apply Otsu threshold to RGBA image+{#fun fcv_otsu_threshold_rgba+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Bool'                        -- ^ use double threshold+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Smart black and white conversion+{#fun fcv_bw_smart+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Bool'                        -- ^ use double threshold+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Resize image with bilinear interpolation+{#fun fcv_resize+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Double'                      -- ^ scale_x+  , `Double'                      -- ^ scale_y+  , fromWord32Ptr `Ptr Word32'    -- ^ out_width pointer+  , fromWord32Ptr `Ptr Word32'    -- ^ out_height pointer+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Binary morphology functions+--------------------------------------------------------------------------------++-- | Binary dilation with disk structuring element+{#fun fcv_binary_dilation_disk+  { fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ width+  , `Int32'                       -- ^ height+  , `Int32'                       -- ^ radius+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Binary erosion with disk structuring element+{#fun fcv_binary_erosion_disk+  { fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ width+  , `Int32'                       -- ^ height+  , `Int32'                       -- ^ radius+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Binary closing (dilation then erosion) with disk+{#fun fcv_binary_closing_disk+  { fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ width+  , `Int32'                       -- ^ height+  , `Int32'                       -- ^ radius+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Binary opening (erosion then dilation) with disk+{#fun fcv_binary_opening_disk+  { fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ width+  , `Int32'                       -- ^ height+  , `Int32'                       -- ^ radius+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Convert to binary+--------------------------------------------------------------------------------++-- | Convert image to binary with custom colors+{#fun fcv_convert_to_binary+  { fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ width+  , `Int32'                       -- ^ height+  , id `CString'                  -- ^ foreground hex color+  , id `CString'                  -- ^ background hex color+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Corner detection+--------------------------------------------------------------------------------++-- | Detect corner peaks in an image+{#fun fcv_corner_peaks+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Word32'                      -- ^ min_distance+  , `Double'                      -- ^ accuracy_thresh+  , `Double'                      -- ^ roundness_thresh+  } -> `CornerPeaksPtr'+#}++--------------------------------------------------------------------------------+-- Crop+--------------------------------------------------------------------------------++-- | Crop a region from an image+{#fun fcv_crop+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Word32'                      -- ^ channels+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Word32'                      -- ^ x offset+  , `Word32'                      -- ^ y offset+  , `Word32'                      -- ^ new_width+  , `Word32'                      -- ^ new_height+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Drawing functions+--------------------------------------------------------------------------------++-- | Draw a circle outline on an image (in-place)+{#fun fcv_draw_circle+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Word32'                      -- ^ channels+  , id `CString'                  -- ^ hex color+  , `Double'                      -- ^ radius+  , `Double'                      -- ^ center_x+  , `Double'                      -- ^ center_y+  , fromWord8Ptr `Ptr Word8'      -- ^ data (modified in place)+  } -> `()'+#}++-- | Draw a filled disk on an image (in-place)+{#fun fcv_draw_disk+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Word32'                      -- ^ channels+  , id `CString'                  -- ^ hex color+  , `Double'                      -- ^ radius+  , `Double'                      -- ^ center_x+  , `Double'                      -- ^ center_y+  , fromWord8Ptr `Ptr Word8'      -- ^ data (modified in place)+  } -> `()'+#}++-- | Add a border around an image+{#fun fcv_add_border+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Word32'                      -- ^ channels+  , id `CString'                  -- ^ hex color+  , `Word32'                      -- ^ border_width+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , fromWord32Ptr `Ptr Word32'    -- ^ output_width pointer+  , fromWord32Ptr `Ptr Word32'    -- ^ output_height pointer+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- EXIF+--------------------------------------------------------------------------------++-- | Get EXIF orientation from a JPEG file+{#fun fcv_get_exif_orientation+  { id `CString'                  -- ^ filename+  } -> `Int32'+#}++--------------------------------------------------------------------------------+-- Document extraction+--------------------------------------------------------------------------------++-- | Extract document with fixed output size+{#fun fcv_extract_document+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Word32'                      -- ^ output_width+  , `Word32'                      -- ^ output_height+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Extract document with automatic size detection+{#fun fcv_extract_document_auto+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , fromWord32Ptr `Ptr Word32'    -- ^ output_width pointer+  , fromWord32Ptr `Ptr Word32'    -- ^ output_height pointer+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Flip and transpose+--------------------------------------------------------------------------------++-- | Flip image horizontally+{#fun fcv_flip_x+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Flip image vertically+{#fun fcv_flip_y+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Transpose image (swap x and y axes)+{#fun fcv_transpose+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Transverse (transpose + flip both axes)+{#fun fcv_transverse+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Rotation+--------------------------------------------------------------------------------++-- | Rotate image 90 degrees clockwise+{#fun fcv_rotate_90_cw+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Rotate image 180 degrees+{#fun fcv_rotate_180+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Rotate image 270 degrees clockwise (90 counter-clockwise)+{#fun fcv_rotate_270_cw+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Foerstner corner detection+--------------------------------------------------------------------------------++-- | Foerstner corner detection+{#fun fcv_foerstner_corner+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Double'                      -- ^ sigma+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Histogram+--------------------------------------------------------------------------------++-- | Generate histogram visualization+{#fun fcv_generate_histogram+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Word32'                      -- ^ channels+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , fromWord32Ptr `Ptr Word32'    -- ^ out_width pointer+  , fromWord32Ptr `Ptr Word32'    -- ^ out_height pointer+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Hex color parsing+--------------------------------------------------------------------------------++-- | Parse hex color string to RGB components+{#fun fcv_parse_hex_color+  { id `CString'                  -- ^ hex color string+  , fromWord8Ptr `Ptr Word8'      -- ^ r output+  , fromWord8Ptr `Ptr Word8'      -- ^ g output+  , fromWord8Ptr `Ptr Word8'      -- ^ b output+  } -> `()'+#}++--------------------------------------------------------------------------------+-- Perspective transform+--------------------------------------------------------------------------------++-- | Calculate perspective transform matrix+{#fun fcv_calculate_perspective_transform+  { `CornersPtr'                  -- ^ source corners+  , `CornersPtr'                  -- ^ destination corners+  } -> `Matrix3x3Ptr'+#}++-- | Apply 3x3 transformation matrix to image+{#fun fcv_apply_matrix_3x3+  { `Int32'                       -- ^ input width+  , `Int32'                       -- ^ input height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ output width+  , `Int32'                       -- ^ output height+  , `Matrix3x3Ptr'                -- ^ transformation matrix+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Detect document corners in an image+{#fun fcv_detect_corners_ptr+  { fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Int32'                       -- ^ width+  , `Int32'                       -- ^ height+  } -> `CornersPtr'+#}++--------------------------------------------------------------------------------+-- RGBA to grayscale+--------------------------------------------------------------------------------++-- | Convert RGBA to single-channel grayscale+{#fun fcv_rgba_to_grayscale+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Single to multichannel+--------------------------------------------------------------------------------++-- | Convert single-channel to RGBA+{#fun fcv_single_to_multichannel+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Sobel edge detection+--------------------------------------------------------------------------------++-- | Apply Sobel edge detection+{#fun fcv_sobel_edge_detection+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , `Word32'                      -- ^ channels+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Trim+--------------------------------------------------------------------------------++-- | Trim border pixels with same color+{#fun fcv_trim+  { fromInt32Ptr `Ptr Int32'      -- ^ width pointer (modified)+  , fromInt32Ptr `Ptr Int32'      -- ^ height pointer (modified)+  , `Word32'                      -- ^ channels+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  } -> `Ptr Word8' toWord8Ptr+#}++-- | Trim border pixels within threshold+{#fun fcv_trim_threshold+  { fromInt32Ptr `Ptr Int32'      -- ^ width pointer (modified)+  , fromInt32Ptr `Ptr Int32'      -- ^ height pointer (modified)+  , `Word32'                      -- ^ channels+  , fromWord8Ptr `Ptr Word8'      -- ^ input data+  , `Double'                      -- ^ threshold_percent+  } -> `Ptr Word8' toWord8Ptr+#}++--------------------------------------------------------------------------------+-- Watershed segmentation+--------------------------------------------------------------------------------++-- | Watershed segmentation+{#fun fcv_watershed_segmentation+  { `Word32'                      -- ^ width+  , `Word32'                      -- ^ height+  , fromWord8Ptr `Ptr Word8'      -- ^ grayscale data+  , `Point2DPtr'                  -- ^ markers array+  , `Word32'                      -- ^ num_markers+  , `Bool'                        -- ^ create_boundaries+  } -> `Ptr Word8' toWord8Ptr+#}

@@ -0,0 +1,250 @@+{-# LANGUAGE OverloadedStrings #-}++module Main where++import Data.Word (Word8)+import qualified Data.Vector.Storable as V+import Test.Hspec++import FlatCV++-- | Create a simple test image (4x4 RGBA)+testImage :: V.Vector Word8+testImage = V.fromList+  [ 255, 0, 0, 255,    100, 0, 0, 255,    50, 0, 0, 255,     0, 0, 0, 255+  , 0, 255, 0, 255,    100, 100, 0, 255,  50, 50, 0, 255,    0, 50, 0, 255+  , 0, 0, 255, 255,    0, 100, 100, 255,  0, 50, 50, 255,    0, 0, 50, 255+  , 255, 255, 255, 255, 200, 200, 200, 255, 100, 100, 100, 255, 0, 0, 0, 255+  ]++-- | Create a white image with a black border for trim testing+borderedImage :: V.Vector Word8+borderedImage = V.fromList+  [ 0, 0, 0, 255,      0, 0, 0, 255,      0, 0, 0, 255,      0, 0, 0, 255+  , 0, 0, 0, 255,      255, 255, 255, 255, 255, 255, 255, 255, 0, 0, 0, 255+  , 0, 0, 0, 255,      255, 255, 255, 255, 255, 255, 255, 255, 0, 0, 0, 255+  , 0, 0, 0, 255,      0, 0, 0, 255,      0, 0, 0, 255,      0, 0, 0, 255+  ]++main :: IO ()+main = hspec $ do+  describe "FlatCV" $ do+    describe "grayscale" $ do+      it "converts an RGBA image to grayscale" $ do+        result <- grayscale 4 4 testImage+        V.length result `shouldBe` 64  -- 4x4x4 channels++      it "produces valid grayscale values" $ do+        result <- grayscale 4 4 testImage+        -- First pixel (255, 0, 0) should become grayish+        -- Gray = 0.299*R + 0.587*G + 0.114*B ≈ 76+        let r = result V.! 0+            g = result V.! 1+            b = result V.! 2+        -- R, G, B should all be equal for grayscale+        r `shouldBe` g+        g `shouldBe` b++    describe "gaussianBlur" $ do+      it "returns an image of the same size" $ do+        result <- gaussianBlur 4 4 1.0 testImage+        V.length result `shouldBe` 64++      it "smooths the image" $ do+        result <- gaussianBlur 4 4 2.0 testImage+        -- The result should be different from the input+        result `shouldNotBe` testImage++    describe "flipX" $ do+      it "returns an image of the same size" $ do+        result <- flipX 4 4 testImage+        V.length result `shouldBe` 64++      it "flipping twice returns the original" $ do+        once <- flipX 4 4 testImage+        twice <- flipX 4 4 once+        twice `shouldBe` testImage++    describe "flipY" $ do+      it "returns an image of the same size" $ do+        result <- flipY 4 4 testImage+        V.length result `shouldBe` 64++      it "flipping twice returns the original" $ do+        once <- flipY 4 4 testImage+        twice <- flipY 4 4 once+        twice `shouldBe` testImage++    describe "rotate180" $ do+      it "returns an image of the same size" $ do+        result <- rotate180 4 4 testImage+        V.length result `shouldBe` 64++      it "rotating twice returns the original" $ do+        once <- rotate180 4 4 testImage+        twice <- rotate180 4 4 once+        twice `shouldBe` testImage++    describe "rotate90CW" $ do+      it "returns an image of the same size for square images" $ do+        result <- rotate90CW 4 4 testImage+        V.length result `shouldBe` 64++      it "rotating four times returns the original" $ do+        r1 <- rotate90CW 4 4 testImage+        r2 <- rotate90CW 4 4 r1+        r3 <- rotate90CW 4 4 r2+        r4 <- rotate90CW 4 4 r3+        r4 `shouldBe` testImage++    describe "crop" $ do+      it "crops a region from the image" $ do+        result <- crop 4 4 4 1 1 2 2 testImage+        V.length result `shouldBe` 16  -- 2x2x4 channels++    describe "otsuThreshold" $ do+      it "returns a binary image" $ do+        result <- otsuThreshold 4 4 False testImage+        V.length result `shouldBe` 64++    describe "bwSmart" $ do+      it "returns an anti-aliased binary image" $ do+        result <- bwSmart 4 4 False testImage+        V.length result `shouldBe` 64++    describe "sobelEdgeDetection" $ do+      it "detects edges in an image" $ do+        result <- sobelEdgeDetection 4 4 4 testImage+        V.length result `shouldBe` 64++    describe "rgbaToGrayscale" $ do+      it "converts to single-channel grayscale" $ do+        result <- rgbaToGrayscale 4 4 testImage+        V.length result `shouldBe` 16  -- 4x4x1 channel++    describe "singleToMultichannel" $ do+      it "converts single-channel to RGBA" $ do+        gray <- rgbaToGrayscale 4 4 testImage+        result <- singleToMultichannel 4 4 gray+        V.length result `shouldBe` 64++    describe "parseHexColor" $ do+      it "parses hex color with hash" $ do+        (r, g, b) <- parseHexColor "#FF0000"+        r `shouldBe` 255+        g `shouldBe` 0+        b `shouldBe` 0++      it "parses hex color without hash" $ do+        (r, g, b) <- parseHexColor "00FF00"+        r `shouldBe` 0+        g `shouldBe` 255+        b `shouldBe` 0++      it "parses blue color" $ do+        (r, g, b) <- parseHexColor "0000FF"+        r `shouldBe` 0+        g `shouldBe` 0+        b `shouldBe` 255++    describe "resize" $ do+      it "scales up an image" $ do+        (outW, outH, result) <- resize 4 4 2.0 2.0 testImage+        outW `shouldBe` 8+        outH `shouldBe` 8+        V.length result `shouldBe` 256  -- 8x8x4++      it "scales down an image" $ do+        (outW, outH, result) <- resize 4 4 0.5 0.5 testImage+        outW `shouldBe` 2+        outH `shouldBe` 2+        V.length result `shouldBe` 16  -- 2x2x4++    describe "trim" $ do+      it "trims border pixels" $ do+        (outW, outH, result) <- trim 4 4 4 borderedImage+        -- The trim function finds the bounding box of non-border pixels+        -- The inner 2x2 white area gets trimmed to remove the black border+        outW `shouldSatisfy` (> 0)+        outH `shouldSatisfy` (> 0)+        V.length result `shouldBe` fromIntegral (outW * outH * 4)++    describe "Corners" $ do+      it "has correct Storable instance" $ do+        let corners = Corners 0 1 2 3 4 5 6 7+        tlX corners `shouldBe` 0+        tlY corners `shouldBe` 1+        trX corners `shouldBe` 2+        trY corners `shouldBe` 3+        brX corners `shouldBe` 4+        brY corners `shouldBe` 5+        blX corners `shouldBe` 6+        blY corners `shouldBe` 7++    describe "prettyShowCorners" $ do+      it "formats corners as a readable string" $ do+        let corners = Corners 0 0 100 0 100 100 0 100+        prettyShowCorners corners `shouldBe` "(0.0,0.0) (100.0,0.0)\n(0.0,100.0) (100.0,100.0)"++    describe "Matrix3x3" $ do+      it "has correct Storable instance" $ do+        let mat = Matrix3x3 1 0 0 0 1 0 0 0 1  -- Identity matrix+        m00 mat `shouldBe` 1+        m11 mat `shouldBe` 1+        m22 mat `shouldBe` 1+        m01 mat `shouldBe` 0++    describe "prettyShowMatrix3x3" $ do+      it "formats matrix as a readable string" $ do+        let mat = Matrix3x3 1 0 0 0 1 0 0 0 1+        -- Just check it contains the values and has newlines+        prettyShowMatrix3x3 mat `shouldContain` "1.00000"+        prettyShowMatrix3x3 mat `shouldContain` "\n"++    describe "calculatePerspectiveTransform" $ do+      it "calculates identity-like transform for same corners" $ do+        let src = Corners 0 0 100 0 100 100 0 100+            dst = Corners 0 0 100 0 100 100 0 100+        mat <- calculatePerspectiveTransform src dst+        -- Diagonal should be close to 1 for identity-like transform+        abs (m00 mat - 1) `shouldSatisfy` (< 0.001)+        abs (m11 mat - 1) `shouldSatisfy` (< 0.001)++    describe "binaryDilationDisk" $ do+      it "returns an image of the same size" $ do+        result <- binaryDilationDisk 4 4 1 testImage+        V.length result `shouldBe` 64++    describe "binaryErosionDisk" $ do+      it "returns an image of the same size" $ do+        result <- binaryErosionDisk 4 4 1 testImage+        V.length result `shouldBe` 64++    describe "binaryClosingDisk" $ do+      it "returns an image of the same size" $ do+        result <- binaryClosingDisk 4 4 1 testImage+        V.length result `shouldBe` 64++    describe "binaryOpeningDisk" $ do+      it "returns an image of the same size" $ do+        result <- binaryOpeningDisk 4 4 1 testImage+        V.length result `shouldBe` 64++    describe "generateHistogram" $ do+      it "generates a histogram image" $ do+        (outW, outH, result) <- generateHistogram 4 4 4 testImage+        outW `shouldSatisfy` (> 0)+        outH `shouldSatisfy` (> 0)+        V.length result `shouldBe` fromIntegral (outW * outH * 4)++    describe "addBorder" $ do+      it "adds a border around the image" $ do+        (outW, outH, result) <- addBorder 4 4 4 "#FF0000" 2 testImage+        outW `shouldBe` 8  -- 4 + 2*2+        outH `shouldBe` 8+        V.length result `shouldBe` 256  -- 8x8x4++    describe "convertToBinary" $ do+      it "converts to binary with custom colors" $ do+        result <- convertToBinary 4 4 "#FFFFFF" "#000000" testImage+        V.length result `shouldBe` 64