c_registration_overlap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4ad664ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import generate_square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ac42e978",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAIABJREFUeJztnW2obNd533//tfbMOfdFinz9VkUSWKVKQQ4lCUIuuBQ3Tms3CbEpuCiQooJAX1ya0EIsNdCSDwa3H0L6oYGKJFQlL6pwEixM2lRWGkIgiWwldmPZcXwTu/ZFwmoqy7pvZ2bvtZ5+WGvv2TPnXN0t3XvumXP1/GAzs9fsPbPm3rP+87yttWRmOI7jXI1w1B1wHOd44GLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTOLQxELSByV9RdJZSY8c1uc4jnNj0GHUWUiKwF8A/xA4B3wW+HEz+9J1/zDHcW4Ih2VZ3A+cNbO/MrMl8ATwoUP6LMdxbgDNIb3vHcA3R+fngPdc6eK5dmyXU4fUFcdxAM7z7b82s7e/0fsPSyx0QNuavyPpYeBhgF1O8h69/5C64jgOwGfsk//nWu4/LDfkHHDX6PxO4IXxBWb2mJndZ2b3zdg5pG44jnO9OCyx+Cxwj6S7Jc2BB4CnDumzHMe5ARyKG2JmnaR/AfwOEIFfNrPnD+OzHMe5MRxWzAIz+23gtw/r/R3HubF4BafjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTOKqYiHplyW9JOmLo7Yzkp6W9NX6+JbRa49KOivpK5I+cFgddxznxjLFsvgvwAc32h4BnjGze4Bn6jmS7gUeAN5d7/kFSfG69dZxnCPjqmJhZr8PvLzR/CHg8fr8ceDDo/YnzGxhZl8DzgL3X6e+Oo5zhLzRmMU7zexFgPr4jtp+B/DN0XXnats+JD0s6XOSPteyeIPdcBznRnG9A5w6oM0OutDMHjOz+8zsvhk717kbjuNcb96oWHxL0u0A9fGl2n4OuGt03Z3AC2+8e47jbAtvVCyeAh6szx8EPjVqf0DSjqS7gXuAZ6+ti47jbAPN1S6Q9OvA+4C3SToH/DvgE8CTkh4CvgF8BMDMnpf0JPAloAM+ambpkPruOM4N5KpiYWY/foWX3n+F6z8OfPxaOuU4zvbhFZyO40zCxcJxnEm4WDiOMwkXC8dxJuFi4TjOJFwsHMeZhIuF4ziTcLFwHGcSLhaO40zCxcJxnEm4WDiOMwkXC8dxJuFi4TjOJFwsHMeZxFWnqDvOgA5aNfEasANXXHS2FBeLNxv9gO8H6vUWAOemxd0Qx3Em4ZbFm4Gx9aD+9yG/xvUH/IbYa1zvvClwsTiuTHEfRoNeYf16y1d+7SDG10/CxeWmw8XiOHOQBbD2slbXbAiCNAou9q9lO/gcUFgNfssTApMKLhg3GS4Wx40Ni2JsFfSDeFMk1N8TRuKi0UDu2/u20bmZIYnxGu2Tdq+1DMTSJxeNmwIXi+PAgTGHDVFgNIh7gQih3FsHv0bvY6O0Zd/etw3nKQ1bzGksNDmv3b+PbOuWRd9HF41jjYvFceFKLscB8YZNoVAM5f7+WgmZ7bdSNt/HmjLwgyCtb/+ibCvByBsiYFYFRSt3Z2xpDOfOccLFYhs5IHi5FoQcuxjxAJ8gxioY1fLYEAtJB1oWwOB29M+HuMW4T2YgQ3XAD+/Ui4ZGnxPYb2kM3yF7YdYxwsXiuLEZizhIWMZCUR+LdbGKXWj9htUAz3k419hiyKuBbWbV0ghghnohiLVv2YqAmJX3iL04xTVLw3I4eCvtnrG4uKgcOS4WNwOjeIIkiHHNmlDTwKxZicXYkghaEwt690TrrofSauAqZ+jPB+FYt0TUdUUQUoaci9tTbl4XDnudWRbnyHCxOI6ElYVACCtLAopANM3KkogRmojNmpWohA3RkLBYLYsuYyFAAHVhZU2ElWWhVGMiY7ejtyJUYxlBQCyGQ+/2bMY2nGOFi8VxYyPTsRa87F2JWASjiEXAmiIYhIANrsjoPUeWhQUglnMLrCwCGKwGo1oaIawskbW3U/mcA1yHVVDUU6rHDReLbUPr6VCowc0ayFTTrKyJGPdbEWOBCAFmDRaFzRssCGvC8DkWVsLRP1c2rL4W2oyyVath5GKkjFK1FDKo7YYMiPpr2q68d9dhXUI5rYKlZtUdGRV/WR6CuO6ObCcuFseIPksxuB1BK6HoU6W9UFRxsfqaNYE8i1hTf/XF4HpYEP24VS6WggmsUREFK+29laHOVoJhBk1xSdSN6i+CynsFoRCwTpANpbSq56huScmYjOtHanvfZnn/bFnnhuNiseX0hVcHV2FqZVHEADFiMcB8hs1isShmkdwE8k4kzQPWVIGIK5GwwCgbYoNwhFTFw0CJNbEIrQ1WRGhzEY82FUvDrAhEzmhRu9rHM2IsgtFbQX0QNWitvHzjH8Fdli3AxWILueLErrH70Wc5mqYIRQgliBkDNoukEzMIIs8DOYq0G0k7Is9WAmGxxCVsrQpco+egZKgXi2woQ0iGNTa0IVBXYikyg1QHfSdocglyxhogrelVpVQCqTUwOgRF1+ajuFuyTbhYbDt98DKEakHEel6fz5piTcyaYk00gbzTkE405KaIRZqL7oTodkSeAYLcqFgX9Rgqq/oMZy8QSeV5Z4QkyBA6iEsjdEZIENqAOiMui2uiLhNhqNUwRunWURGZrLw2VJTmvO6SVGtCYX1uinM0uFgcA9aKr6pQqBeIJg7WhM2b4nLMI2knYFF0uypisSu6k0UsLEBuelcErDF636OfS1asinpuEFoN56G1EgDtREhGaETorKRbU722Jkksx1Vp+Ti9G6yKgrsXx4WrioWku4D/CvwNyv/sY2b2HyWdAf4b8C7g68A/NbNv13seBR4CEvAvzex3DqX3bwZGLonGghFrtqOJWB/AnEfyLJB2ypFnVSx2oDsh2lOQ5xTLYmbVFQGLBpTYxCAag2VRBWI5FotilYTWUBKxhdCChUDorFg6QAhFNErhVlODoKlYEB2gjCwU6yLnVWVpVawhW5LSKEPkJeJHxRTLogP+tZn9iaRbgOckPQ38c+AZM/uEpEeAR4CPSboXeAB4N/DdwGckfY+ZG5JXRaNf3qGpCIMGFySWGaAhYDEUoZjF4nrsliBmOhFoT2hwP/IOtCehO2XFDQlGnhn0QjEf/bonrT1XUrEs9lTiElUs4qI8hg7SAmILqS0WRlyWLExoBAFCFCEE6GpqNaXyfbquzkUJ0GdJDircUqD87jhHyVXFwsxeBF6sz89L+jJwB/Ah4H31sseB3wM+VtufMLMF8DVJZ4H7gT+83p1/U7BZddlbFbFYFTQBm5UjzQN55HKkuehOQpoXoehOGTbPxQWYGYqZ2BjNrBs+LqfyedlE7gKWBCZSDGgZCC3F3chaxTlKoqPYnSYUDWvAsoaaDZNKUiMI9RFVi7VsPBWLJrnFsM28rpiFpHcB3w/8MfDOKiSY2YuS3lEvuwP4o9Ft52qbM4V+PkdvZcQ4qsgMNQMSsSZiMRarYh5JO5E8rwKxUyyKNIe0C2m3CsWpjrCbCMGYzTuaJjGLiRMjsVimSMqiS5G2i+QschZdbLAYSTOR2zLYLQhrKZPB+loNGTKR2yImeR6QGdbVLEgoboegiFYoYlTSwzZkSIBVlsTZCiaLhaTTwG8AP2Vmr46nNW9eekDbvv9xSQ8DDwPscnJqN94UHFhTAauS7v7o059NSY/mph69SMygO2HkHcNOJOLJjvm8I8bMvOmYN4ndpuPUbEmQETD2UkOywOV2xqJrWHaRbGVAt1lAwMjkZSwWRBY5ltgGNQ1bLz/4L2H8d1NF0bQR6PTtCbaSSWIhaUYRil81s9+szd+SdHu1Km4HXqrt54C7RrffCbyw+Z5m9hjwGMCtOuM/H5uM537UeIVCWFVozmcl+zGLJVZxIgzp0faUaE8beQ7dqYzNM7NbF3zX6T1u3d3jRNOyG1t2Y8etsz1uafaGj13khi5HLqY5F7s5l7o5bYq8urvDxcWcxWJG10ZSAJsFwrJmQ2vVZ6xp0FwzIrkR6soUeeVaZ2FlCvtQM7JW9Tn6ZQmhpl5tVHMR8Dzq0TAlGyLgl4Avm9nPjV56CngQ+ER9/NSo/dck/RwlwHkP8Oz17PSbDu1fJm+wKOaB3Ig0q8dOyX4U9yNjJxNhJ3HqxJLbTlzmrbsXuaVZcKpZcDIsecvsIm9vzpNqHGHPZizyjPNpl293J7nY7XA5zZjFxCxkLsbE5cWcvSQyTbEKsoYUqwxSEmFWzsO8BElDG8iAmprxSMKCoTbURXJUUqkhr6rE8miy2lAOnuhnuTo3limWxXuBfwb8maTP17Z/QxGJJyU9BHwD+AiAmT0v6UngS5RMykc9EzKBAyaQXY2hXLuvmWhWz/PMsJmhWSY2iZ1Zx4mm5ZZmwW2zS3xXc5nTcY+3N+d5e3yVtpRRsZfn7NmM03GP3dByPu5yIe0Mnxnq4jVdF2hTiT9YF8i1gMsSxf3oq0RHLlNp789Lgdeaa1WTHsNKXdf6b+pcV6ZkQ/6Ag71PgPdf4Z6PAx+/hn45mwwDqgRArc4DKSXbGqyLPKMUXs0MZplm3jGfJ07Nl9w2v8SZ+UXeNjvPmXiRW+Jl3hovcFu4TGtVLLRkz2bM1JFNBGWiMl2OdBbIJtocuNxkuiZjGfJMhFbkBtRVsWoAExatxFFiIAA2qxWcSUh5VZmac6mniJEhTWp1gtp4eT7nyPAKzm1gFNAb5oX0JrjCfotjtCaFxV4oRhZGrHM3Zpn5PHF6d8Fbdy9y++6r3D5/hbtmL3NbvMht4TJnwpIzMbKwlgBcsstczIHv5B121fJyOs2uSrakre5BmyMXmkTbZFJWqdcIYMGKePVl5KlaP7EWgFkRtuJS1arRzZXL3QjdWlwstpXRlPNVkLMUYpX2Wr8QSzaCWo2ZqztCNEJjzGKiCZnd2HIyLDkZlpwKC24Je5xUxy1BnNScSFs/OEHItLScsgWZPpbRcDru0OXIXjPjxLxdpVabTJ4JZZUlLiLkKBSMHMtaOgStTxyVKFWjB1gLU7ZYdG44r89Jdm4sm6nT6orY6HE8Y7SfQVomZxmSEUIe4gwZ0Vok1cc6H5RcB2UY/TnMyOyqHY6ZEjuhYyd2NKEIUAiGgg3xiL7Ooo9ZrGIXWqVSa4HW2neC0SLD+z3etVXGnSPDxWKLGFyQ3qoYXhg979e6rIvgqq+etNUksHKdMCDnEmfocqTNkUWeDUHM1kqGYiwSEQ1nESNUIYkYiUCy/YOZWr0pG1V2bhx9HxktonNF+oV9+hiNsxX4/8S20Jvem/NCFNbFY1yUBasBOdzUH1YKreqqUxmVwV6tiTgh15AQLZFltUZ6chUiG008W/sqm289LtI6aK2OsXXxOjNCzo3DYxZbjJmVIdrv5VEaV49GsS6GX21GloZKQLEO5oARyUSMmbrqWmRmiNlo89JMZo4xV2amxIzEXInY319FaG3I9+4GNchZl+0b2l+jINOChhXAh+82nl3r80W2BheLbcPyKhOSczm31XyJgbp8HTYajHVg9jGLniZkmpAGyyIT2LMZu9axZy0La2ktkTAWlrmYjYvWcCnvcMl2uJh32LOGRW64nGbspRmX24a2jWWeSCvCUoS2zkRNlEVwRmtijN2QInTGsBAwrH+/8QrgztbgYrGF2LCTVymFHs5TRiGXtS1zrgvnWlmAJpX1JkIqK1nRBVIb2FvOiMG41M252O2wG1peTbsEMkGZk9axm1tajGTGnsEli5zPc17Nu7yadotopB0u5zmL1NCmSM4BywHLQoNC9V/ggO8U6r6nG0FMC6rrjB5gWThbhYvFtrMxcIbl+7VakXuIW2ya/NUNyQZdjixyw16ecSnvMFdxL/bCHheto63jdM8CexZrALSpR8mcJBOZetS4xdBNRtbDKKAp6y2g6jaNthAg5dUWA4N75YvbbCsuFtuMjQeRrZvndY+O3rJQH6/Ipey6lF6L1AXaFNlLDRe7ObOQeDmeJlVXBGAvXiBbcVH2bMaezXglneJ82uXl7jTf7k7yarfL/907zfl2l1cun+Dy5TnpUoMWkXBZhEXvhtT1Odti4ZQ+Quj6vUfqXiM1o0OqK2eNv2e/mZHvYrZVuFhsC3WjYNUVri30A6qsLGUKiA51JY+hLqAuD2tihs5WC+wmysraXRWLtuFSO2c3lkrM25pLRGWyBU6GBYFMSyRZ4FLeYc9mnE8n+E46wXe6E7y8PMXFNOd8u8uF5ZzLyzLzlC6gZV2Ls6WKQhWJ0Wrgw8ZDZqNYi62Eof+ezlbjYrHtWKmEWJ2vD6o+yBk2BmroyuxPS4GUynyOvdQQlLmQdoZCrfPpBNmKlZHqYz/r9ELa4ZX2JN9pd9lLMy4s51xazFkuGmwZ0CKUwOayfmZb5oaoK5ZE6ctoy4BUYy9VJLQpEL43yFbjYrEtKFx5v5C161Yp1HExVjHxNQhGWVQ3kKPRqeHC3g5m4lKcE2S80pzgdLPk/GyXk2E5xCXaGtt4tTvB5TTj1XaX88sdlimu1rNYRmjDYE2oq5/bC1bqYxW1UGxI88I4E7J29Iz3QM0ev9gmXCy2HEsZ0YECRoO6VFeXSqhNhGXAgogL0SwMixCbEkMwUWamCvYuzzETs5iIIXOpmXMhtixSw07sSvDSAp0FLqcZl7o5l7uVNdGlwGJvTlpEWATipUDcKynTuIC46K2Lup9Ia6sUaj2Ky0FZd7OPR9Rsjwc2tx8Xiy2i33lLIUMOq7iFVBeFGQ+yfjcvG1yRoa5hfHTCmhK7WC4jNoPL7QwzkXKgCZlFbkpJuAW6HFjmhovtnL2u4eJizt7ejJxDEYplv3DvuLaiuh69y1FTuL1IlABs3R+1txauFK8wj19sKy4W24rltc2CS1sJdpbl5mr9RbKyW1iNEYQhXiCsA2vLzFRbBpIimLgUc50vEoh1RmpvVXQ5sEgNe13Dom1o24bUxbLKd6vB5RjmgqylR0cFWKw/d44/LhbbzlCUBWhUoJUzWnRFOJpAaDNxWRaZsVhcg1L9YCgHugA5iRyNbGLZNMSYWXQNTUyjjxNtCnQpslg2tIuGvKiZj7amSDsRe6tiuTqUIbarTZN7IVvtup6hS2WjoX67wjy2mMqxnkZ1S2NbcLHYBqyf2FGX2M+2L9g5zBMZ31MHk2pxk7oyWHPTiwVgJZVqIZB2BNFKfGIeCI1hJmLMo7cVKQVSF+jaiC1DcT06ERZhVdJdBSIu+31Pi1j0ls2qBqT2r49brKVN8+q7jAObzlbiYrHtDDGLsPLnc6YsT1V/vfPKFcmp7AoW6q5gfUVnaMr0cQsq8cUkUjD2kgjR6LcjtUzZHCgLFrGIxF6/6G7ZxjB0K6EIbV+AVXdY7y2JbKvnXYYur7Ib4/UpNqyGtUKsQUxcRLYBF4ttxjJmqx3GV6M5rE/ISsXMD8tM2WY0YLFPqQpbAqgsuydKBqUuhZd2IznUikkYVuouaVChdrUpchGH6m4siiD1O6rHZQm0hqXR7KWVgC26QSxUraA+GzK4IO5mHAtcLI4rNRNiNs6GVNO/BjjLZsWG5TLwM/0cDtUl+Kyqx2oi2FC70e+G3pV9TEtmpReLIhSxLe5HbI24LBt+hGVGbe8W5TLFvE58G1Km41iEjdyP3moapVB9R7LtwcXiOGNWBmCXCV3Gogix7McRl1YXnCmrbvd1ELK6+rfKLudlo2RWM0VtUzDqUVOh4wrRUvxVA5pdKRIrQlViFHQZpVSep5EbMrgjeQhqlq9zgEvCKqVcTw77X9W5Ai4Wx4VRhmBlTVit5CwmvVIYYhihNUJ1NZQEKovnQq3uFMW6qBsdj8UCajp0JAir2oliWcQapwitFfdnWa2JNqM2DZPFinWRh0AsqWZAchUQG00iG4sJXNk9cWvjSHCx2HZyXzMNltJq7Yd6qC7cK4lQJ2kpRVIysEhIxRUJSeRYF+id2Uosmn7KOyOhqPNNqhURWxsyHf3ksN6aCKmIRVimck+bikBYH9hMK6Ho60T6zMdawJZVOtXckthGXCyOI2Proq+5yBlymWkqidAI64zVzmCD5gBl2bvczwBde+/VYxhXZeb1WaTD0Zdz9+5EqgHXzdWwepdpXNp9gOWwtpJ3LxRjwXCr4shwsThO9ANJ/cAs+34qZKzthjqMACV+EepqOCoZjtz0lkfdSqBuHQAMKVatiUW1IoYybhsmhoWuL7oq8RItUhGILqO2W6sBYRTctK5bdz/6NS04OJhpninZGlwston+F1RhbZ5Iv+iuqomuasJbokwyo4iAwSAYkRqQrK6HzcKwX0fTaAh+DvuOjFfd6pvqYrlDkVU/zdwYCq2GIGZbq0C7VMQCavC1GywhUirneXTOSCQ2LQpnq3CxOA70gyj2voPKQKxbBEpl0I3jF8XiEBCKEHR1QyBK0dWw+U8vJmPLon5eqGKxWoti5W4MqdFxENNsVUcB67NKR0vneTr0eOJisS30JZSw9staVs/qA5ygEFbWRaIIRq7L6ceAujTEM0J9T5OgKVPVZcVF6ZfrH6/p2fdD1sc4qlh0uYhIl0dxElAfrMx5vzUBJeuR+nkgJVZhvYsC+4OZ+767i8o24WJx3KiWhUFZLbsGCy3EIhiNIatL7/UCFMIgRjJDbRGM8WZFa5ZFP4AHy2KU3YC19SiGmERX3ZCUsMGySEP2YxyfuJLbsU8chgCni8Y24GJxDNiMX6j6+6aS6hjCDR1DFqIEGxMW6+uhbrBcUVjFMMo5rC3Hn+HAtTJhlQbt6yVyLoFLKOLQi0VKWI1LXNEF2YxP7Dt3odgWXCy2nbXBEzdes1UMo65hN4iJjeINgMlWg7gy7ARWLoDxrubjOEMvQONy7Y3Vrla1E7ayFMbvBb5S9zHHxeK4MSrS6jcNLovgjNakYCQEsYrGWkzEsM2FdXqXZDSg18RhXDy1GbjMaeVSbFgTNurX0P+hoxvuh1sVW42LxTZx0OAY796VDY2Ni9E2hxZCEQyrJdUpQSxTzIGyu9k4iCnte/+yt+hGAdRYGMYDfxy07OMWsB6TqOdrX/GgoKUXXR0LXCxuZvr9R6S1QWu9WKwJRaiL6R5Q+9BbDptWQn/daKboWgVmf83UiWAuFFvNVfe3l7Qr6VlJX5D0vKSfre1nJD0t6av18S2jex6VdFbSVyR94DC/wJsNq9mG1ZFKcLFtsa7D2q6cjw7rurV7+oyF5dVRApO1vetGBVSruRxWLYzxsTb564DFdz39efMwxbJYAD9oZhckzYA/kPTfgX8CPGNmn5D0CPAI8DFJ9wIPAO8Gvhv4jKTvMbP9P0vOG2NU6Uk2CGJYdq+fjEWfWh1ZDwqr+xWK29G/Zf9kbd+OsbWwKsvuYyXlmo1g5rh/B5279XBsuapYWPkruFBPZ/Uw4EPA+2r748DvAR+r7U+Y2QL4mqSzwP3AH17PjjusBn0VjCFuMIpjDLUSEqWKq3dDrD9df8vxFPHSsBKEPp4xHvBTayWcY8+kmIWkCDwH/C3gP5nZH0t6p5m9CGBmL0p6R738DuCPRrefq22b7/kw8DDALiff+De42Vn7JV5VcgLVUqil3iZssBxqTGJ0r8F6NmRjd3ZgEBnHOYhJYlFdiO+TdBvwW5K+9zUuP2gPvn22p5k9BjwGcKvOuG16IxitF6GDxGLTqhiabd09Ocgi8djETc9VA5xjzOwVirvxQeBbkm4HqI8v1cvOAXeNbrsTeOGae+q8foag4+i42vWvE8u2P9sxPpybhinZkLdXiwJJJ4AfAv4ceAp4sF72IPCp+vwp4AFJO5LuBu4Bnr3eHXcOYLT61JVmdvZWwoGZjfr6Qe2OM8UNuR14vMYtAvCkmX1a0h8CT0p6CPgG8BEAM3te0pPAl4AO+KhnQq4T/cAdz04dbUzUM94rdR+9O7FZwQnYfm9x9OJV0qBuRdz0aBt+OW7VGXuP3n/U3Th+DFWY6wO/rNP5GkbjlQa2wmsO+qtWXx5409H/fTmFz9gnnzOz+97o/V7BebMysjr2tb/GPa/benAxeNPgYnEzsDGI+0li/aI5+y9/AwO8n1HqvGl5XdkQZ8u40uC9UoGU41wDLhaO40zCxcJxnEl4zOK4c6U1MGqA83W7Ip4Cda6Ai8XNSL8qlg985zribojjOJNwy+JmxdOcznXGLQvHcSbhYuE4ziRcLBzHmYSLheM4k3CxcBxnEi4WjuNMwsXCcZxJuFg4jjMJFwvHcSbhYuE4ziRcLBzHmYSLheM4k3CxcBxnEi4WjuNMwsXCcZxJuFg4jjMJFwvHcSbhYuE4ziRcLBzHmYSLheM4k3CxcBxnEi4WjuNMwsXCcZxJuFg4jjMJFwvHcSYxWSwkRUl/KunT9fyMpKclfbU+vmV07aOSzkr6iqQPHEbHHce5sbwey+IngS+Pzh8BnjGze4Bn6jmS7gUeAN4NfBD4BUnx+nTXcZyjYpJYSLoT+BHgF0fNHwIer88fBz48an/CzBZm9jXgLHD/9emu4zhHxVTL4ueBnwbyqO2dZvYiQH18R22/A/jm6Lpztc1xnGPMVcVC0o8CL5nZcxPfUwe07dvSW9LDkj4n6XMti4lv7TjOUdFMuOa9wI9J+mFgF7hV0q8A35J0u5m9KOl24KV6/TngrtH9dwIvbL6pmT0GPAZwq87sExPHcbaLq1oWZvaomd1pZu+iBC5/18x+AngKeLBe9iDwqfr8KeABSTuS7gbuAZ697j13HOeGMsWyuBKfAJ6U9BDwDeAjAGb2vKQngS8BHfBRM0vX3FPHcY4UmR29B3Crzth79P6j7obj3NR8xj75nJnd90bv9wpOx3Em4WKjUSwbAAAE5UlEQVThOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zCRcLx3Em4WLhOM4kXCwcx5nEJLGQ9HVJfybp85I+V9vOSHpa0lfr41tG1z8q6aykr0j6wGF13nGcG8frsSz+gZl9n5ndV88fAZ4xs3uAZ+o5ku4FHgDeDXwQ+AVJ8Tr22XGcI+Ba3JAPAY/X548DHx61P2FmCzP7GnAWuP8aPsdxnC2gmXidAf9TkgH/2cweA95pZi8CmNmLkt5Rr70D+KPRvedq2xqSHgYerqcXPmOf/H/AX7+B73AUvI3j01c4Xv31vh4ef/tabp4qFu81sxeqIDwt6c9f41od0Gb7GorgPDbcJH1u5OJsNcepr3C8+ut9PTz6eOMbZZIbYmYv1MeXgN+iuBXfknR77cTtwEv18nPAXaPb7wReuJZOOo5z9FxVLCSdknRL/xz4R8AXgaeAB+tlDwKfqs+fAh6QtCPpbuAe4Nnr3XHHcW4sU9yQdwK/Jam//tfM7H9I+izwpKSHgG8AHwEws+clPQl8CeiAj5pZmvA5j139kq3hOPUVjld/va+HxzX1V2b7wgmO4zj78ApOx3EmceRiIemDtdLzrKRHjro/AJJ+WdJLkr44atvKilVJd0n6X5K+LOl5ST+5rf2VtCvpWUlfqH392W3t6+jzo6Q/lfTpY9DXw620NrMjO4AI/CXwN4E58AXg3qPsU+3X3wd+APjiqO0/AI/U548A/74+v7f2ewe4u36feAP7ejvwA/X5LcBf1D5tXX8pafXT9fkM+GPg725jX0d9/lfArwGf3ua/g9qHrwNv22i7bv09asvifuCsmf2VmS2BJygVoEeKmf0+8PJG81ZWrJrZi2b2J/X5eeDLlCK4reuvFS7U01k9bBv7CiDpTuBHgF8cNW9lX1+D69bfoxaLO4Bvjs4PrPbcEtYqVoFxxepWfAdJ7wK+n/KLvZX9rWb95yl1OU+b2db2Ffh54KeBPGrb1r7CqtL6uVohDdexv1MrOA+LSdWeW85WfAdJp4HfAH7KzF6tqe4DLz2g7Yb110oa/fsk3UZJyX/va1x+ZH2V9KPAS2b2nKT3TbnlgLYb/Xdw3Sutxxy1ZXGcqj23tmJV0owiFL9qZr9Zm7e2vwBm9grwe5SZydvY1/cCPybp6xT3+Acl/cqW9hU4/ErroxaLzwL3SLpb0pwytf2pI+7TldjKilUVE+KXgC+b2c9tc38lvb1aFEg6AfwQ8Ofb2Fcze9TM7jSzd1H+Ln/XzH5iG/sKN6jS+kZGa68Qwf1hSgT/L4GfOer+1D79OvAi0FIU+CHgrZR1O75aH8+Mrv+Z2v+vAP/4Bvf171HMx/8NfL4eP7yN/QX+DvCnta9fBP5tbd+6vm70+32ssiFb2VdKRvEL9Xi+H0vXs79ewek4ziSO2g1xHOeY4GLhOM4kXCwcx5mEi4XjOJNwsXAcZxIuFo7jTMLFwnGcSbhYOI4zif8PC2fb1/NHzZcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initialize 1000 squares with center (x,y) and size (x,y) following a gaussian distribution\n",
    "\n",
    "np.random.seed(1234)\n",
    "shape = (500, 500)\n",
    "arr_smooth = np.zeros(shape)\n",
    "center = np.random.normal(200, 20, size=2000).reshape((1000, 2))\n",
    "size = np.random.normal(50, 20, size=2000).reshape((1000, 2))\n",
    "\n",
    "for i in range(1000):\n",
    "    arr_smooth += generate_square(shape, center[i], size[i])\n",
    "\n",
    "plt.imshow(arr_smooth)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "80324944",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAHFJJREFUeJztnX+MHOd53z8Pf5i8syXUZ56U05GGVJguIhmtbQh0ABeFaqW16hiR/5GhFClUgAD/UeEEbRGRDdAifwhQWkBwgTZFicQog0SRiSaBBCOIK6sVhAK2aUm2E1OKIiZSpZMO4lFyYBonMjry6R87czu7Oz/e2ZnZmdn9foDD7c7N7D63O/Od532e531ec3eEEKKIPW0bIIToBxILIUQQEgshRBASCyFEEBILIUQQEgshRBCNiYWZ3WNmL5vZBTM72dT7CCFmgzVRZ2Fme4G/BP4JsAF8D/gld3+x9jcTQsyEpjyLY8AFd/9rd/9b4HHg3obeSwgxA/Y19LrrwBuJ5xvAZ7J2/oAd8IN8sCFThBAAl/nxJXdfnfb4psTCUraNjHfM7ARwAuAgy3zG7m7IFNF19h053Or777yx0er7z4pv+f/8f1WOb2oYsgEcSTw/DLyV3MHdT7v7ne5+534ONGSGEKIumvIsvgccNbPbgDeB+4F/3tB7AcO70zR3ibbvbOMsyp1O1EPw+ft6xfepdng67r5jZv8K+CawF/iau59v4r2EELOhKc8Cd/8T4E+aen0hxGxRBacQIojGPIt5YGd9pZX3TX4pil+IriCxaIDttaVKxy8zFKmiL2haQdv35rtTHReCBG4+kVgIUUAb2bIuCq5iFkKIIORZTEHVYYYYpYt3UTGJPAshRBASCyFEEBILIUQQEgshRBASCyFEEL3IhoTkuePipGn+oZ31lUaLlEQ16qxzUOZleuRZCCGC6IVn0RSxN7K9tjRSYg2Dcugsb2N837qRlyO6SKfFYrehTcD8h+21JZY332vapEaJi73KiNG0BWLxe7QlTGWGFnVO6CtzwmvIMkqnxaJvVK3s3F6NR4XhrzM8Jp3lresVLBJiyNzELJY335P7LkSDzI1YCCGaRcMQMRV5cYS0k6rs+D/5+nVO3BvpFSJPtBTyLMTM6FoXdVEOiYUQIggNQxpg+hTuUunjlzfz/17kak+bHtSJs3h0+juPT+ROGxlA6Jh7e3WPUp2is/T9OhQCKK43GTIU7qLit31UL8wqKijrU5C102JRpoJzGqaNsheVgY//Pbwis74q1Lxq0D6doKI7dFosFoGkYA3ujvWkCfOqQW98s/rrz2Mp9PjNY5mVgXfR4Poxaa/dVTGXWHSI5a3rNc5vyQ6WzuOFLppHYiFmhkSqGWb1uarOQggRhDyLiLTx6rSvscxKSixilDhFOvm3/JhF36fhd5nxz3bfm++y88aGLpIIfQ5iaqqUbxcFDfPEVrUo7SCxaJDw3H/+Mbo4RBdYaLGQS1+NKoG1ohNvdBjY7nKRCswOUIBTCBGExEIIEUThMMTMvgZ8Ebjo7p+Itq0AXwduBV4DvuzuP47+dgo4DlwDvuLu32zE8o4zGOK05z7HQ6y0asCdWRsj5oKQmMX/AP4L8LuJbSeBp939ETM7GT1/yMxuB+4H7gBuAb5lZh9392uVjAwsf61SlpuMviuWIcQkhcMQd38WGL9a7wXORI/PAF9KbH/c3a+6+6vABeBYTbYKIVpk2mzIze6+CeDum2Z2U7R9HfhOYr+NaNsEZnYCOAFwkOWpjBj3JKr0ahxNWZZvQpPcP/aEwgq7im1O2tHVSUZtE5peLvNZKgsySt0BTkvZ5mk7uvtpd7/T3e/cz4GazRBC1M20nsXbZrYWeRVrwMVo+wZwJLHfYeCtKgaKcoTEbab90nWnXWym9SyeBB6IHj8APJHYfr+ZHTCz24CjwLlqJgohukBI6vQPgLuAQ2a2AfwH4BHgrJkdB14H7gNw9/NmdhZ4kUGG7sGqmZA0kgsaJylbXr28dX33mCuHkiOo/GUEi2IZWfalob6boi8UioW7/1LGn+7O2P9h4OEqRgmRpO+p7M6sl/J6tcMXem6ImA1pF0uZWaeh9F1Uuo7EosOk9dEcT8d2PZW678jhVGEoEoNpZuyWTUOLcmhuiBAiCHkWHSN5Rx0GXScLxmLGPY2Jjl8t3EmTw46d9ZVUL6LIcxgNOIcScu8b2lK2y/m0sYcmu4PD7LxLeRZCiCAkFkKIIDQMEaIjTDu3KWslvLqRWAhBxspgAfuEUCQCZeNKWR3hCycuVqyz0DBECBGEPAsxt6RnVIb3x+WEp5B29y/KNIUSewJ9L+uXWAjRMyZFMH8uU11oGCKECEKehWiFcZd8cvZt+fvYpJufv2BTMnswzXKVi0avxWIyilzODRvtwJ0Yy/Z8bBnCNJH90JNlZywWEDLPY+BaDwVjmu8g5Hyoo6K17LyV8Urc8WHE9uoHU7cfvJTaZI4rq1mi2OxAoddikUfYF9ruSldtMk1OPrRTVvKkGtyxQz5n9fXoOnMrFvPA8A6SFL78VdeLVmXXrMv6CJ2/Mv49ZnkG49sPXppmfkxzKMAphAhCnoWYGXlDwzhmUY2QNouTvUG63hOkK3RaLLLGyFlGl115e9Qlb2dFsuWt66l2TNf8JYxpy5bjz11dvhcTDUOEEEF02rMoS9m8ufLsQoQzV2IhJusTxodUVcbn8fCl6KSZZpiTVVMgusPCi8Xu+LuBcfiNU7RhG3o4S1w5tBfozoUUT6Qq8sKy1nOZrlVeGcImbI03Ql4mvfVfHbGrK6vX4dDVweMpl+kc/3/i8yImPj+arlNZeLFomuRdVlH3flEkepP1EikcusrqRy4DsAVwqb/r+kosRGsEXWwx0d15l0sHOLjVj/j8z668vft4K+Xvux7H2P945dDkVISYpLe5WyLfcBZPYjFG2Q7OO+srmR7DeGfrkCDqtD0Tshi3rUrac7nhLtXzykvv3gzA1js3pHoWseiND1OyhxftiKTEogRZi+VkiUA8iUpzHsQ80A8/TgjROr30LEJc6S78Y/EU+HgsGRLgTHop8dTleWIkTjEeh8ghDhLGbJGWXRidqp3VVm9y8aalqVrfhcRMrnBgGKdIxFlC4jXpi0y11ymrC9eU6Amh2ZwRwZsiBjMuDDAaJIwZDxYOxWP04poMBo4KQ7KvSX6QMDvgmM2eqVOmXUNisQCMpG/reh2lgRcOiYUQJSkXsN5D0hMZejn527pIt60TQnQGeRYdRilX0SUKxcLMjgC/C/wMcB047e7/2cxWgK8DtwKvAV929x9Hx5wCjgPXgK+4+zcbsV6IHpBWbZneIDrd0S8qyupSBecO8G/c/QUzuwF43syeAv4l8LS7P2JmJ4GTwENmdjtwP3AHcAvwLTP7uLtfa+ZfEE0xXoS2vbY01bqcI3MpEunStKwHpGc+Pn3DcKHOFy5/tJQNoh4KxcLdN4HN6PFlM3sJWAfuBe6KdjsDPAM8FG1/3N2vAq+a2QXgGPDtuo3vK3lTuOsu9057b2UyqhPaySyrJmJ0e9Gs3G6EFkvFLMzsVuBTwHeBmyMhwd03zeymaLd14DuJwzaibY2RNp8jtKdC6ErZusBE02S1IiiKXc2qDWSwZJnZh4A/BH7V3X+St2vKtolPwcxOmNlzZvbc+4RX8gnRFNurexr37PpMkGdhZvsZCMXvu/sfRZvfNrO1yKtYAy5G2zeAI4nDDwNvjb+mu58GTgPcaCuF3V3GvQc1jQ1n2pmn4yfHMs0NYZJxijg+8ZnlC7vbjh3Yv/v43PIF/hufGzl+650bGrFLDAnJhhjwO8BL7v5o4k9PAg8Aj0S/n0hsf8zMHmUQ4DwKnKvT6Hkg6w5Wd1fvtGnleV96UkhSA5wlp9nX1R3r3NX3Afju9sd44fJHefaVj43u0OOmMn0hxLP4LPAvgD83sx9E2/4dA5E4a2bHgdeB+wDc/byZnQVeZJBJeVCZEFGGvGzHC5c/OugPMSYOyUldTbUhzFo6Ip/0NGeZtXW7sopcSDbk/5IehwC4O+OYh4GHK9hVSHJYktVnooiufAnzRtZFtX1pD3Bw9/nlQ8O/Pbt6aOJ1nj30sYlt8czNW14avf8sb/5093HaUGn0HMk+NwaeXcikslGy960/BtJWsF0VnCKT8ZNSyyUsNhKLjjFvJd75rnvegs/p3bD70ndzHtEnL4QIQp5FzeTFSuJxc9ybczakrIfRYAp0UQhdvCl76Ba+tm5XviuJxZTEopAMiA2fjzO5/kTzC+50Z3EiMR9ILMbIbOv/xsbIh5VcOUxMx2h8Jk1k89fLqNWWlNZ6Ze7oWYVuWRdYH4PFvRGL8S+j7PoesyItQDnpZi4xi3DR8tb1VBe3K26t6Be9EYt5Jjkk6fLQoarI5N9N0zy0SUHN6t1QxbbRIeXQljJ3/9ALKb7p9fHCUzZECBFEHwUuiJDqO7njQoTTW7FIxjB6+08I0SN0nYmZkefJpccHUmpEMjIVhdPuE3+/cbzdwfpoZmuabMgiILGg/VZzXQ5qivYI7vY2o3N3LsSiaiOcfTW8RhsUF3aNFozFJO/ieSda6Mkak1a9GjpzMz5WM4G7y1yIhRBlmKjZiX7HItq14UdRu4XgFO/rxbvkIbEQC8dEi8bEnB1otroyvuDKem1d8Lh6JRZNVW3urK/U9kF04UsVogl6JRZCLCpd6DousRCi43TFW5VYVGSaOQqDMXEdd4r8av1567ol2kVi0RL13C0KouQZ75EUtD6mjEU7SCx6wPRrjKQfd+ObFQ0SC4nEQnSWvF4c8ohmj8RCLDyxAM26KGuaDEebwU6JBd2r2Atl2pWsdFcOp2zxVCjJps19CUR3Uiyyiq+SX1xfL3CRTuhMzzqFLm0tV1D38yw6KRain0zjIo+74kWl1n2d9DcPSCyEaIG0buIhtOnxSCyE6BB1zDBtSlDUsFcIEUTvPYuy0eq6VDdvMaLMY2a81knWXWq8rRwoDtAGaUHdLi8+JM9CCBFErzyLtDkNZf8B3UGFmI5eicW8UHYo1GXXVCwOGoYIIYIoFAszO2hm58zsh2Z23sx+I9q+YmZPmdkr0e8PJ445ZWYXzOxlM/t8k/+AEGI2hHgWV4HPufs/AD4J3GNmPwecBJ5296PA09FzzOx24H7gDuAe4LfMbG8TxgshZkdhzMLdHfhp9HR/9OPAvcBd0fYzwDPAQ9H2x939KvCqmV0AjgHfrtPwWVEl3TkeTA1J83ZpTsIsbBnvrVEUgN5p0JautK/rKkEBzsgzeB74GPBf3f27Znazu28CuPummd0U7b4OfCdx+Ea0bfw1TwAnAA6yPP1/IBaWaYV8vPV/GnFQuY5GufMiQkEBTne/5u6fBA4Dx8zsEzm7py2TNbE+n7ufdvc73f3O/RwIs1YI0RqlUqfu/jdm9gyDWMTbZrYWeRVrwMVotw3gSOKww8BbdRgrxKyp6hWEeCbJ4WnW/tPaUWddUaFYmNkq8H4kFEvAzwO/CTwJPAA8Ev1+IjrkSeAxM3sUuAU4CpyrzWIxM1TAJpKEeBZrwJkobrEHOOvu3zCzbwNnzew4g1UU7wNw9/NmdhZ4kUE86kF3v9aM+UKIWRGSDfkz4FMp298B7s445mHg4crWCSE6g8q9hZhD4jjIyAWuVdRFU6SlJhXHWFwkFoFM0+V5/MP9SV5ePyfaXaU4Ks3utG1pJ0LoflUIEZ+QBs5l2F5bYnnzvbmpf5gVnRQL3b1EF8kT7azO82U7hXd5hnEnxSImpEJvZ73clxFa9ZcmWGWr+br8xQtRFk1RF0IE0WnPQnSPulfoCjkBm1oVrItME5/azXw0PPFPnoUQIgh5FgtKyEpgdcy4HEcZiP4isRBzT1nRm8ad71IfkqbQMEQIEYTEQggRhIYhYqFJG6Ik4zeLMLwIpXNikSyaCkmZba8t1Vr8pJNDzBMj19CiTyRb3nyv1AW+SDl7Ieqk92IhhMimTk9ZYiEySauJaKL2QkxSxgOOv5PkcDxt28IPQ/rC9toS26vD5NPy1vUWrWmHwf9fXAw2cVxFgUp+7kOWJv62iN9JGSQWMyR5Mo7ftWN3scz0/LwZtEWrzGtGrCjLwohFyMIyMLyI4g8medz26h7dfXpI2nc2FOullG0ijYURiy6QPBnHA09q+CO6jsRC7BISG0gf/+cjb2w+kFgEMvAKlkq5qnnt1ebNk7hyKG3Vypg9iX1GxWZ5c/r3rDJsiL+LomrNqt/TPF1gmhsihAhinoRPiIVkODRcyt1WlV6JReaisWNpQM3v6CYHL/lE/CLku8pK84Z22w5l/BidR6P0SizapvQ8lMR4d14+6OxgZfaINq++JGjdkCDL+kWRmMVrm3SJefweRAvkZTxyU8Yp+48Xm03jJUy0DezYhddHJBZiockdysxZxqoqyoYIIYKQWAghgtAwpACVZYuuE8eLRuMySynbqtErscj6x6ukuCbLl1P6ACTfa8r3mUjLIeER/aJXYlGFtPLeAelR87pz7MrZ95MqC2nPG8ExCzPba2bfN7NvRM9XzOwpM3sl+v3hxL6nzOyCmb1sZp9vwnAhRDbLm+/t1gXFP1Up41n8CvAScGP0/CTwtLs/YmYno+cPmdntwP3AHcAtwLfM7OPufq2ytQtG0V1tZ30lu6o1Gvb8JKO4J+24tBmlfZwxOt6VbED5Dl37CK/xSLuQqngbZQoAd3uwNDARLkmQZ2Fmh4FfAH47sfle4Ez0+AzwpcT2x939qru/ClwAjtVjrhCiLUKHIV8Ffg1I3mZudvdNgOj3TdH2deCNxH4b0TYhRI8pFAsz+yJw0d2fD3zNtMYGnvK6J8zsOTN77n2uBr60EKItQmIWnwV+0cy+ABwEbjSz3wPeNrM1d980szXgYrT/BnAkcfxh4K3xF3X308BpgBttZUJMhJgHsmIGIRdeaLxht19swxmZQpvd/RRwCsDM7gL+rbv/spn9J+AB4JHo9xPRIU8Cj5nZowwCnEeBc/WbLkS3ye2+HhA4Lbo4Z52urVJn8Qhw1syOM1i+5D4Adz9vZmeBFxlMKnywq5mQvKh5XrQ8NEqdd0KMv0baiRUyjbmJDMZoi7w6ZgQkOmiP9x5J2XtnfUV1KR2klFi4+zPAM9Hjd4C7M/Z7GHi4om2NMN4XM6soa17ISsGli2H8WdS78E7eFPUQyh6zCN9rG2gimRAiiLkr9y67RmRR8U4aoat5JQufdGcTTTGr2MXciYUII39oUq/IKf4wIHkjK3MjSR7X5mfZabHQSTYkbwGg8biCvJhuMhkvq+E1Aye6AfO3ivosF3VZZoXt1Q+ObAtbcSusvfroa03OTZhVflzMlpChcMjqb7GgjHghyUWRprCtCp0TC9FPcj2fglW/RD9QNkQIEYQ8izlne22pdPYmb/gU0/W4yLh98miqI7EQMyNUuMoIXNbx0H1B6xsSC7HQTLOA0aIisRCiBfro9UgsFpS8tnphE8nqW517nsnLEoWl6Ue/jzjVv7x1vfxQbd7qLIQouwD1xPE5PSlnTb4HESq4kxP74s9oljU6EgvRW7LiDbsBzsBgah0sM//T6lVnIYQIQp6FEDn0MRDZFPIshBBByLMQI1w5ZFxZTc5izbqflLnPJCtD1cGqLF1Z6GmhxWLfm++yXFOAax4YnJR7SArBwUvpjdfTTuCsVOBoSjas2VCXMhoxRcHQtFmiua8XmDrdff2WRWOhxWKeGZ3nofJnUZ2FF4uqF1CV1Fuc+staU7PotbdX92TercvetYQoYu7EomyRStUPID+XryFOG8Q3gJAhTFEtRpnuVvG+ZYYjo9WyRYzOBk42UApCFZyiLxy85MEtAOMLL0/8dfLOFn3ePSf7ApRXI+pFYtER8iZ2ZTF0YUNjFsP3yHrtotRpVnZElGeYfSqzf3tILHKYpvHqtAHT9OOKvIPBiVb3SXRwKzt1mpzIlGTfm++yPEVviC6lRtumOH4xeVMInXlax+fcCbGwD+xn38+UaGleI3mCUJcQtEH6iben4O+MeRWjx4w+X+rV5yGq0wmxiNGCuKIMdUzPLrNSeaculhZY6P9/Z32lZHOSkLhCcVxgfL+6FpwJJRmXSA45qjKN0Le9Zkrb7z86hMz/LtIC2aGf+SyEde4p05wkJK4wuk/OWhpy4UXPWHixEKIO0tYxLVtFW0eAs8lhvMSiJUKq+5InW17GY/ykvHLIuPKzk57LlUsHBg8OXWH1I5cB2Dp0w8g+8XaArXdu4AoHxl5lT+J3+aGUYlL9RWLRIHlDjaC1LiOBmGqex6Xxizxs/60qryHmGomFEHPIeEnAPqg8N0RTE4UQQQR5Fmb2GnAZuAbsuPudZrYCfB24FXgN+LK7/zja/xRwPNr/K+7+zdotp/xqUuOuf0g59WgFY1jqtE+Zjt3YRUvvv+9IfjFe26lNMaSMZ/GP3f2T7n5n9Pwk8LS7HwWejp5jZrcD9wN3APcAv2Vme0PeQMEvIbpLlWHIvcCZ6PEZ4EuJ7Y+7+1V3fxW4AByr8D5CiA4QGuB04H+ZmQP/3d1PAze7+yaAu2+a2U3RvuvAdxLHbkTbRjCzE8CJ6OlP//T1r74DXCplfcWATQUOUdbWdpmw97V27Aihu5/t5Pk2tLW9czGddHv+XpWXDBWLz7r7W5EgPGVmf5Gzb1plycS85khwTu8eZPZcYojTafpkK/TLXtnaHGb2XJXjg4Yh7v5W9Psi8McMhhVvm9laZMQacDHafQM4kjj8MPBWFSOFEO1TKBZm9kEzuyF+DPxT4EfAk8AD0W4PAE9Ej58E7jezA2Z2G3AUOFe34UKI2RIyDLkZ+GMzi/d/zN3/1My+B5w1s+MMRkj3Abj7eTM7C7wI7AAPuvu1gPc5XbxLZ+iTrdAve2Vrc1Sy19zVJk0IUYwqOIUQQbQuFmZ2j5m9bGYXzOxk2/YAmNnXzOyimf0osW3FzJ4ys1ei3x9O/O1UZP/LZvb5Gdt6xMz+j5m9ZGbnzexXumqvmR00s3Nm9sPI1t/oqq2J999rZt83s2/0wNbXzOzPzewHceajVnvdvbUfYC/wV8DfBT4A/BC4vU2bIrv+EfBp4EeJbf8ROBk9Pgn8ZvT49sjuA8Bt0f+zd4a2rgGfjh7fAPxlZFPn7GWQVv9Q9Hg/8F3g57poa8Lmfw08Bnyjy+dBZMNrwKGxbbXZ27ZncQy44O5/7e5/CzzOoAK0Vdz9WWC89ryTFavuvunuL0SPLwMvMSiC65y9PuCn0dP90Y930VYAMzsM/ALw24nNnbQ1h9rsbVss1oE3Es9Tqz07wkjFKpCsWO3E/2BmtwKfYnDH7qS9kVv/AwZ1OU+5e2dtBb4K/BqQ7DzUVVthWGn9fFQhDTXa23Y/i6Bqz47Tif/BzD4E/CHwq+7+kyjVnbpryraZ2euDNPonzezvMEjJfyJn99ZsNbMvAhfd/XkzuyvkkJRtsz4Paq+0TtK2Z9Gnas/OVqya2X4GQvH77v5H0ebO2gvg7n8DPMNgZnIXbf0s8ItRe4bHgc+Z2e911Fag+UrrtsXie8BRM7vNzD7AYGr7ky3blEUnK1Zt4EL8DvCSuz/aZXvNbDXyKDCzJeDngb/ooq3ufsrdD7v7rQzOy//t7r/cRVthRpXWs4zWZkRwv8Aggv9XwK+3bU9k0x8Am8D7DBT4OPARBn07Xol+ryT2//XI/peBfzZjW/8hA/fxz4AfRD9f6KK9wN8Hvh/Z+iPg30fbO2frmN13McyGdNJWBhnFH0Y/5+NrqU57VcEphAii7WGIEKInSCyEEEFILIQQQUgshBBBSCyEEEFILIQQQUgshBBBSCyEEEH8f3NgwaK7jJ3yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initialize 1000 squares with center (x,y) and size (x,y) following an uniform distribution\n",
    "np.random.seed(1066)\n",
    "arr = np.zeros(shape)\n",
    "center = np.random.random(400).reshape((200, 2)) * arr.shape[0]\n",
    "size = np.random.random(400).reshape((200, 2)) * arr.shape[0] // 5\n",
    "\n",
    "for i in range(200):\n",
    "    arr += generate_square(shape, center[i], size[i])\n",
    "\n",
    "plt.imshow(arr)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "293181b8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAD4NJREFUeJzt3VuMXVd9x/HvL8ZxKGASN8R1bEtxI7fSBLUBuQYpVZUS2rgBEV5SGYnWlSL5JZVArQR2kdryECn0IeIpUi1AdQXBGnFRrAi1OAaEKlEcBxKIbUwmJE1GMbGAhstDTS7/Ppw95Niey/LMuQ18P9Lo7LO89jn/I41/s/baa++TqkKSlnLZuAuQtDoYFpKaGBaSmhgWkpoYFpKaGBaSmgwtLJLsSnI6yUySfcN6H0mjkWGss0iyBvg+8GfALPAw8L6qOjnwN5M0EsMaWewEZqrqB1X1S+AQcPuQ3kvSCLxmSK+7GXi27/ks8LaFOl+edXUFrxtSKZIAfs7//qiq3rTc/YcVFpmn7bzjnSR7gb0AV/BbvC23DKkUSQAP1ef+ZyX7D+swZBbY2vd8C/Bcf4eqOlBVO6pqx1rWDakMSYMyrLB4GNieZFuSy4HdwOEhvZekERjKYUhVvZTkb4H/BNYAn6qqE8N4L0mjMaw5C6rqS8CXhvX6kkbLFZySmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkpoYFpKaGBaSmhgWkposGRZJPpXkbJLH+9o2JDmS5Inu8aq+f9ufZCbJ6SS3DqtwSaPVMrL4N2DXBW37gKNVtR042j0nyRSwG7ih2+e+JGsGVq2ksVkyLKrq68BPLmi+HTjYbR8E3tvXfqiqzlXVU8AMsHNAtUoao+XOWWysqjMA3eM1Xftm4Nm+frNd20WS7E1yPMnxFzm3zDIkjcqgJzgzT1vN17GqDlTVjqrasZZ1Ay5D0qAtNyyeT7IJoHs827XPAlv7+m0Bnlt+eZImxXLD4jCwp9veAzzQ1747ybok24DtwLGVlShpErxmqQ5JPgvcDFydZBb4J+AeYDrJncAzwB0AVXUiyTRwEngJuKuqXh5S7ZJGaMmwqKr3LfBPtyzQ/27g7pUUJWnyLBkW0kpdduPUSN7nlUdPjuR9flMZFhqqy26c4oWp9SN5rytZOpQMlOXz2hBJTRxZTJhBD9nH8Ze0/zO8MLWen14/qr9JS49g+kcfjjIujSMLSU0MC0lNDAtJTQwLSU0MC0lNDAtJTQwLSU0MC0lNDAtJTQwLSU0MC0lNvDZkgg3ias0Lr8T0eggtlyMLSU0cWWhk/u/6YX/lQ8td4l8drc2NuhxttTEsJtTgLu0+/1Cm5QYxS/E/128mD0MkNTEsJDUxLCQ1cc5CzeZO5S4179F/yven11/2q4nNqeuG++V0J7kWgCue9Oswh8GRhaQmjiw0Nn/5Ow8ve9/pH/7RACtRC8NihQZ9N+5Bf8fGlSd/dt7zlZz2HMRpV61eHoZIamJYSGpiWEhqYlhIauIE52+A/knTlUxSnvc6F0yc6tefIwtJTQwLSU0MC0lNlgyLJFuTfDXJqSQnknyga9+Q5EiSJ7rHq/r22Z9kJsnpJLcO8wNIGo2WCc6XgL+vqm8leQPwSJIjwN8AR6vqniT7gH3Ah5NMAbuBG4BrgYeS/F5VvTycj6BhuXA16U+vv4w3PvnKmKrRuC0ZFlV1BjjTbf88ySlgM3A7cHPX7SDwNeDDXfuhqjoHPJVkBtgJfGPQxa9Wiy3pHszdsRZ/b89kaDku6TczyXXAW4BvAhu7IJkLlGu6bpuBZ/t2m+3aJK1izWGR5PXA54EPVtVif5oyT1vN83p7kxxPcvxFhn0jV0kr1bQoK8laekHxmar6Qtf8fJJNVXUmySbgbNc+C2zt230LcNFdT6rqAHAAYH02XBQmGpzlXnl68QKu4R3C/PX6H11S/+kfDqUMLWLJsEgS4JPAqaq6t++fDgN7gHu6xwf62u9Pci+9Cc7twLFBFv3rYKG5iUHfLn+++ZHFVnH2B8l8E5wX3i18Pv2fbdCf599/djUAJ5++dqCvq6W1jCxuAv4K+G6SR7u2f6AXEtNJ7gSeAe4AqKoTSaaBk/TOpNzlmRBp9Ws5G/JfzD8PAXDLAvvcDdy9grq0ivWfXn3jk2t/tf2D67ed1++fr391dDDdeH/OuRHFxi+vXbDPlSdfaHotXRovJNOCLp6fGOxdvLS6GBYTZvB3pnYRlQbDa0MkNTEsJDXxMGSZ5u7qvZy7cfd/8c4w+WU7GiTDQmPTH2Zz3ybWvo9zMaPmYYikJo4s1GzlS70XO2RrO2SaW8PhlbOj58hCUhPDQlITw0JSE+csNDKLzzO0nYKee42VfMGzlseRhaQmjixWsaklrtTsrV2Y7yxD2zeLLWfB2XLN3QPDGwJPLkcWkpoYFpKaGBaSmhgWkpoYFpKaeDZkhZZ3jcJ6Wq+FWMxSV2p6iboGybAYk8GcIlw8DBZ6j/6Ac3GTWhkWq8Dyv2Nk/jC50nzQMjhnIamJYSGpiYchq9hSE5gLzVk4T6HlMCw0EVrvgGXQjY+HIZKaGBaSmhgWkpoYFpKaOMG5TMuZaJv7FrNRWWgx15Xz1OHEoZbiyEJSE8NCUhMPQ8bg0q9UHd29MKWFOLKQ1GTJsEhyRZJjSR5LciLJR7v2DUmOJHmie7yqb5/9SWaSnE5y6zA/gKTRaBlZnAPeUVV/CNwI7ErydmAfcLSqtgNHu+ckmQJ2AzcAu4D7kqwZRvGSRmfJOYuqKuAX3dO13U8BtwM3d+0Hga8BH+7aD1XVOeCpJDPATuAbgyxcwzeObyr3FO7kapqzSLImyaPAWeBIVX0T2FhVZwC6x2u67puBZ/t2n+3aLnzNvUmOJzn+IkvdxEXSuDWFRVW9XFU3AluAnUnevEj3zPcS87zmgaraUVU71g7gfpSShuuSzoZU1Qv0Djd2Ac8n2QTQPZ7tus0CW/t22wIs/j17kibeknMWSd4EvFhVLyR5LfBO4GPAYWAPcE/3+EC3y2Hg/iT3AtcC24FjQ6hdQ+b8gfq1LMraBBzszmhcBkxX1YNJvgFMJ7kTeAa4A6CqTiSZBk4CLwF3VdXLwylf0qi0nA35DvCWedp/DNyywD53A3evuDpJE8MVnJKaGBaSmngh2SowjsVR0oUcWUhqYlhIamJYSGpiWEhqYlhIamJYSGpiWEhqYlhIamJYSGpiWEhqYlhIamJYSGpiWEhqYlhIamJYSGri/SxGyBvgajVzZCGpiWEhqYlhIamJYSGpiWEhqYlhIamJYSGpiWEhqYlhIamJYSGpiWEhqYlhIamJYSGpiWEhqYlhIalJc1gkWZPk20ke7J5vSHIkyRPd41V9ffcnmUlyOsmtwyhc0mhdysjiA8Cpvuf7gKNVtR042j0nyRSwG7gB2AXcl2TNYMqVNC5NYZFkC/Au4BN9zbcDB7vtg8B7+9oPVdW5qnoKmAF2DqZcSePSOrL4OPAh4JW+to1VdQage7yma98MPNvXb7Zrk7SKLRkWSd4NnK2qRxpfM/O01TyvuzfJ8STHX+Rc40tLGpeWG/beBLwnyW3AFcD6JJ8Gnk+yqarOJNkEnO36zwJb+/bfAjx34YtW1QHgAMD6bLgoTCRNliVHFlW1v6q2VNV19CYuv1JV7wcOA3u6bnuAB7rtw8DuJOuSbAO2A8cGXrmkkVrJVwHcA0wnuRN4BrgDoKpOJJkGTgIvAXdV1csrrlTSWKVq/EcA67Oh3pZbxl2G9GvtofrcI1W1Y7n7u4JTUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlKTprBI8nSS7yZ5NMnxrm1DkiNJnuger+rrvz/JTJLTSW4dVvGSRudSRhZ/WlU3VtWO7vk+4GhVbQeOds9JMgXsBm4AdgH3JVkzwJoljcFKDkNuBw522weB9/a1H6qqc1X1FDAD7FzB+0iaAK9p7FfAl5MU8K9VdQDYWFVnAKrqTJJrur6bgf/u23e2aztPkr3A3u7pLx6qz/0Y+NEyPsM4XM3qqRVWV73WOjy/v5KdW8Pipqp6rguEI0m+t0jfzNNWFzX0AufAr3ZKjvcd4ky01VQrrK56rXV45uYbl6vpMKSqnusezwJfpHdY8XySTV0Rm4CzXfdZYGvf7luA51ZSpKTxWzIskrwuyRvmtoE/Bx4HDgN7um57gAe67cPA7iTrkmwDtgPHBl24pNFqOQzZCHwxyVz/+6vqP5I8DEwnuRN4BrgDoKpOJJkGTgIvAXdV1csN73Ng6S4TYzXVCqurXmsdnhXVm6qLphMk6SKu4JTUZOxhkWRXt9JzJsm+cdcDkORTSc4mebyvbSJXrCbZmuSrSU4lOZHkA5Nab5IrkhxL8lhX60cntda+91+T5NtJHlwFtQ53pXVVje0HWAM8CfwucDnwGDA1zpq6uv4EeCvweF/bvwD7uu19wMe67amu7nXAtu7zrBlhrZuAt3bbbwC+39U0cfXSO63++m57LfBN4O2TWGtfzX8H3A88OMm/B10NTwNXX9A2sHrHPbLYCcxU1Q+q6pfAIXorQMeqqr4O/OSC5olcsVpVZ6rqW932z4FT9BbBTVy91fOL7una7qcmsVaAJFuAdwGf6GueyFoXMbB6xx0Wm4Fn+57Pu9pzQpy3YhXoX7E6EZ8hyXXAW+j9xZ7Ierth/aP01uUcqaqJrRX4OPAh4JW+tkmtFV5daf1It0IaBlhv6wrOYWla7TnhJuIzJHk98Hngg1X1s+5U97xd52kbWb3VO41+Y5Ir6Z2Sf/Mi3cdWa5J3A2er6pEkN7fsMk/bqH8PBr7Sut+4RxarabXnxK5YTbKWXlB8pqq+0DVPbL0AVfUC8DV6VyZPYq03Ae9J8jS9w+N3JPn0hNYKDH+l9bjD4mFge5JtSS6nd2n74THXtJCJXLGa3hDik8Cpqrp3kutN8qZuREGS1wLvBL43ibVW1f6q2lJV19H7vfxKVb1/EmuFEa20HuVs7QIzuLfRm8F/EvjIuOvpavoscAZ4kV4C3wn8Nr37djzRPW7o6/+Rrv7TwF+MuNY/pjd8/A7waPdz2yTWC/wB8O2u1seBf+zaJ67WC+q+mVfPhkxkrfTOKD7W/ZyY+780yHpdwSmpybgPQyStEoaFpCaGhaQmhoWkJoaFpCaGhaQmhoWkJoaFpCb/D537Tg/xUE0rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# One by one, generate the uniform distribution of square and compute their overlap using 3 metrics\n",
    "# Dice, Weigthed-Dice and Correlation, those with good (enough) score should be kept\n",
    "\n",
    "np.random.seed(1066)\n",
    "arr = np.zeros(shape)\n",
    "center = np.random.random(400).reshape((200, 2)) * arr.shape[0]\n",
    "size = np.random.random(400).reshape((200, 2)) * arr.shape[0] // 5\n",
    "\n",
    "def compute_dice_voxel(density_1, density_2):\n",
    "    overlap_idx = np.nonzero(density_1 * density_2)\n",
    "    numerator = 2 * len(overlap_idx[0])\n",
    "    denominator = np.count_nonzero(density_1) + np.count_nonzero(density_2)\n",
    "\n",
    "    if denominator > 0:\n",
    "        dice = numerator / float(denominator)\n",
    "    else:\n",
    "        dice = np.nan\n",
    "\n",
    "    overlap_1 = density_1[overlap_idx]\n",
    "    overlap_2 = density_2[overlap_idx]\n",
    "    w_dice = np.sum(overlap_1) + np.sum(overlap_2)\n",
    "    denominator = np.sum(density_1) + np.sum(density_2)\n",
    "    if denominator > 0:\n",
    "        w_dice /= denominator\n",
    "    else:\n",
    "        w_dice = np.nan\n",
    "\n",
    "    return dice, w_dice\n",
    "\n",
    "\n",
    "def compute_correlation(density_1, density_2):\n",
    "    indices = np.where(density_1 + density_2 > 0)\n",
    "    if np.array_equal(density_1, density_2):\n",
    "        density_correlation = 1\n",
    "    elif (np.sum(density_1) > 0 and np.sum(density_2) > 0) \\\n",
    "            and np.count_nonzero(density_1 * density_2):\n",
    "        density_correlation = np.corrcoef(density_1[indices],\n",
    "                                          density_2[indices])[0, 1]\n",
    "    else:\n",
    "        density_correlation = 0\n",
    "\n",
    "    return max(0, density_correlation)\n",
    "\n",
    "arr = np.zeros(shape)\n",
    "for i in range(200):\n",
    "    curr = generate_square(shape, center[i], size[i])\n",
    "    dice, w_dice = compute_dice_voxel(arr_smooth, curr)\n",
    "    corr = compute_correlation(arr_smooth, curr)\n",
    "    # This value of 0.1 control how strict the search is\n",
    "    if dice + w_dice + corr > 0.1:\n",
    "        arr += curr\n",
    "\n",
    "plt.imshow(arr)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "600b08d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAEBxJREFUeJzt3V2MXVd5xvH/E8exKWASNyQY21LSyK3koDYgY5BSVSmhjQuIcJPKSLSuFMk3qQRqJXCK1JaLSKEXEVeRagGqEYTI4kOxItTiGBCqREkcSCB2CHFImlg2sYCGjwtMPt5enD3J8cyZmTUz52vi/08anX2W9z7nPdL4mbXWXnufVBWStJgLJl2ApNXBsJDUxLCQ1MSwkNTEsJDUxLCQ1GRkYZFkV5LHkpxIsm9U7yNpPDKKdRZJ1gA/Bv4COAk8AHywqo4P/c0kjcWoehY7gRNV9ZOq+h1wN3DjiN5L0hhcOKLX3Qw80/f8JPCO+Xa+KOtqPa8dUSmSAH7N//2sqt643ONHFRYZ0HbOeCfJXmAvwHp+j3fk+hGVIgngvvrS/67k+FENQ04CW/uebwFO9e9QVfurakdV7VjLuhGVIWlYRhUWDwDbklyZ5CJgN3BoRO8laQxGMgypqheS/D3wX8Aa4LNVdWwU7yVpPEY1Z0FVfQ342qheX9J4uYJTUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlKTRcMiyWeTnEnySF/bxiSHkzzePV7S92+3JjmR5LEkN4yqcEnj1dKz+A9g16y2fcCRqtoGHOmek2Q7sBu4ujvmziRrhlatpIlZNCyq6tvAL2Y13wgc6LYPAB/oa7+7qs5W1ZPACWDnkGqVNEHLnbO4vKpOA3SPl3Xtm4Fn+vY72bXNkWRvkqNJjj7P2WWWIWlchj3BmQFtNWjHqtpfVTuqasda1g25DEnDduEyj3s2yaaqOp1kE3Cmaz8JbO3bbwtwaiUFamkuuGb7yF77pYeOj+y1Nf2W27M4BOzptvcA9/S1706yLsmVwDbg/pWVKGkaLNqzSPJF4Drg0iQngX8BbgcOJrkZeBq4CaCqjiU5CBwHXgBuqaoXR1S7pDFaNCyq6oPz/NP18+x/G3DbSoqSNH1cwSmpyXInODUGLZOVTjpqXOxZSGpiz2LK9Pcmntu+YdH9L+bc3sfsYy4+/qvhFKbznj0LSU0MC0lNDAtJTQwLSU0MC0lNDAtJTQwLSU0MC0lNDAtJTVzBeZ5qWR062+zVoqPgtS7Ty7A4Ty1nGbj/kc9vhsUU++VVLaPEc3sIc4/xWhENh3MWkprYs9DYjPJmwivh8KqNYTGlntu+gd9e1fJ9Kud+jcLcY2Z/zcLSJzZnzExw+p/r/OQwRFITw0JSE4chrzLrnzh32PGGJ1465/lKzobMrM1YynoLhyyvHoaFmr1yWrZ93qM/WJazEGw5PD08Gg5DJDWxZ6GJaVt01jP7LM/6J9bNGWJptAyLFRr22oFhd9Vnd8lXMofwhu3vXGk5WsUchkhqYlhIamJYSGrinIWaeUry/GZYnAf6J01XcgObc17H4DjvOAyR1MSwkNTEYcgUm32dhzRJi/YskmxN8s0kjyY5luTDXfvGJIeTPN49XtJ3zK1JTiR5LMkNo/wAksajpWfxAvCPVfW9JK8HHkxyGPg74EhV3Z5kH7AP+FiS7cBu4GrgzcB9Sf6wql4czUfQqMxeTfrLqy5Y8hLr/teYvby77eY+PduvOPXy9vGn3rykGjQci4ZFVZ0GTnfbv07yKLAZuBG4rtvtAPAt4GNd+91VdRZ4MskJYCfwnWEXv1ottKR7KddLLPe9PZOh5VjSb2aSK4C3At8FLu+CZCZQLut22ww803fYya5N0irWHBZJXgd8GfhIVS30pykD2mrA6+1NcjTJ0edp745KmoymsyFJ1tILii9U1Ve65meTbKqq00k2AWe69pPA1r7DtwCnmKWq9gP7ATZk45ww0fAs98rTuQu4RjeE6Z+TmPHXb3rg5e2/3fCzl7c/96ZL+VdupP9mxF6uPnqLhkWSAJ8BHq2qO/r+6RCwB7i9e7ynr/2uJHfQm+DcBtw/zKJfDeabm1jKpF+LQfMjC63i7A+SQROcS707eP/nnP3Z+gOiPxhm9AfEjM/96lIO/vTtnlaegJaexbXA3wA/TPJQ1/ZP9ELiYJKbgaeBmwCq6liSg8BxemdSbvFMiLT6tZwN+W8Gz0MAXD/PMbcBt62gLklTxhWcmtfc+Ynx3HBX08mwmDLDH4tPduLv3InHcz/bcV5ZXHWQtw84eu48xsGfDtpP4+CFZJKa2LMYspYb7v7yqguGftZjfoN6Kq7i1NIZFss0c1fvwacXF/bbq84OXFcwbF5DoWEyLDQx/fMz/fMXMwbNY8wE4PrRlaV5GBZD1raScN3A/xzDNt8X8TgE0XIYFlNgtVx+vfKQWWg+Z+7cyqBAnemNDPMLn9XGsyGSmhgWkpoYFpKaOGehsVl4XmHQfMbceYyZuYphfuGz2hgWU2CaJzWlGYbFKrbYwq7e2YTBKzhnLPTXvmU1ar9BC9Jab0ozc6w3sZlezllIamJYSGpiWEhqYlhIamJYSGri2ZAVWt41CRsYfJZiaRa7GM07YGuYDIsJGc4pwoXDYL736A84FzOplWGxCiz/O0YGh8nF5oOWwbDQ1FroXhz2iMbPCU5JTexZrGKLTWDON2fhX2Uth2GhqTDf1aSzGXSTY1hoaJZzhmfu5O3CF69dzHYDY0Kcs5DUxLCQ1MSwkNTEOYtlWsm4+WK2L/EIv71ck2fPQlITw0JSE4chE7DUIczShy2LvP9QX03nC3sWkposGhZJ1ie5P8nDSY4l+UTXvjHJ4SSPd4+X9B1za5ITSR5LcsMoP4Ck8WjpWZwF3lVVfwJcA+xK8k5gH3CkqrYBR7rnJNkO7AauBnYBdyZZM4riJY3PonMWVVXAb7qna7ufAm4EruvaDwDfAj7Wtd9dVWeBJ5OcAHYC3xlm4Rq9cXwz+ex7ayw2n+N8y+Q0zVkkWZPkIeAMcLiqvgtcXlWnAbrHy7rdNwPP9B1+smub/Zp7kxxNcvR5FruJi6RJawqLqnqxqq4BtgA7k7xlgd0z6CUGvOb+qtpRVTvWDuF+lJJGa0lnQ6rqOXrDjV3As0k2AXSPZ7rdTgJb+w7bAiz8PXuSpt6icxZJ3gg8X1XPJXkN8G7gk8AhYA9we/d4T3fIIeCuJHcAbwa2AfePoHaNmJeCq1/LoqxNwIHujMYFwMGqujfJd4CDSW4GngZuAqiqY0kOAseBF4BbqurF0ZQvaVxazob8AHjrgPafA9fPc8xtwG0rrk7S1HAFp6QmhoWkJl5Itgo40ahpYM9CUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUpPmsEiyJsn3k9zbPd+Y5HCSx7vHS/r2vTXJiSSPJblhFIVLGq+l9Cw+DDza93wfcKSqtgFHuuck2Q7sBq4GdgF3JlkznHIlTUpTWCTZArwX+HRf843AgW77APCBvva7q+psVT0JnAB2DqdcSZPS2rP4FPBR4KW+tsur6jRA93hZ174ZeKZvv5Ndm6RVbNGwSPI+4ExVPdj4mhnQVgNed2+So0mOPs/ZxpeWNCkXNuxzLfD+JO8B1gMbknweeDbJpqo6nWQTcKbb/ySwte/4LcCp2S9aVfuB/QAbsnFOmEiaLov2LKrq1qraUlVX0Ju4/EZVfQg4BOzpdtsD3NNtHwJ2J1mX5EpgG3D/0CuXNFYtPYv53A4cTHIz8DRwE0BVHUtyEDgOvADcUlUvrrhSSROVqsmPADZkY70j10+6DOlV7b760oNVtWO5x7uCU1ITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUhPDQlITw0JSk6awSPJUkh8meSjJ0a5tY5LDSR7vHi/p2//WJCeSPJbkhlEVL2l8ltKz+POquqaqdnTP9wFHqmobcKR7TpLtwG7gamAXcGeSNUOsWdIErGQYciNwoNs+AHygr/3uqjpbVU8CJ4CdK3gfSVPgwsb9Cvh6kgL+var2A5dX1WmAqjqd5LJu383A//Qde7JrO0eSvcDe7ulv7qsv/Rz42TI+wyRcyuqpFVZXvdY6On+0koNbw+LaqjrVBcLhJD9aYN8MaKs5Db3A2f/yQcnRviHOVFtNtcLqqtdaR2dmvnG5moYhVXWqezwDfJXesOLZJJu6IjYBZ7rdTwJb+w7fApxaSZGSJm/RsEjy2iSvn9kG/hJ4BDgE7Ol22wPc020fAnYnWZfkSmAbcP+wC5c0Xi3DkMuBryaZ2f+uqvrPJA8AB5PcDDwN3ARQVceSHASOAy8At1TViw3vs3/xXabGaqoVVle91jo6K6o3VXOmEyRpDldwSmoy8bBIsqtb6Xkiyb5J1wOQ5LNJziR5pK9tKlesJtma5JtJHk1yLMmHp7XeJOuT3J/k4a7WT0xrrX3vvybJ95PcuwpqHe1K66qa2A+wBngC+APgIuBhYPska+rq+jPgbcAjfW3/BuzrtvcBn+y2t3d1rwOu7D7PmjHWugl4W7f9euDHXU1TVy+90+qv67bXAt8F3jmNtfbV/A/AXcC90/x70NXwFHDprLah1TvpnsVO4ERV/aSqfgfcTW8F6ERV1beBX8xqnsoVq1V1uqq+123/GniU3iK4qau3en7TPV3b/dQ01gqQZAvwXuDTfc1TWesChlbvpMNiM/BM3/OBqz2nxDkrVoH+FatT8RmSXAG8ld5f7Kmst+vWP0RvXc7hqpraWoFPAR8FXuprm9Za4ZWV1g92K6RhiPW2ruAclabVnlNuKj5DktcBXwY+UlW/6k51D9x1QNvY6q3eafRrklxM75T8WxbYfWK1JnkfcKaqHkxyXcshA9rG/Xsw9JXW/Sbds1hNqz2ndsVqkrX0guILVfWVrnlq6wWoqueAb9G7Mnkaa70WeH+Sp+gNj9+V5PNTWisw+pXWkw6LB4BtSa5MchG9S9sPTbim+UzlitX0uhCfAR6tqjumud4kb+x6FCR5DfBu4EfTWGtV3VpVW6rqCnq/l9+oqg9NY60wppXW45ytnWcG9z30ZvCfAD4+6Xq6mr4InAaep5fANwO/T+++HY93jxv79v94V/9jwF+NudY/pdd9/AHwUPfznmmsF/hj4PtdrY8A/9y1T12ts+q+jlfOhkxlrfTOKD7c/Ryb+b80zHpdwSmpyaSHIZJWCcNCUhPDQlITw0JSE8NCUhPDQlITw0JSE8NCUpP/B/ZUkwzd3UqgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# One by one, generate the uniform distribution of square and compute their overlap using a distance to barrycenter\n",
    "\n",
    "np.random.seed(1066)\n",
    "arr = np.zeros(shape)\n",
    "center = np.random.random(400).reshape((200, 2)) * arr.shape[0]\n",
    "size = np.random.random(400).reshape((200, 2)) * arr.shape[0] // 5\n",
    "\n",
    "pos_smooth = np.where(arr_smooth)\n",
    "val_smooth = arr_smooth[pos_smooth]\n",
    "center_smooth = np.average(pos_smooth, axis=1, weights=val_smooth)\n",
    "for i in range(200):\n",
    "    curr = generate_square(shape, center[i], size[i])\n",
    "    pos = np.where(curr)\n",
    "    if len(pos[0]) > 1:\n",
    "        pos = np.average(pos, axis=1)\n",
    "        # This value represent a distance in pixel\n",
    "        if np.linalg.norm(pos - center_smooth) < 100:\n",
    "            arr += curr\n",
    "\n",
    "plt.imshow(arr)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7e83ffa",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}