Racunska geometrija.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4c3f380b-2024-446c-b112-45b01e53a4a6",
   "metadata": {},
   "source": [
    "# Računska geometrija\n",
    "\n",
    "Računska geometrije je področje algoritmov in podatkovnih struktur, ki se ukvarja z učinkovitim reševanjem geometrijskih problemov. Ti vključujejo delo s točkami, daljicami, večkotniki in drugimi geometrijskimi objekti ter relacijami med njimi, kot sta npr. razdalja ali vsebovanost. Računska geometrija je očiten del računalniške grafike, vida, robotike. Manj očitno pa je prisotna tudi v številnih drugih problemih, ki dopuščajo geometrijsko formulacijo.\n",
    "\n",
    "Omejili se bomo na reševanje problemov v ravnini, saj se problemi v višjih dimenzijah običajno dodatno zakomplicirajo. Poleg tega je reševanje ravninskih problemov enostavno za vizualizacijo. Kljub temu pa moramo biti pri reševanju pozorni na številne *posebne primere*, kot so kolinearne točke, sovpadanje točk, vzporedne daljice, ... Pomembna ovira je tudi *računska natančnost*. Če imamo opravka s celoštevilskimi objekti, želimo rešiti problem z uporabo celih števil, da ne vpeljemo računske napake, ki bi lahko povzročila povsem napačen rezultat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "60bd11ab-fda6-401a-afd0-b2b7f5ea1ec8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <iostream>\n",
    "#include <cmath>\n",
    "#include <vector>\n",
    "#include <algorithm>\n",
    "using namespace std;\n",
    "\n",
    "typedef pair<int,int> PII;\n",
    "typedef vector<PII> VII;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6a36633e-5315-417d-9ffe-f674a9ea88b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "template<class A, class B> \n",
    "ostream& operator<<(ostream& os, pair<A,B> &p) {\n",
    "    os << \"(\" << p.first << \", \" << p.second << \")\";\n",
    "    return os;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4240b3f3-93ff-45b6-96c4-e803438f85a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "template<typename T>\n",
    "void print(const vector<T> &sez) {\n",
    "    for (T x : sez) cout << x << \" \";\n",
    "    cout << endl;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deeba58f-502e-4498-a45c-453b0f6c05b6",
   "metadata": {},
   "source": [
    "## Razdalje in presečišča\n",
    "\n",
    "Razdalje in presečišča so najbolj osnovni koncepti, ki jih moramo obvladati. Ne vključujejo kakšnih novih algoritmičnih prijemov, vendar služijo kot ponovitev geometrije in linearne algebre. Za naštete probleme seveda obstaja več formul. Ogledali si bomo najbolj enostavne, ki jih lahko izpeljemo brez prepisovanja iz kakšnega učbenika. Glede na predstavitev premic imamo različne pristope. Predpostavili bomo, da imamo premico $P$ predstavljeno s točko $P_0$ in smernim vektorjem $V_P$. \n",
    "\n",
    "Razdalje:\n",
    "\n",
    "- Točki $S$ in $T$: Pitagorov izrek poznamo še iz osnovne šole.\n",
    "- Točka $S$ in premica $P$: Izračunamo projekcijo $S'$ točke $S$ na premico $P$ in izračunamo razdaljo med $S$ in projekcijo $S'$. **Projekcijo** točke na premico izračunamo s pomočjo **skalarnega produkta**: $\\operatorname {proj} _b a = \\frac{a \\cdot b}{||b||^{2}}b$. Če delamo projekcijo na enotski smerni vektor premice, je dolžina projekcije kar enaka skalarnemu produktu.\n",
    "- Točka $S$ in daljica $AB$: Izračunamo projekcijo točke na nosilko (premico) daljice. Če je projekcija izven daljice, bo najkrajša razdalja do krajišča $A$ ali $B$, sicer pa do projekcije $S'$.\n",
    "- Daljici $AB$ in $CD$: Predpostavimo, da se daljici ne sekata, sicer je odgovor 0. Najkrajša razdalja med daljicama bo enaka razdalji med enim izmed krajišč in drugo daljico. Izberemo najmanjšo izmed štirih možnosti.\n",
    "\n",
    "Presečišča:\n",
    "\n",
    "- Točka $S$ in premica $P$: Če je vektorski produkt vektorja $P_0 S$ in smernega vektorja premice $P$ enak 0, leži točka na premici.\n",
    "- Točka $S$ in daljica $AB$: Preverimo, ali točka leži na nosilki daljice in znotraj očrtanega pravokotnika (*bounding box*).\n",
    "- Premici $P$ in $R$: Če sta premici vzporedni, imamo neskončno ali nobenega presečišča. Sicer rešimo sistem enačb $P_0 + a V_P = R_0 + b V_R$ za obe koordinati.\n",
    "- Premica in daljica: Izračunamo presečišče premice in nosilke daljice ter preverimo, ali leži presečišče na daljici.\n",
    "- Daljici $AB$ in $CD$: Ugotoviti moramo, ali se daljici sploh sekata, nato pa uporabimo rešitev za izračun presečišča nosilk obeh daljic. Daljici se sekata natanko takrat, ko sta krajišči prve daljice $A$ in $B$ na nasprotnih straneh nosilke daljice $CD$ in obratno. **Stran/smer** ugotovimo s pomočjo **vektorskega produkta**. Točka $A$ je na levi strani vektorja $CD$ (v pozitivni smeri oz. nasprotni smeri urinega kazalca), če je vektorski produkt $CD \\times CA$ pozitiven (na drugi strani bi bil negativen). Posebej pozorni moramo biti na primere, ko se daljici dotikata, kjer je vektorski produkt lahko 0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dffc1cf7-4022-497f-a2f3-0e551dc5ead0",
   "metadata": {},
   "source": [
    "## Površina večkotnika\n",
    "\n",
    "Začnimo s trikotnikom $ABC$. Če imamo podane koordinate oglišč, si lahko izberemo oglišče $A$ za izhodišče in izračunamo polovico absolutne vrednost vektorskega produkta $p = \\frac{1}{2} |AB \\times AC|$.\n",
    "\n",
    "Konveksen večkotnik lahko enostavno razbijemo na trikotnike in uporabimo prejšnji rezultat.\n",
    "\n",
    "<img alt=\"triangulacija (wikipedija)\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/2/23/Convex_polygon_trivial_triangulation.svg/269px-Convex_polygon_trivial_triangulation.svg.png\" width=\"150px\">\n",
    "\n",
    "Na težave naletimo pri večkotnikih, ki niso nujno konveksni. Uporabimo formulo s predznačenimi vsotami trapezov $p=|\\sum _{i=1}^{n}\\frac {1}{2}(y_{i}+y_{i+1})(x_{i}-x_{i+1})|$. Predpostavimo lahko, da se večkotnik v celoti nahaja nad x-osjo (formula deluje tudi brez te predpostavke). Postavimo se nekam na x-os in opazujmo ozek vertikalen stolpec. Vsakič, ko bomo pri obhodu večkotnika prečkali stolpec v desno stran, bomo območje pod njim odšteli, pri prehodu v levo pa prišteli. Marsikaj se bo izničilo in ostala bodo samo območja, ki imajo nad seboj liho število prečkanj (ta se izmenjujejo v levo in desno), kar je ravno notranjost večkotnika. Če naredimo obhod v drugo smer, bo rezultat negativen, po absolutni vrednosti pa enak.\n",
    "\n",
    "<img alt=\"trapezna formula (wikipedija)\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/2/26/Trapez-formel-einf.svg/289px-Trapez-formel-einf.svg.png\">\n",
    "\n",
    "Omenimo še, da deluje enak argument, če si izberemo poljubno izhodišče (npr. (0,0)) in seštevamo predznačene površine trikotnikov, ki jih z izbranim izhodiščem formirajo stranice na robu večkotnika."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "837d2d95-c97a-4d74-bbfb-edecf4924b77",
   "metadata": {},
   "source": [
    "## Vsebovanost točke\n",
    "\n",
    "Začnimo z najenostavnejšim primerom točke $T$, ki se nahaja v trikotniku $ABC$ ali pač ne. Točka se nahaja v trikotniku, če se pri obhodu trikotnika ves čas nahaja na isti strani, kar preverimo z vektorskim produktom. Ali je to pozitivna ali negativna stran, je odvisno od smeri obhoda. Sledeči vektorski produkti morajo imeti enak predznak: $AB \\times AT$, $BC \\times BT$ in $CA \\times CT$. Pozorni moramo biti, kaj problem zahteva v primeru, da se točka nahaja točno na robu trikotnika.\n",
    "\n",
    "Naslednji primer je vsebovanost točke v konveksnem večkotniku. Enostavno ga lahko razbijemo na trikotnike (ki imajo skupno izbrano oglišče) in prevedemo problem na vsebovanost točke v trikotniku. Deluje pa tudi prej omenjeni pristop z lokacijo točke na isti strani obhoda večkotnika.\n",
    "\n",
    "Kako pa rešimo problem za poljuben večkotnik (*point in polygon*), ki ni nujno konveksen? V tem primeru uporabimo tehniko metanja žarka (*ray casting*). Če sledimo poltraku iz točke $T$ v poljubno smer, se ob vsakem križanju z robom večkotnika spremeni lokacija znotraj/zunaj. Če je število križanj liho, je točka znotraj večkotnika, sicer je izven. Pomembna podrobnost je, kaj se zgodi, če žarek seka večkotnik v enem od oglišč. Sprememba je namreč odvisna od sosednjih oglišč. Če več sosednjih oglišč leži na žarku, nas zanima prvo oglišče, ki ne. Če sta obe na isti strani žarka, ni spremembe, sicer pa je. Prikladna izbira smeri je npr. $z=(-1,0)$. Lahko pa se tej komplikaciji izognemo s tako izbiro smeri (naključno), da do tega ne pride.\n",
    "\n",
    "<img alt=\"ray casting (wikipedia)\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/c/c9/RecursiveEvenPolygon.svg/306px-RecursiveEvenPolygon.svg.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b6c9a3c-8d24-4122-9e4b-b513c794ffcd",
   "metadata": {},
   "source": [
    "V spodnji implementaciji bomo predpostavili, da se točka ne nahaja na robu večkotnika. Če to ni res, bi lahko to posebej preverili. Pretvarjali se bomo, da ima točka za $\\epsilon$ večjo y koordinato. To ne spremeni rešitve, vendar poenostavi razmislek, ker so vsa oglišča nad ali pod njo, ne pa na isti višini (oglišča z enako višino bodo obravnavana kot nižja). Problem bi z malo več truda lahko rešili tudi v celih številih, vendar zaradi preglednosti ne bomo dodatno komplicirali."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d641ef5f-a48e-4c04-ba92-b182cd39ff01",
   "metadata": {},
   "outputs": [],
   "source": [
    "#define OPERATOR_SUBTRACT operator- // workaround for a bug of cling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b0db1177-2905-4960-9f1a-268240255a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "PII OPERATOR_SUBTRACT(PII a, PII b) {\n",
    "    return {a.first-b.first, a.second-b.second};\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5f8ecdb2-795c-4b47-bf21-5644da94a3bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "int point_in_polygon(vector<PII> poly, PII t) {\n",
    "    int n=poly.size(), cnt=0;\n",
    "    auto [x,y] = t;\n",
    "    for (int i=0;i<n;i++) {\n",
    "        int j=(i+1)%n;\n",
    "        if ((poly[i].second<=y) != (poly[j].second<=y)) {  // stranica seka vodoravno premico\n",
    "            PII s = poly[j]-poly[i];  // vektor stranice: i -> i+1\n",
    "            PII v = t-poly[i];  // vektor do tocke: i -> t\n",
    "            double k = (double)v.second/s.second;\n",
    "            double xp = poly[i].first + k*s.first;  // presecisce z vodoravno premico\n",
    "            if (xp < x) cnt++;\n",
    "        }\n",
    "    }\n",
    "    return cnt%2;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b0396533-87bf-42ed-8411-afa65b2cdef5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "vector<PII> poly = {{0,0}, {1,1}, {3,1}, {4,2}, {5,1}, {6,2}, {7,0}, {8,1}, {9,0}, {10,1}, {10,3}, {0,3}};\n",
    "cout << point_in_polygon(poly, {9,1}) << endl;\n",
    "cout << point_in_polygon(poly, {9,2}) << endl;\n",
    "cout << point_in_polygon(poly, {6,1}) << endl;\n",
    "cout << point_in_polygon(poly, {5,0}) << endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d708402f-1a2b-4385-a86d-a8899c627cf7",
   "metadata": {},
   "source": [
    "## Konveksna ovojnica\n",
    "\n",
    "Konveksna ovojnica/ogrinjača/lupina (*convex hull*) množice točk v ravnini je najmanjša konveksna množica, ki vsebuje vse podane točke. Običajno nas zanima rob konveksne ovojnice, ki je najkrajša sklenjena črta, ki vsebuje vse točke. Včasih tudi lomljeni črti, ki predstavlja rob konveksne ovojnice, rečemo kar konveksna ovojnica. Predstavljamo si jo lahko kot elastiko, ki se skrči okoli množice točk.\n",
    "\n",
    "<img alt=\"konveksna ovojnica\" src=\"convex-hull.svg\" width=\"250px\">\n",
    "\n",
    "Iščemo ekstremne (robne) točke na robu ovojnice, ki jo definirajo. V primeru več kolinearnih točk na robu, je stvar definicije problema, ali želimo poročati samo oglišča ali tudi točke vzdolž stranic konveksne ovojnice. V nadaljevanju se bomo omejili na primere, kjer ni treh kolinearnih točk.\n",
    "\n",
    "Ogledali si bomo par najbolj klasičnih algoritmov, obstaja pa jih še veliko več. Problem seveda postane težji, če ga rešujemo v treh ali še več dimenzijah."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8f3b57da-7a57-4e3e-a36d-bda8e997eecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "vector<PII> points = {{4,0}, {2,3}, {5,2}, {6,1}, {8,4}, {6,6}, {5,4}, {4,5}, {2,6}, {1,1}, {1,5}, {3,2}};"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a2c3738-8414-44cf-aad0-64d292ead4b0",
   "metadata": {},
   "source": [
    "### Identifikacija stranic\n",
    "\n",
    "Rob konveksne ovojnice je sestavljen iz daljic med pari točk. Če lahko za posamezen par točk oz. daljico med njima ugotovimo, ali je del konveksne ovojnice, lahko zgradimo konveksno ovojnico. Daljica je del roba konveksne ovojnice, če se vse ostale točke nahajajo na isti strani (recimo na levi/pozitivni). Časovna zahtevnost takega postopka je $O(n^3)$, kjer je $n$ število točk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "75879ce5-d366-43b8-979d-ad97538c5cb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "int cross(PII u, PII v) {\n",
    "    return u.first*v.second - u.second*v.first;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7f3b6bb1-41e0-4dec-9144-ddb89487a5dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A D\n",
      "D E\n",
      "E F\n",
      "F I\n",
      "I K\n",
      "J A\n",
      "K J\n"
     ]
    }
   ],
   "source": [
    "int n=points.size();\n",
    "for (int i=0;i<n;i++) {\n",
    "    for (int j=0;j<n;j++) if (i!=j) {  // vektor daljice i-j\n",
    "        PII d=points[j]-points[i];\n",
    "        int ok=1;\n",
    "        for (int k=0;k<n;k++) if (k!=i && k!=j) {\n",
    "            PII v=points[k]-points[i];\n",
    "            if (cross(d,v)<0) ok=0;\n",
    "        }\n",
    "        if (ok) cout << char('A'+i) << \" \" << char('A'+j) << endl;\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "449403ab-ff0a-488f-b565-a4b6d73d7b45",
   "metadata": {},
   "source": [
    "Stranice seveda niso izpisane v vrstnem redu, kot si sledijo na konveksni ovojnici, vendar bi jih lahko uredili, če bi bilo treba."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d9c3bc-c06f-47c8-b0af-2256ca5aa2fd",
   "metadata": {},
   "source": [
    "### Zavijanje darila\n",
    "\n",
    "Pri iskanju konveksne ovojnice smo lahko bolj učinkoviti. Naraven pristop zavijanja darila (*gift wrapping, Jarvis march*) začne z izbiro točke, ki je gotovo del konveksne ovojnice. V ta namen lahko izberemo npr. najbolj levo točko $A$ (najnižjo med najbolj levimi) in raztegnemo ovojni papir navzgor. Papir ovijamo v smeri urinega kazalca dokler se ne dotakne naslednje točke. Postopek ovijanja ponavljamo, dokler ne pridemo do začetne točke.\n",
    "\n",
    "<img alt=\"gift wrapping (wikipedia)\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/d/de/Jarvis_march_convex_hull_algorithm_diagram.svg/266px-Jarvis_march_convex_hull_algorithm_diagram.svg.png\">\n",
    "\n",
    "Naslednjo točko, ki se jo dotakne papir pri ovijanju, lahko poiščemo na več načinov. Ker so med gradnjo konveksne ovojnice vedno vse točke na isti strani zadnje točke A (del neke polravnine skozi A), lahko med njimi poiščemo najbolj levo z uporabo vektorskega produkta $AC \\times AB$ za primerjavo, ali je točka C bolj levo (oz. v nasprotni smeri urinega kazalca) od točke B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "38217878-1508-4850-8843-af6bb71df2b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "VII gift_wrapping(VII points) {\n",
    "    int n=points.size();\n",
    "    PII start=*min_element(points.begin(), points.end());\n",
    "    vector<PII> hull;\n",
    "    PII a=start;\n",
    "    while (a!=start || hull.empty()) {\n",
    "        hull.push_back(a);\n",
    "        PII b = (a!=points[0])?points[0]:points[1];  // katerakoli tocka, ki ni a\n",
    "        for (PII c : points) if (c!=a) {\n",
    "            PII ac=c-a, ab=b-a;  // vektorja AC, AB\n",
    "            if (cross(ab,ac)>0) b=c;\n",
    "        }\n",
    "        a = b;\n",
    "    }\n",
    "    return hull;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "db2920da-095c-48f7-93fd-cd04d56287a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 1) (1, 5) (2, 6) (6, 6) (8, 4) (6, 1) (4, 0) \n"
     ]
    }
   ],
   "source": [
    "auto hull = gift_wrapping(points);\n",
    "print(hull);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39e4a2af-d20d-414e-9619-7024a4803d59",
   "metadata": {},
   "source": [
    "Časovno zahtevnost algoritma lahko analiziramo v odvisnosti od velikosti rezultata (*output-sensitive*) - število točk $h$ na konveksni ovojnici. V tem primeru je časovna zahtevnost $O(hn)$. Če pa velikosti rezultata ne upoštevamo, so lahko v najslabšem primeru vse točke na robu ovojnice, zato je časovna zahtevnost $O(n^2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2057a7e-7210-430a-80e0-d3ff3cf67274",
   "metadata": {},
   "source": [
    "### Grahamov pregled\n",
    "\n",
    "Konveksno ovojnico točk v ravnini lahko poiščemo bolj učinkovito kot v kvadratnem času in sicer z uporabo Grahamovega pregleda (*Graham scan*). Ponovno si izberimo neko ekstremno točko $T$, ki je zagotovo del konveksne ovojnice (npr. najnižjo med najbolj levimi točkami). Uredimo preostale točke glede na kote vektorjev iz točke $T$ (od tistih, ki kažejo navzdol, proti vodoravnim in tistim, ki kažejo navzgor). Naj bo ta urejen seznam točk $P_1, P_2, \\ldots$ Če jih povežemo, dobimo ovojnico, ki vsebuje vse točke (kar v ogliščih), vendar ni konveksna. Vemo tudi, da bo konveksna ovojnica neka podmnožica tega seznama točk. Vrstni red točk je že pravilen, samo izbrati moramo prave.\n",
    "\n",
    "Algoritem gradi konveksno ovojnico postopno z dodajanjem novih točk v izbranem vrstnem redu po kotih. Po vsaki dodani točki popravi konveksnost zgrajene ovojnice, če je nova točka podrla konveksnost z obratom v napačno smer. To naredi z odstranjevanjem točk s konca zgrajene ovojnice, dokler zaključek ovojnice z novo točko ni konveksen.\n",
    "\n",
    "Grahamov pregled pravzaprav gradi vedno večjo konveksno ovojnico z dodajanjem posameznih točk po kotih. V $i$-tem koraku doda točko $P_i$ in iz konveksne ovojnice točk $T, P_1, P_2, \\ldots, P_{i-1}$ izračuna konveksno ovojnico točk $T, P_1, P_2, \\ldots, P_i$.\n",
    "\n",
    "Časovna zahtevnost algoritma je zaradi urejanja $O(n \\log n)$. Preostanek algoritma vključuje dodajanje in odstranjevanje točk z ovojnice, vendar je vsaka točka lahko dodana in odstranjena kvečjemu enkrat. Zato je ta del algoritma linearen v odvisnosti od števila točk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ab64b0df-7a10-4610-adf6-4b56b3f9ce8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "VII graham_scan(VII points) {\n",
    "\tint n=points.size();\n",
    "    PII t=*min_element(points.begin(), points.end());\n",
    "\n",
    "    vector<pair<double,PII>> angles;\n",
    "    for (PII p : points) if (p!=t) {\n",
    "        PII v = p-t;\n",
    "        angles.push_back({atan2(v.second, v.first), p});\n",
    "    }\n",
    "    sort(angles.begin(), angles.end());\n",
    "\n",
    "    vector<PII> hull = {t};  // stack\n",
    "    for (auto [_,c] : angles) {\n",
    "        while (hull.size()>=2) {  // restore convexity\n",
    "            PII a=hull[hull.size()-2], b=hull[hull.size()-1];\n",
    "            PII ab=b-a, ac=c-a;\n",
    "            if (cross(ab,ac)>0) break;\n",
    "            hull.pop_back();\n",
    "        }\n",
    "        hull.push_back(c);\n",
    "    }\n",
    "    return hull;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "67b4ba0a-2966-4eaf-ab9a-90a92aa97402",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 1) (4, 0) (6, 1) (8, 4) (6, 6) (2, 6) (1, 5) \n"
     ]
    }
   ],
   "source": [
    "auto hull2 = graham_scan(points);\n",
    "print(hull2);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "C++17",
   "language": "C++17",
   "name": "xcpp17"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".cpp",
   "mimetype": "text/x-c++src",
   "name": "c++",
   "version": "17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}