solveKaggle.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a standalone notebook to solve a jigsaw puzzle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install matplotlib opencv-python scipy tqdm;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import cv2\n",
    "import math\n",
    "import random\n",
    "import itertools\n",
    "import json\n",
    "import tqdm.notebook as tqdm\n",
    "from functools import cache\n",
    "from scipy.signal import savgol_filter\n",
    "from scipy.signal import find_peaks\n",
    "import scipy.optimize\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Add utilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def imshow(img, title=None):\n",
    "    plt.title(title)\n",
    "    plt.imshow(img[0:1680, 0:1700])\n",
    "    plt.show()\n",
    "\n",
    "class Item():\n",
    "    def __init__(self, **kwargs):\n",
    "        self.update(**kwargs)\n",
    "\n",
    "    def update(self, **kwargs):\n",
    "        self.__dict__.update(kwargs)\n",
    "\n",
    "class LoopingList(list):\n",
    "    def __getitem__(self, i):\n",
    "        if isinstance(i, int):\n",
    "            return super().__getitem__(i % len(self))\n",
    "        else:\n",
    "            return super().__getitem__(i)\n",
    "\n",
    "def plot_contour(contour, **kwargs):\n",
    "    plt.plot(contour[:, :, 0], contour[:, :, 1], **kwargs)\n",
    "\n",
    "def fill_contour(contour, **kwargs):\n",
    "    plt.fill(contour[:, :, 0], contour[:, :, 1], **kwargs)\n",
    "\n",
    "def transform_contour(contour, idx, params):\n",
    "    x, y, degrees = params\n",
    "    matrix = cv2.getRotationMatrix2D(contour[idx][0], degrees, 1)\n",
    "    translate = (x, y) - contour[idx][0]\n",
    "    return cv2.transform(contour, matrix) + translate\n",
    "\n",
    "def sub_contour(c, idx0, idx1):\n",
    "    if idx1 > idx0:\n",
    "        return c[idx0:idx1]\n",
    "    else:\n",
    "        return np.concatenate([c[idx0:], c[:idx1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detect pieces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = []\n",
    "pieces = []\n",
    "for filename in tqdm.tqdm(sorted(glob.glob(\"scans/kaggle/*.jpg\"))):\n",
    "    img_original = cv2.imread(filename)\n",
    "    h, w = img_original.shape[:2]\n",
    "    img_gray = cv2.cvtColor(img_original, cv2.COLOR_BGR2GRAY)\n",
    "    _, img_binary = cv2.threshold(img_gray, 30, 255, cv2.THRESH_BINARY)\n",
    "    # no blur, it causes more harm than good for these images\n",
    "    raw_contours, _ = cv2.findContours(img_binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "    raw_contours = [c for c in raw_contours if cv2.contourArea(c) > 100e3]\n",
    "    # draw fully filled contours\n",
    "    img_contours = np.zeros((h, w), dtype=np.uint8)\n",
    "    for contour in raw_contours:\n",
    "        cv2.drawContours(img_contours, [contour], 0, (255, 255, 255), -1)\n",
    "    # remove small connected dirt\n",
    "    img_corrected = cv2.morphologyEx(img_contours, cv2.MORPH_OPEN, np.ones((9, 9), dtype=np.uint8))\n",
    "    # get the clean contour\n",
    "    contours, _ = cv2.findContours(img_corrected, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "    contours = [c for c in contours if cv2.contourArea(c) > 100e3]\n",
    "    img_masked = img_gray & img_corrected\n",
    "    img_pieces = []\n",
    "    for contour in contours:\n",
    "        x, y, w, h = cv2.boundingRect(contour)\n",
    "        piece = Item(\n",
    "            contour=contour - (x, y),\n",
    "            img=img_masked[y:y+h, x:x+w],\n",
    "            filename=filename,\n",
    "            filepos=(x, y)\n",
    "        )\n",
    "        img_pieces.append(piece)\n",
    "    \n",
    "    pieces.extend(img_pieces)\n",
    "    \n",
    "    image = Item(\n",
    "        filename=filename,\n",
    "        img_gray=img_gray,\n",
    "        pieces=img_pieces,\n",
    "    )\n",
    "    images.append(image)\n",
    "\n",
    "for id, piece in enumerate(pieces):\n",
    "    piece.update(id=id)\n",
    "\n",
    "print(f\"Detected {len(pieces)} pieces\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show the smallest and biggest pieces by area\n",
    "pieces.sort(key= lambda piece: cv2.contourArea(piece.contour))\n",
    "\n",
    "for piece in pieces[:1] + pieces[-1:]:\n",
    "    img = np.copy(piece.img)\n",
    "    plt.title(f\"Area={int(cv2.contourArea(piece.contour))}\")\n",
    "    plt.imshow(img)\n",
    "    plot_contour(piece.contour, c=\"red\", ls=\"--\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detect piece corners"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find corners via peak distance from center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for piece in tqdm.tqdm(pieces):\n",
    "    (cx, cy), r = cv2.minEnclosingCircle(piece.contour)\n",
    "    centered_contour = piece.contour - np.array([cx, cy])\n",
    "    # ensure peaks are not at start or end of the distances array\n",
    "    distances = np.sum(centered_contour**2, axis=2)[:, 0]\n",
    "    distance_offset = np.argmin(distances)\n",
    "    distances = np.concatenate([distances[distance_offset:], distances[:distance_offset]])\n",
    "\n",
    "    # find peak distances\n",
    "    peak_indices = [(distance_idx + distance_offset) % len(distances) for distance_idx in find_peaks(distances, prominence=10000)[0]]    \n",
    "    peak_indices.sort()\n",
    "    piece.update(peak_indices=LoopingList(peak_indices))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the pieces having the smallest / highest number of peak indices\n",
    "pieces.sort(key= lambda piece: len(piece.peak_indices))\n",
    "\n",
    "for piece in pieces[:1] + pieces[-1:]:\n",
    "    img = np.copy(piece.img)\n",
    "    plt.title(f\"Number of peaks={len(piece.peak_indices)}\")\n",
    "    for idx in piece.peak_indices:\n",
    "        cv2.circle(img, piece.contour[idx, 0], 10, (255, 0, 0), 3)\n",
    "    plt.imshow(img)\n",
    "    plot_contour(piece.contour, c=\"red\", ls=\"--\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter corners by angle geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for piece in tqdm.tqdm(pieces):\n",
    "    def is_right_angle(idx):\n",
    "        N = len(piece.contour)\n",
    "        angle_contour = np.concatenate([piece.contour[idx:], piece.contour[:idx]]) - piece.contour[idx]\n",
    "        contour0 = angle_contour[math.ceil(N*0.005):int(N*0.04)]\n",
    "        contour1 = angle_contour[-int(N*0.04):-math.ceil(N*0.005)]\n",
    "        angles0 = np.arctan2(contour0[:, 0, 1], contour0[:, 0, 0]) % (2 * np.pi)\n",
    "        angles1 = np.arctan2(contour1[:, 0, 1], contour1[:, 0, 0]) % (2 * np.pi)\n",
    "        mean0 = scipy.stats.circmean(angles0)\n",
    "        mean1 = scipy.stats.circmean(angles1)\n",
    "        std0 = scipy.stats.circstd(angles0)\n",
    "        std1 = scipy.stats.circstd(angles1)\n",
    "        angle = (mean0 - mean1) % np.pi\n",
    "        return -0.3 < angle - np.pi/2 < 0.3 and std0 < 0.1 and std1 < 0.1\n",
    "\n",
    "    # find peaks that look like right angles\n",
    "    angle_indices = [idx for idx in piece.peak_indices if is_right_angle(idx)]\n",
    "    piece.update(angle_indices=LoopingList(angle_indices))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the pieces having the smallest / highest number of angle indices\n",
    "pieces.sort(key= lambda piece: len(piece.angle_indices))\n",
    "\n",
    "for piece in pieces[:1] + pieces[-1:]:\n",
    "    img = np.copy(piece.img)\n",
    "    plt.title(f\"Number of peaks={len(piece.angle_indices)}\")\n",
    "    for idx in piece.angle_indices:\n",
    "        cv2.circle(img, piece.contour[idx, 0], 10, (255, 0, 0), 3)\n",
    "    plt.imshow(img)\n",
    "    plot_contour(piece.contour, c=\"red\", ls=\"--\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter corners by rectangle geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for piece in tqdm.tqdm(pieces):\n",
    "    def compute_rectangle_error(indices):\n",
    "            # get coordinates of corners\n",
    "            corners = LoopingList(np.take(piece.contour, sorted(list(indices)), axis=0)[:, 0, :])\n",
    "            # compute the side lengths and diagonal lengths\n",
    "            lengths = [math.sqrt(np.sum((corners[i0] - corners[i1])**2)) for i0, i1 in [(0, 1), (1, 2), (2, 3), (3, 0), (0, 2), (1, 3)]]\n",
    "            def f_error(a, b):\n",
    "                return abs(b - a) / (a + b)\n",
    "            return sum([f_error(lengths[i], lengths[j]) for i, j in [(0, 2), (1, 3), (4, 5)]])\n",
    "\n",
    "    # form a good rectangle with corner indices and missing indices\n",
    "    rectangles = []  # list of (score, [indices])\n",
    "    for indices in itertools.combinations(piece.peak_indices, 4):\n",
    "        if set(indices).issuperset(piece.angle_indices) or len(piece.angle_indices) > 4:\n",
    "            error = compute_rectangle_error(indices)\n",
    "            rectangles.append((error, indices))\n",
    "\n",
    "    if len(rectangles) >= 1:\n",
    "        error, indices = sorted(rectangles)[0]\n",
    "    else:\n",
    "        error, indices = [100, piece.angle_indices]\n",
    "    piece.update(rectangle_error=error)\n",
    "    piece.update(corner_indices=LoopingList(indices))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the pieces having the best / worst rectangle\n",
    "pieces.sort(key= lambda piece: piece.rectangle_error)\n",
    "\n",
    "for piece in pieces[:1] + pieces[-1:]:\n",
    "    img = np.copy(piece.img)\n",
    "    plt.title(f\"Rectangle error={piece.rectangle_error}\")\n",
    "    for idx in piece.corner_indices:\n",
    "        cv2.circle(img, piece.contour[idx, 0], 10, (255, 0, 0), 3)\n",
    "    plt.imshow(img)\n",
    "    plot_contour(piece.contour, c=\"red\", ls=\"--\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save/Restore snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(pieces) == 1998\n",
    "\n",
    "np.savez_compressed(\"kaggle_piece.contour.npz\", **dict([(str(idx), piece.contour) for idx, piece in enumerate(pieces)]))\n",
    "np.savez_compressed(\"kaggle_piece.img.npz\", **dict([(str(idx), piece.img) for idx, piece in enumerate(pieces)]))\n",
    "\n",
    "jpieces = []\n",
    "for idx, piece in enumerate(pieces):\n",
    "    jpieces.append(dict(\n",
    "        id=piece.id,\n",
    "        filename=piece.filename,\n",
    "        filepos=[int(x) for x in piece.filepos],\n",
    "        corner_indices=[int(x) for x in piece.corner_indices],\n",
    "    ))\n",
    "\n",
    "with open(\"kaggle_piece.json\", \"w\") as f:\n",
    "    json.dump(jpieces, f)\n",
    "\n",
    "print(f\"Dumped {len(jpieces)} pieces to files\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pieces = []\n",
    "\n",
    "piece_contours = np.load(\"kaggle_piece.contour.npz\")\n",
    "piece_imgs = np.load(\"kaggle_piece.img.npz\")\n",
    "\n",
    "with open(\"kaggle_piece.json\", \"r\") as f:\n",
    "    jpieces = json.load(f)\n",
    "    for idx, jpiece in enumerate(jpieces):\n",
    "        contour = piece_contours[str(idx)]\n",
    "        img = piece_imgs[str(idx)]\n",
    "        piece = Item(\n",
    "            id=jpiece['id'],\n",
    "            filename=jpiece['filename'],\n",
    "            filepos=tuple(jpiece['filepos']),\n",
    "            corner_indices=LoopingList([np.int64(x) for x in jpiece['corner_indices']]),\n",
    "            img=img,\n",
    "            contour=contour,\n",
    "        )\n",
    "        pieces.append(piece)\n",
    "\n",
    "print(f\"Loaded {len(pieces)} pieces from files\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyze edges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract edges & detect sign"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for piece in pieces:\n",
    "    edges = LoopingList()\n",
    "    contour = piece.contour.astype(np.float64)  # convert to float for rotation\n",
    "    for quarter in range(4):\n",
    "        idx0 = piece.corner_indices[quarter]\n",
    "        idx1 = piece.corner_indices[quarter+1]\n",
    "        p0 = contour[idx0][0]\n",
    "        p1 = contour[idx1][0]\n",
    "        # normalize the contour: first point at (0, 0), last point at (X, 0)\n",
    "        dx, dy = p1 - p0\n",
    "        length=math.sqrt(dx**2 + dy**2)\n",
    "        angle_radians = math.atan2(dy, dx)\n",
    "        matrix = cv2.getRotationMatrix2D(p0, math.degrees(angle_radians), 1)\n",
    "        normalized_piece_contour = cv2.transform(contour, matrix) - p0\n",
    "        normalized_edge_contour = sub_contour(normalized_piece_contour, idx0, idx1 + 1)\n",
    "\n",
    "        # compute the sign of the edge\n",
    "        heights = normalized_edge_contour[:, 0, 1]\n",
    "        if np.max(np.abs(heights)) > 20:\n",
    "            sign = 1 if np.max(heights) > - np.min(heights) else -1\n",
    "        else:\n",
    "            sign = 0\n",
    "\n",
    "        edge = Item(\n",
    "            idx0=idx0,\n",
    "            idx1=idx1,\n",
    "            normalized_piece_contour=normalized_piece_contour,\n",
    "            sign=sign,\n",
    "            length=length,\n",
    "        )\n",
    "        edges.append(edge)\n",
    "    \n",
    "    for idx, edge in enumerate(edges):\n",
    "        edge.update(\n",
    "            prev=edges[idx-1],\n",
    "            next=edges[idx+1]\n",
    "        )\n",
    "    \n",
    "    piece.update(\n",
    "        edges=edges,\n",
    "        nb_flats=len([edge for edge in edges if edge.sign == 0])\n",
    "    )\n",
    "\n",
    "print(\"edge sign:\", Counter([edge.sign for piece in pieces for edge in piece.edges]))\n",
    "print(\"nb of flats:\", Counter([piece.nb_flats for piece in pieces]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the pieces having the smallest / highest number of flats\n",
    "pieces.sort(key= lambda piece: piece.nb_flats)\n",
    "\n",
    "sign2color = {-1: \"red\", 0: \"green\", 1: \"yellow\"}\n",
    "\n",
    "for piece in pieces[:1] + pieces[-1:]:\n",
    "    img = np.copy(piece.img)\n",
    "    plt.title(f\"Nb of flats={piece.nb_flats}\")\n",
    "    for edge in piece.edges:\n",
    "        plot_contour(sub_contour(piece.contour, edge.idx0, edge.idx1), c=sign2color[edge.sign])\n",
    "    plt.imshow(img)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the pieces having the smallest / longest edge\n",
    "edge_pieces = [(edge, piece) for piece in pieces for edge in piece.edges]\n",
    "edge_pieces.sort(key= lambda ep: ep[0].length)\n",
    "\n",
    "for edge, piece in edge_pieces[:1] + edge_pieces[-1:]:\n",
    "    img = np.copy(piece.img)\n",
    "    plt.title(f\"Edge length={edge.length}\")\n",
    "    plot_contour(sub_contour(piece.contour, edge.idx0, edge.idx1), c='red')\n",
    "    plt.imshow(img)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute puzzle size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_size(area, perimeter):\n",
    "    # perimeter = 2 * (H+W)\n",
    "    # area = H*W\n",
    "    # H**2 - perimeter/2 * H + area = 0\n",
    "    a = 1\n",
    "    b = -perimeter/2\n",
    "    c = area\n",
    "    delta = b**2 - 4*a*c\n",
    "    h = int((-b - math.sqrt(delta)) / (2*a))\n",
    "    w = int((-b + math.sqrt(delta)) / (2*a))\n",
    "    return (min(h, w), max(h, w))\n",
    "\n",
    "nb_flats = Counter([piece.nb_flats for piece in pieces])\n",
    "assert nb_flats[2] == 4\n",
    "area = len(pieces) + 2   # there are 2 missing pieces in the kaggle dataset\n",
    "perimeter = nb_flats[1] + 2*nb_flats[2]\n",
    "w, h = compute_size(area, perimeter)\n",
    "print(f\"Size of puzzle grid: {w} x {h}\")\n",
    "assert w * h == area\n",
    "assert 2 * (w + h) == perimeter\n",
    "\n",
    "grid_size = (w, h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NB_SAMPLES = 19\n",
    "\n",
    "for piece in tqdm.tqdm(pieces):\n",
    "    for edge in piece.edges:\n",
    "        # compute the distance from the first point, this is not exactly edge.arc_length\n",
    "        edge_contour = sub_contour(edge.normalized_piece_contour, edge.idx0, edge.idx1)\n",
    "        deltas = edge_contour[1:] - edge_contour[:-1]\n",
    "        distances = np.cumsum(np.sqrt(np.sum(deltas**2, axis=2)))\n",
    "        distance = distances[-1] / (NB_SAMPLES - 1)  # distance between 2 sample points\n",
    "        # get N equidistant points\n",
    "        sample_indices = (np.array([np.argmax(distances >= i*distance - 0.0001) for i in range(NB_SAMPLES)]) + edge.idx0) % len(piece.contour)\n",
    "\n",
    "        edge.update(\n",
    "            sample_indices=sample_indices,\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show some problematic pieces\n",
    "pieces.sort(key= lambda piece: piece.id)\n",
    "\n",
    "for piece in [pieces[514], pieces[818], pieces[998], pieces[946]]:  # pieces[:1]:\n",
    "    plt.axis('equal')\n",
    "    plt.title(f\"id {piece.id} {piece.filename}\")\n",
    "    plt.imshow(piece.img)\n",
    "    plot_contour(piece.contour, c=\"yellow\", ls=':')\n",
    "    for edge, marker in zip(piece.edges, 'x^+v'):\n",
    "        plot_contour(piece.contour[edge.sample_indices], marker=marker, ls='', c=\"red\")\n",
    "    plt.axvline(x=0, c=\"gray\", ls=\":\")\n",
    "    plt.axhline(y=0, c=\"gray\", ls=\":\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Start the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = Item(\n",
    "    grid = {},  # key=(i, j), value=piece\n",
    "    grid_size = grid_size,\n",
    ")\n",
    "\n",
    "def plot_solution(solution):\n",
    "    for piece in solution.grid.values():\n",
    "        plot_contour(piece.placed_contour)\n",
    "\n",
    "def img_solution(solution):\n",
    "    width, height = 1000, 1000\n",
    "    for piece in solution.grid.values():\n",
    "        x, y, w, h = cv2.boundingRect(piece.placed_contour.astype(int))\n",
    "        width = max(width, x+w)\n",
    "        height = max(height, y+h)\n",
    "\n",
    "    img = np.zeros((height, width), dtype=np.uint8)\n",
    "    for piece in solution.grid.values():\n",
    "        cv2.fillPoly(img, [piece.placed_contour.astype(int)], 255)\n",
    "    return img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute the border"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_pieces = [piece for piece in pieces if piece.nb_flats > 0]\n",
    "after_flat_features = {}  # key=piece, value=features\n",
    "before_flat_features = {}\n",
    "\n",
    "for piece in flat_pieces:\n",
    "    for edge in piece.edges:\n",
    "        if edge.sign == 0 and edge.prev.sign != 0:\n",
    "            first_flat = edge\n",
    "        if edge.sign == 0 and edge.next.sign != 0:\n",
    "            last_flat = edge\n",
    "    piece.update(\n",
    "        first_flat = first_flat,\n",
    "        last_flat = last_flat,\n",
    "        before_flat = first_flat.prev,\n",
    "        after_flat = last_flat.next,\n",
    "    )\n",
    "    before_flat_features[piece] = Item(\n",
    "        sign=piece.before_flat.sign,\n",
    "        length=piece.before_flat.length,\n",
    "        points=sub_contour(piece.first_flat.normalized_piece_contour, piece.before_flat.idx0, piece.before_flat.idx1),\n",
    "        sample_points=piece.first_flat.normalized_piece_contour[piece.before_flat.sample_indices][::-1],        \n",
    "    )\n",
    "    after_flat_features[piece] = Item(\n",
    "        sign=piece.after_flat.sign,\n",
    "        length=piece.after_flat.length,\n",
    "        points=sub_contour(piece.last_flat.normalized_piece_contour, piece.after_flat.idx0, piece.after_flat.idx1),\n",
    "        sample_points=piece.last_flat.normalized_piece_contour[piece.after_flat.sample_indices],\n",
    "    )\n",
    "\n",
    "\n",
    "def plot_border_pieces(ordered_border):\n",
    "    offset = np.array([0., 0.])\n",
    "    plt.axis('equal')\n",
    "    for idx, piece in enumerate(ordered_border):\n",
    "        contour = (piece.last_flat if idx==0 else piece.first_flat).normalized_piece_contour\n",
    "        plot_contour(contour + offset)\n",
    "        offset += contour[piece.last_flat.idx1, 0]\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "class BorderMatcher():\n",
    "    def __init__(self, max_error=np.inf):\n",
    "        self.max_error = max_error\n",
    "\n",
    "    def eval_features(self, features0, features1):\n",
    "        raise NotImplemented()\n",
    "\n",
    "    @cache\n",
    "    def eval_pieces(self, piece0, piece1):\n",
    "        features0 = after_flat_features[piece0]\n",
    "        features1 = before_flat_features[piece1]\n",
    "        if features0.sign == -features1.sign and abs(features0.length - features1.length) < 25:\n",
    "            return self.eval_features(features0, features1)\n",
    "        else:\n",
    "            return np.inf\n",
    "\n",
    "    def next_matches(self, piece0):\n",
    "        matches = []\n",
    "        for piece1 in flat_pieces:\n",
    "            error = self.eval_pieces(piece0, piece1)\n",
    "            if error < self.max_error:\n",
    "                matches.append((error, piece1))\n",
    "        return sorted(matches, key=lambda result: (result[0], result[1].id))\n",
    "\n",
    "    def next_candidates(self, piece0):\n",
    "        return [piece1 for error, piece1 in self.next_matches(piece0)]\n",
    "\n",
    "    def prev_matches(self, piece1):\n",
    "        matches = []\n",
    "        for piece0 in flat_pieces:\n",
    "            error = self.eval_pieces(piece0, piece1)\n",
    "            if error < self.max_error:\n",
    "                matches.append((error, piece0))\n",
    "        return sorted(matches, key=lambda result: (result[0], result[1].id))\n",
    "\n",
    "    def prev_candidates(self, piece1):\n",
    "        return [piece0 for error, piece0 in self.prev_matches(piece1)]\n",
    "    \n",
    "    def evaluate(self):\n",
    "        mismatches = []  # list of piece\n",
    "        best_matches = set()  # list of (error, piece0, piece1)\n",
    "        for piece in tqdm.tqdm(flat_pieces):\n",
    "            error, prev_piece = self.prev_matches(piece)[0]\n",
    "            best_matches.add((error, prev_piece, piece))\n",
    "            error, next_piece = self.next_matches(piece)[0]\n",
    "            best_matches.add((error, piece, next_piece))\n",
    "            error, prev_next_piece = self.prev_matches(next_piece)[0]\n",
    "            if piece != prev_next_piece:\n",
    "                mismatches.append(piece)\n",
    "        print(f\"{len(mismatches)} mismatches: {[piece.id for piece in mismatches]}\")\n",
    "        best_matches = list(best_matches)\n",
    "        best_matches.sort(key=lambda match: (match[0], match[1].id, match[2].id))\n",
    "        print(\"Show 1 best match and 3 worst matches\")\n",
    "        for error, piece, next_piece in best_matches[:1] + best_matches[-3:]:\n",
    "            plt.title(f\"Candidate {piece.id}-{next_piece.id}, error={error}\")\n",
    "            plot_border_pieces([piece, next_piece])\n",
    "            plt.show()\n",
    "\n",
    "    def plot_candidates(self, piece=None):\n",
    "        if piece is None:\n",
    "            piece = random.sample(flat_pieces, 1)\n",
    "        for error, candidate in self.next_matches(piece):\n",
    "            plt.title(f\"Candidate {piece.id}-{candidate.id}, error={error}\")\n",
    "            plot_border_pieces([piece, candidate])\n",
    "            plt.show()\n",
    "\n",
    "    def plot_histograms(self):\n",
    "        nb_candidates = [len(self.prev_candidates(piece)) for piece in flat_pieces]\n",
    "        plt.title(\"Number of candidates before flat\")\n",
    "        plt.hist(nb_candidates, bins=20)\n",
    "        plt.show()\n",
    "\n",
    "        nb_candidates = [len(self.next_candidates(piece)) for piece in flat_pieces]\n",
    "        plt.title(\"Number of candidates after flat\")\n",
    "        plt.hist(nb_candidates, bins=20)\n",
    "        plt.show()\n",
    "\n",
    "        min_error = [self.prev_matches(piece)[0][0] for piece in flat_pieces]\n",
    "        plt.title(\"Min error before flat\")\n",
    "        plt.hist(min_error, bins=30)\n",
    "        plt.show()\n",
    "\n",
    "        min_error = [self.next_matches(piece)[0][0] for piece in flat_pieces]\n",
    "        plt.title(\"Min error after flat\")\n",
    "        plt.hist(min_error, bins=30)\n",
    "        plt.show()\n",
    "    \n",
    "    def get_issues(self):\n",
    "        issues = Counter()\n",
    "        for piece in flat_pieces:    \n",
    "            best_prev_candidates = self.prev_candidates(piece)\n",
    "            if len(best_prev_candidates) == 0:\n",
    "                issues['no_prev_candidate'] += 1\n",
    "            best_next_candidates = self.next_candidates(piece)\n",
    "            if len(best_next_candidates) == 0:\n",
    "                issues['no_next_candidate'] += 1\n",
    "                continue\n",
    "            best_next_piece = best_next_candidates[0]\n",
    "            best_prev_next_candidates = self.prev_candidates(best_next_piece) + [None]\n",
    "            best_prev_next_piece = best_prev_next_candidates[0]\n",
    "            if piece != best_prev_next_piece:\n",
    "                issues['mismatch_piece'] += 1\n",
    "        return(issues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Place the first corner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "piece0 = [piece for piece in pieces if piece.nb_flats == 2][0]\n",
    "\n",
    "h = int(piece0.last_flat.length)\n",
    "w = int(piece0.first_flat.length)\n",
    "\n",
    "# Binary mask of active edges to match\n",
    "PAD = 50\n",
    "img_mask = np.zeros((600, 600), dtype=np.uint8)\n",
    "cv2.rectangle(img_mask, (0, 0), (2*PAD, h), 255, -1)\n",
    "cv2.rectangle(img_mask, (0, 0), (w, 2*PAD), 255, -1)\n",
    "\n",
    "img_edges = np.zeros_like(img_mask)\n",
    "cv2.rectangle(img_edges, (0, 0), (600, PAD), 255, -1)\n",
    "cv2.rectangle(img_edges, (0, 0), (PAD, 600), 255, -1)\n",
    "# plt.title('mask')\n",
    "# plt.imshow(img_mask)\n",
    "# plt.show()\n",
    "# plt.title('edges')\n",
    "# plt.imshow(img_edges)\n",
    "# plt.show()\n",
    "\n",
    "def eval_first_corner(params, debug=False):\n",
    "    img_piece = np.zeros_like(img_mask)\n",
    "    contour0 = transform_contour(piece0.first_flat.normalized_piece_contour, piece0.first_flat.idx1, params)\n",
    "    cv2.fillPoly(img_piece, [contour0.astype(int)], 255)\n",
    "    img_xor = cv2.bitwise_xor(img_edges, img_piece)\n",
    "    img_and = cv2.bitwise_and(cv2.bitwise_not(img_xor), img_mask)\n",
    "    error = np.sum(img_and>0) / np.sum(img_mask>0)\n",
    "    if debug:\n",
    "        print(error)\n",
    "        plt.title('piece')\n",
    "        plt.imshow(img_piece)\n",
    "        plt.show()\n",
    "        plt.title('xor')\n",
    "        plt.imshow(img_xor)\n",
    "        plt.show()\n",
    "        plt.title('and')\n",
    "        plt.imshow(img_and)\n",
    "        plt.show()\n",
    "    return error\n",
    "\n",
    "res = scipy.optimize.minimize(eval_first_corner, [PAD, PAD, 180], method='Powell', tol=0.1)\n",
    "print(res.fun, res.x)\n",
    "piece0.update(\n",
    "    placed_contour=transform_contour(piece0.first_flat.normalized_piece_contour, piece0.first_flat.idx1, res.x),\n",
    "    top_edge=piece0.first_flat,\n",
    ")\n",
    "\n",
    "solution.grid.clear()\n",
    "solution.grid[(0, 0)] = piece0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img_solution(solution))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Border match via sample points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SampleBorderMatcher(BorderMatcher):\n",
    "    def eval_features(self, features0, features1):\n",
    "        diff = features1.sample_points - features0.sample_points\n",
    "        offset = np.mean(diff, axis=0)\n",
    "        error = np.sum((diff - offset)**2)\n",
    "        return error\n",
    "\n",
    "sampleBorderMatcher = SampleBorderMatcher()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleBorderMatcher.evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Border match via points & sample points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PointBorderMatcher(BorderMatcher):\n",
    "    def eval_features(self, features0, features1):\n",
    "        offset = features0.sample_points[0]\n",
    "        error = 0\n",
    "        for point1 in features1.sample_points[::-1]:\n",
    "            distance = np.min(np.sum((features0.points - point1 - offset)**2, axis=2))  # squared distance between sample point0 and closest point1\n",
    "            error += distance\n",
    "        return error\n",
    "\n",
    "pointBorderMatcher = PointBorderMatcher()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pointBorderMatcher.evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Border match via optimized offset on points & sample points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class OffsetPointBorderMatcher(BorderMatcher):\n",
    "    def eval_offset(self, features0, features1, offset):\n",
    "        error = 0\n",
    "        for point1 in features1.sample_points[::-1]:\n",
    "            distance = np.min(np.sum((features0.points - point1 - offset)**2, axis=2))  # squared distance between sample point0 and closest point1\n",
    "            error += distance\n",
    "        return error\n",
    "\n",
    "    def eval_features(self, features0, features1):\n",
    "        offset = features0.sample_points[0][0]\n",
    "        res = scipy.optimize.minimize(lambda offset: self.eval_offset(features0, features1, offset), offset, method='Powell', tol=1)\n",
    "        return res.fun\n",
    "\n",
    "offsetPointBorderMatcher = OffsetPointBorderMatcher()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offsetPointBorderMatcher.evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Composite border matcher"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CompositeBorderMatcher(BorderMatcher):\n",
    "    def __init__(self):\n",
    "        self.max_error = 1000\n",
    "        self.matchers = [SampleBorderMatcher(10000), PointBorderMatcher(2000), OffsetPointBorderMatcher(700)]\n",
    "\n",
    "    def eval_features(self, features0, features1):\n",
    "        error = np.inf\n",
    "        for matcher in self.matchers:\n",
    "            error = matcher.eval_features(features0, features1)\n",
    "            if error > matcher.max_error:\n",
    "                break\n",
    "        return error\n",
    "        \n",
    "compositeBorderMatcher = CompositeBorderMatcher()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compositeBorderMatcher.evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Match border pieces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for piece in tqdm.tqdm(flat_pieces):\n",
    "    next_candidates = compositeBorderMatcher.next_candidates(piece)\n",
    "    piece.after_flat.update(\n",
    "        piece_candidates = next_candidates,\n",
    "        edge_candidates = [piece.before_flat for piece in next_candidates],\n",
    "    )\n",
    "    prev_candidates = compositeBorderMatcher.prev_candidates(piece)\n",
    "    piece.before_flat.update(\n",
    "        piece_candidates = prev_candidates,\n",
    "        edge_candidates = [piece.after_flat for piece in prev_candidates],\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_group(piece0):\n",
    "    group = [piece0]\n",
    "    piece = piece0\n",
    "    while len(piece.after_flat.piece_candidates) == 1:\n",
    "        piece = piece.after_flat.piece_candidates[0]\n",
    "        group.append(piece)\n",
    "    piece = piece0\n",
    "    while len(piece.before_flat.piece_candidates) == 1:\n",
    "        piece = piece.before_flat.piece_candidates[0]\n",
    "        group.insert(0, piece)\n",
    "    return group\n",
    "\n",
    "def compute_groups():\n",
    "    used_pieces = set()\n",
    "    groups = []  # list of pieces\n",
    "    for piece in flat_pieces:\n",
    "        if piece not in used_pieces:\n",
    "            group = compute_group(piece)\n",
    "            used_pieces.update(group)\n",
    "            groups.append(group)\n",
    "    return groups\n",
    "\n",
    "groups = compute_groups()\n",
    "groups.sort(key=len, reverse=True)\n",
    "print(f\"Computed {len(groups)} groups, best group has {len(groups[0])} pieces.\")\n",
    "print(f\"{sum([len(group)-1 for group in groups])}/{len(flat_pieces)} matched edges along the border.\")\n",
    "for group in groups:\n",
    "    if len(group) == 1: break\n",
    "    print(len(group), ' '.join([('K' if piece.nb_flats==2 else '') + str(piece.id) for piece in group]))\n",
    "    plt.title(f\"Group has {len(group)} pieces and {len([piece for piece in group if piece.nb_flats==2])} corner piece\")\n",
    "    plot_border_pieces(group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "piece0 = solution.grid[(0, 0)]\n",
    "y0, x0 = piece0.placed_contour[piece0.first_flat.idx0][0]\n",
    "\n",
    "# Binary mask of active edges to match\n",
    "img_mask = np.zeros((1000, 1000), dtype=np.uint8)\n",
    "cv2.rectangle(img_mask, (0, int(y0)+PAD), (2*PAD, int(y0)+PAD+200), 255, -1)  # Why y0+PAD ???\n",
    "cv2.polylines(img_mask, [sub_contour(piece0.placed_contour, piece0.after_flat.idx0, piece0.after_flat.idx1)[PAD:-PAD].astype(int)], False, 255, PAD)\n",
    "plt.title('mask')\n",
    "plt.imshow(img_mask)\n",
    "plt.show()\n",
    "\n",
    "img_edges = np.zeros_like(img_mask)\n",
    "cv2.fillPoly(img_edges, [piece0.placed_contour.astype(int)], 255)\n",
    "cv2.rectangle(img_edges, (0, int(y0)+PAD), (PAD, int(y0)+h), 255, -1)\n",
    "plt.title('edges')\n",
    "plt.imshow(img_edges)\n",
    "plt.show()\n",
    "\n",
    "def eval_border_piece(piece, params, debug=False):\n",
    "    x, y, degrees = params\n",
    "    img_piece = np.zeros_like(img_mask)\n",
    "    contour = transform_contour(piece.first_flat.normalized_piece_contour, piece.first_flat.idx0, (x, y), degrees)\n",
    "    cv2.fillPoly(img_piece, [contour.astype(int)], 255)\n",
    "    img_xor = cv2.bitwise_xor(img_edges, img_piece)\n",
    "    img_and = cv2.bitwise_and(cv2.bitwise_not(img_xor), img_mask)\n",
    "    error = np.sum(img_and>0) / np.sum(img_mask>0)\n",
    "    if debug:\n",
    "        print(error)\n",
    "        plt.title('piece')\n",
    "        plt.imshow(img_piece)\n",
    "        plt.show()\n",
    "        plt.title('xor')\n",
    "        plt.imshow(img_xor)\n",
    "        plt.show()\n",
    "        plt.title('and')\n",
    "        plt.imshow(img_and)\n",
    "        plt.show()\n",
    "    return error\n",
    "\n",
    "# sample matching: after 129 [(572, 140), (1534, 1379), (2014, 1862)]  # (score, piece)\n",
    "# xor matching with default params: slow and inaccurate\n",
    "# xor matching with optimized params after 129: slow and very accurate (0.0339, 140), (0.0521, 1379), (0.0667, 1862)  # (score, piece)\n",
    "\n",
    "matches = []\n",
    "for piece1 in candidates_after_flat(piece0):\n",
    "    matches.append(scipy.optimize.minimize(lambda params: eval_border_piece(piece1, params), [PAD, PAD+y0, -90], method='Powell', tol=0.1)\n",
    "    print(piece)\n",
    "# piece1 = pieces[998]\n",
    "# eval_border_piece(piece1, [PAD, PAD+y0, -90], True)\n",
    "# eval_border_piece(pieces[140], [PAD, PAD+y0, -90])\n",
    "# eval_border_piece(pieces[1705], [PAD, PAD+y0, -90])\n",
    "# res = scipy.optimize.minimize(lambda params: eval_border_piece(piece1, params), [PAD, PAD+y0, -90], method='Powell', tol=0.1)\n",
    "# eval_border_piece(piece1, res.x, True)\n",
    "\n",
    "# best_pieces = [(eval_border_piece(piece, [PAD, PAD+y0, -90]), piece.id) for piece in tqdm.tqdm(flat_pieces)]\n",
    "# best_pieces.sort()\n",
    "# print(best_pieces)\n",
    "# print(best_pieces[:10], best_pieces[-1:])\n",
    "# [(0.09865159254807693, 1705), (0.10021033653846154, 874), (0.10109299879807693, 1379), (0.10372220552884616, 295), (0.10691481370192307, 1158)]\n",
    "# res = scipy.optimize.minimize(place_corner, [PAD, PAD, 180], method='Powell')\n",
    "# # print(res)\n",
    "# x, y, degrees = res.x\n",
    "# piece0.update(\n",
    "#     placed_contour=transform_contour(piece0.first_flat.normalized_piece_contour, piece0.first_flat.idx1, (x, y), degrees),\n",
    "# )\n",
    "# \n",
    "# solution.grid.clear()\n",
    "# solution.grid[(0, 0)] = piece0\n",
    "# for piece in solution.grid.values():\n",
    "#     plot_contour(piece.placed_contour)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# matrix = cv2.getRotationMatrix2D(p0, math.degrees(angle_radians), 1)\n",
    "# normalized_piece_contour = cv2.transform(contour, matrix) - p0\n",
    "\n",
    "# img = np.pad(p0.img, PAD)\n",
    "# plt.imshow(img)\n",
    "# plot_contour(sub_contour(p0.contour, e0.idx0, e0.idx1) + (PAD, PAD), c='red')\n",
    "# plt.show()\n",
    "\n",
    "# Thick binary mask following the edge\n",
    "# img_mask = np.zeros_like(img)\n",
    "# cv2.polylines(img_mask, [sub_contour(p0.contour, e0.idx0+10, e0.idx1-10)], False, 255, 20)\n",
    "# plt.imshow(img_mask)\n",
    "# plt.show()\n",
    "\n",
    "# cv2.fillPoly(img, [contour0.astype(int) + PAD], 255)\n",
    "# plt.imshow(img)\n",
    "# plt.show()\n",
    "# \n",
    "# plt.axis('equal')\n",
    "# plot_contour(contour0)\n",
    "# plot_contour(contour0[[flat0.idx1]], ls='', marker='x', c='red')\n",
    "# \n",
    "# point0 = contour0[flat0.idx1]\n",
    "# plot_contour(contour1 + point0)\n",
    "# plt.show()\n",
    "# \n",
    "# # img = np.pad(p1.img, PAD)\n",
    "# \n",
    "# img_xor = cv2.bitwise_xor(img_p1, img_p0)\n",
    "# plt.imshow(img_xor)\n",
    "# plt.show()\n",
    "\n",
    "# cv2.polylines(img_mask, [sub_contour(p0.contour, e0.idx0+10, e0.idx1-10)], False, 255, 20)\n",
    "# plt.imshow(img_mask)\n",
    "# plt.show()\n",
    "# \n",
    "# cv2.bitwise_and(img_mask, img)\n",
    "\n",
    "# \n",
    "# plt.imshow(img)\n",
    "# plot_contour(sub_contour(p1.contour, e1.idx0, e1.idx1) + (PAD, PAD), c='red')\n",
    "# plt.show()\n",
    "\n",
    "# def draw_piece(img, piece, idx, xy, degrees):\n",
    "#     h, w = piece.img.shape()\n",
    "#     piece_img = cv2.copyMakeBorder(piece.img, PAD, PAD, PAD, PAD, cv2.BORDER_CONSTANT)\n",
    "#     matrix = cv2.getRotationMatrix2D((PAD + w/2, PAD + h/2), degrees, 1)  # rotate around the center\n",
    "#     pt = cv2.transform(piece.contour + PAD, matrix)\n",
    "#     cv2.warpAffine(piece_img, matrix, (w+2*PAD, h+2*PAD), img)\n",
    "# \n",
    "# piece0 = pieces[7]\n",
    "# flat0 = piece0.first_flat\n",
    "# contour0 = flat0.normalized_piece_contour\n",
    "# ref_contour = sub_contour(contour0.astype(int), flat0.next.idx0, flat0.next.idx1)\n",
    "# piece1 = pieces[1579]\n",
    "# flat1 = piece1.last_flat\n",
    "# contour1 = flat1.normalized_piece_contour\n",
    "# \n",
    "# # plt.axis('equal')\n",
    "# # plot_contour(contourX, c='red')\n",
    "# # plot_contour(contourY, c='green')\n",
    "# # plot_contour(sub_contour(contourX, flat0.next.idx0, flat0.next.idx1), marker='*')\n",
    "# # plt.show()\n",
    "# \n",
    "# # Binary mask of the edge\n",
    "# pad = 20  # add padding for visualization\n",
    "# x, y, w, h = cv2.boundingRect(ref_contour)\n",
    "# img0 = np.zeros((h+2*pad, w+2*pad), dtype=np.uint8)\n",
    "# cv2.fillPoly(img0, [contour0.astype(int) + (-x + pad, -y + pad)], 255)\n",
    "# imgMask = np.zeros_like(img0)\n",
    "# cv2.polylines(imgMask, [ref_contour[10:-10] + (-x + pad, -y + pad)], False, 255, 10)\n",
    "# img1 = np.zeros_like(img0)\n",
    "# # x, y, w, h = cv2.boundingRect(sub_contour(contourY.astype(int), flat1.prev.idx0, flat1.prev.idx1))\n",
    "# # contourY = transform_contour(contour1, flat1.idx0, (0, 0), 0)\n",
    "# # cv2.fillPoly(img1, [contourY.astype(int) + (-x + 1 + pad, -y + -4 + pad)], 255)\n",
    "# # imgXor = img0.copy()\n",
    "# # cv2.bitwise_xor(img0, img1, imgXor)\n",
    "# # imgAnd = img0.copy()\n",
    "# # cv2.bitwise_and(cv2.bitwise_not(imgXor), imgMask, imgAnd)\n",
    "# # print(f\"{np.sum(imgAnd>0)}/{np.sum(imgMask>0)} = {np.sum(imgAnd>0)/np.sum(imgMask>0)}\")\n",
    "# \n",
    "# plt.imshow(img0)\n",
    "# plt.show()\n",
    "# plt.imshow(imgMask)\n",
    "# plt.show()\n",
    "# plt.imshow(img1)\n",
    "# plt.show()\n",
    "# plt.imshow(imgXor)\n",
    "# plt.show()\n",
    "# plt.imshow(imgAnd)\n",
    "# plt.show()"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}