Dynamické programování Martin Dolejší, Ondřej Drbohlav 29. října 2013
1
Dynamické programování
Převzato od Ondřeje Drbohlava Dynamické programování je algoritmus pro hledání nejkratší cesty ve vrstveném grafu (layered graph). Vezměme si následující příklad. Jsou celkem čtyři města, která navštívíme v předem známém pořadí. V každém z těchto měst jsou dva obchody. Naším cílem je v každém městě vybrat právě jeden obchod tak, aby byly minimalizovány náklady. Náklady jsou definovány jako cena nákupu ve vybraných obchodech, plus cena přesunu mezi vybranými obchody. Situaci si znázorníme následujícím schématem (grafem):
Náklady na nákup v jednotlivých obchodech jsou zobrazeny v uzlech grafu. Například: v obchodě x ve městě C pořídíme nákup za 5 jednotek. Náklady na přesun mezi obchody v jednotlivých městech jsou na hranách grafu. Například: přesun mezi (x,B) a (y,C) nás bude stát 5 jednotek. Řešením pro tyto konkrétní hodnoty nákladů je následující (cesta vyznačena červeně):
1
Celkové náklady, tj. součet hodnot v uzlech a na hranách cesty, je 27. Řekněme, že nás bude zajímat i to, jakou bychom měli zvolit cestu, pokud bychom chtěli skončit v posledním městě ne v obchodě x, ale v obchodě y. Důvod tohoto zájmu vyplyne z dalších úvah. Řešení je následující:
Náklady na tuto cestu jsou 30 jednotek. Jedná se o neoptimální cestu, ale o nejlepší takovou, která končí v (y,D). Jak úlohy tohoto typu řešit? Je to jednoduché, když si všimneme, že optimální cesty pro N měst nám pomohou k snadnému nalezení cest pro N+1 měst. Řekněme, že do naší úlohy přidáme ještě jedno, páté město E, takže úloha bude popsána následujícím grafem:
Jak najdeme nejlevnější cestu do (x,E) a do (y,E)? Z předchozího víme, že náklady na nejlevnější cesty, které končí v čtvrtém městě, jsou 27 pro (x,D) a 30 pro (y,D). Náš problém s pěti městy si tedy můžeme zredukovat následovně:
Zde 27 a 30 jsou hodnoty nejlevnějších cest končících v (x,D) resp. (y,D). Z jakého obchodu v D půjdeme do (x,E), pokud chceme zachovat podmínku nejmenších nákladů? Můžeme tam jít buď z (x,D) za celkový náklad 27+7+5 = 39, nebo z (y,D) za celkový náklad 30+3+5 = 38. Cesta z (y,D) je levnější a je tedy řešením pro cestu končící v (x,E):
2
Stejnou úvahou pro cestu končící v (y,E) dojdeme k tomu, že k celkovým nižším nákladům vede cesta z (x,D) a její cena je 37:
Tato cesta je zároveň nejlevnější (optimálním) řešením výběru obchodů přes všech pět měst:
1.1
Pseudokód
init: total_cost(x,A) := cost_node(x, A) total_cost(y,A) := cost_node(y, A); for TOWN = {B, C, ... E}: PREVIOUS_TOWN := town preceding TOWN for SHOP = {x,y} find S (which will either be x or y) such that: total_cost(S, PREVIOUS_TOWN) + cost_transition([S, PREVIOUS_TOWN], [SHOP, TOWN]) + cost(SHOP, TOWN) is minimal put the value of this S to: previous_node(SHOP, TOWN) put the value of this minimal cost to: total_cost(SHOP, TOWN)
3
Obrázek 1: Angiografie ledvin. endfor SHOP endfor TOWN Výsledkem tohoto pseudokódu budou minimální náklady (total_cost) pro všechny obchody ve všech městech a pro každý obchod v každém městě (krom prvního) obchod v předchozím městě, odkud vede cesta s minimálními náklady. Abychom zrekonstruovali nejlevnější výběr obchodů přes všechna města, stačí vybrat v posledním městě obchod s nejmenšími celkovými náklady a vytrasovat minimální cestu. Rekonstrukci cesty řeší následující pseudokód: SHOP := argmin_x(total_cost(x, lastTown)) //(shop in last town with minimal total cost) path := [SHOP] for TOWN = {last town, ..., C, B} SHOP := previous_node(SHOP, TOWN) path := [SHOP, path] // append shop to path from the left endfor TOWN
1.2
Automatické označení objektu
Vezměte si obrázek 1. Naším úkolem bude najít v tomto obrázku cévu. Takovou úlohu je možné vyřešit pomocí dynamického programování: • města—jednotlivé řádky obrazu, řazené od shora dolů • obchody—pixely v jednotlivých řádcích • náklady na nákup—náklady na výběr daného pixelu 4
• náklady na přesun mezi obchody—náklady na přesun mezi pixely na sousedních řádcích. Pro jednoduchost zvolíme následující možnost: protože je céva světlá, náklady na výběr pixelu definujeme jako 255-intenzita. Náklady na světlé pixely budou tedy nízké, na tmavé pixely vysoké. Náklady na přesun mezi pixely na sousedních řádcích: 0 v případě, že x-ová souřadnice pixelů se liší maximálně o 1. Nekonečno v případě, že se liší o víc. Dovolené přechody do daného pixelu jsou tedy jen ve třech případech: pixel P přímo nad ním, pixel zleva sousedící s P a pixel zprava sousedící s P. Hledáme tedy takovou posloupnost P bodů [xi , yi ], kde yi = i ∈ 1, 2, 3..., n a |xi − xi+1 | ≤ 1, která minimalizuje i φ(xi , yi ), kde φ(x, y) = fmax − f (x, y) a f je intenzita obrázku. 1.2.1
Implementační poznámky—návrh řešení
Implementaci si rozdělte na dvě funkce—spočítání celkové ceny a přechodů mezi jednotlivými vrstvami, a na vytrasování nejlepší cesty. Výběr pokračování cesty s minimální cenou se pro daný řádek obrazu dá udělat pro všechny pixely řádku najednou. Potřebujeme k tomu funkci min. Hodnoty celkové ceny pro řádek předcházející současnému řádku označíme pro úsporu místa tc. Vytvoříme matici, která bude odpovídat následujícímu schématu: Inf tc(1) tc(2)
tc(1) tc(2) tc(3)
... ... ...
tc(N-2) tc(N-1) tc(N)
tc(N-1) tc(N) Inf
a ke každému řádku ještě navíc přičteme náklady na pixely současného řádku. Když pak aplikujeme funkci min, první výstupní parametr bude odpovídat celkovým nákladům pro současný řádek, v druhém parametru budou čísla 1, 2 a 3, která kódují jestli optimální cesta vede zleva (1), přímo nad daným pixelem (2), nebo zprava (3). Trik s nekonečnem v levém horním a pravém dolním rohu matice souží k tomu, aby nedošlo k výběru levého resp. pravého pixelu pro první resp. poslední pixel řádku (protože levý od prvního neexistuje, stejně tak pravý od posledního).
1.3
Hledání rozhraní
Podívejte se na obr. 2. Naším cílem bude najít pomocí dynamického programování rozhraní mezi tmavou částí obrázku nahoře (sklivec) a světlejší vrstevnatou částí uprostřed (sítnice). Nejdříve bude nutné získat cenovou mapu. Jako první snížíme úroveň šumu mediánovou filtrací im=medfilt2(im,[5,5],’symmetric’); pak vypočítáme gradient intenzity filtrovaného obrázku ∂f (x) ∂f (x) ∇f (x) = , , ∂x1 ∂x2 5
Obrázek 2: Řez sítnicí lidského oka získaný pomocí optical coherence tomography (OCT). Ve středu řezu je vidět macula a subchoroidální výpotek. v diskrétním obraze odhadneme gradient pomocí diferencí. Pro výpočet diferencí použijte konvoluci s jádry [1, −1] a [1, −1]T , spočítat délku gradientu už je snadné. Pokud proměnná gradient obsahuje pro každý pixel velikost gradientu získáme cenovou mapu jako max(gradient(:))-gradient;. Tenhle přepočet je nutný proto, abychom, stejně jako v minulé úloze, mohli hledat minimum. Násobení velikosti gradientu -1 převede hledání maxima na minimum a posunutí můžeme zvolit libovolně. Zde ve zloleno maximum délky gradientu, aby cena nebyla záporná. Formálně hledáme posloupnost P bodů [xi , yi ], kde yi = i ∈ 1, 2, 3..., n a |xi − xi+1 | ≤ 1, která maximalizuje i ||∇f (xi , yi )||2 , kde f je intenzita obrázku. Zkuste si najít nejlevnější cestu zleva doprava a potom nalezněte nejlevnější cestu která začíná v bodě [1,85] a končí v [200,92]. Počáteční bod vnutíme cestě při výpočtu celkové ceny. Stačí v první řadě pixelů nastavit ceny všech bodů kromě počátečního na inf. Koncový bod algoritmu vnutíme při rekonstrukci cesty (rozmyslete si jak).
1.4
Bonus—zobecnění 1
Jistě vám přijde nepraktické, že se námi hledaná křivka může na každé řadě posunout jen o 1px. Bylo by vhodnější tento parametr (jde vlastně o hladkost cesty) volit v závislosti na obrázku. Pokud jste využili náš návrh jak programovat výpočet celkové ceny a předchůdce, stačí pro skok o k pixelů rozšířit schéma pro výpočet celkových nákladů takto:
6
Inf .. . tc(1) .. . tc(k+1)
... .. . ... .. . ...
tc(1) .. . tc(k) .. . tc(2k)
tc(N − 2k) .. . tc(N − k) .. . tc(N)
... .. . ... .. . ...
... .. . ... .. . ...
tc(N − k) .. . tc(N) .. . Inf
Každý z řádků tabulky reprezentuje současný řádek dat posunutý o -k až k pozic. Dohromady s nulovým posunem dostaneme tedy tabulku o k + 1 řádcích. V rozšíření si můžete všimnout další nedokonalosti. Ať je skok mezi pixely libovolně velký neprojeví se jeho velikost v ceně cesty. Náprava je snadná, stačí k matici vzdáleností přičíst následující matici vynásobenou α: k .. . 0 .. . k
. . .k .. . ... .. . ...
... .. . 0 .. . k
k .. . ... .. . ...
... .. . 0 .. . k
... .. . ... .. . ...
k .. . 0 .. . k
Cena se pak bude skládat z ceny pixelů a ceny hran mezi pixely. Cena pixelů odpovídá datům, cena hran délce cesty. Zvolíme-li parametry k a α dost velké, neuplatní se už omezení počtu sousedů ve výsledku (maximální skok mezi vrstvami bude příliš dlouhý na to aby byl vybrán), k jen omezí velikost matice ve které hledáme maxima. Parametrem hladkosti cesty se stává α.
1.5
Bonus—zobecnění 2
Jistě vám přijde nepraktické, hledat hrany od kraje do kraje obrázku. Nepraktické je také hledat cestu shora dolů či zleva doprava. Logickým zobecněním je hledání hran mezi dvěma libovolnými body (stále ještě ve vrstveném grafu). Zde se nabízejí dvě řešení. 1. Otočit obrázek tak, aby měl počáteční a koncový bod stejnou x-ovou popřípadě y-ovou souřadnici, následně obrázek oříznout tak aby byly oba body na krajích, najít nejlevnější cestu a výsledek transformovat zpět. 2. Oříznout obrázek tak, aby byl počáteční i koncový bod na kraji obrazu ve smyslu osy x, spočítat nejlevnější cestu, oříznout původní obrázek tak aby byl počáteční i koncový bod na kraji obrazu ve smyslu osy y, spočítat nejlevnější cestu, a porovnáním obou cest vybrat tu s nižšími náklady.
1.6
Bonus—Dijkstrův algoritmus
Jistě vám přijde omezující hledat body vždy jen v jednom směru bez možnosti návratu. Toto omezení je v dynamickém programování nepřekonatelné, jeho zavedením však získáme na rychlosti. Optimální algoritmus na hledání cesty v
7
obecném grafu který neobsahuje hrany záporných délek navrhl Dijkstra. Algoritmus prochází graf od počátečního bodu a snaží se pokračovat k sousedům vrcholu s nejnižší doposud zjištěnou cenou. Více například zde.
2
Zadání 1. Stáhněte si potřebná data 2. Naprogramujte funkci [cost,prec]=DPcost(image,direction,mstep,start) která vypočítá celkovou cenu nejlevnější cesty do každého pixelu a ukazatel na předchozí pixel. Vstupní parametry: • image [n × m] vstupní obraz odpovídající cenám jednotlivých pixelů. • direction string ’lr’ pro hledání cesty z levého okraje obrázku do pravého (i opačně) nebo ’tb’ pro cestu shora dolů. • mstep [1 × 1] BONUS počet pixelů o který smí cesta zahnout v každé z vrstev (standardně 1, pokud nebudete řešit bonusové úlohy, tento parametr nebude mít vliv na běh funkce, nechte ho však v hlavičce funkce). • start [1 × 1] číslo řádku či sloupce (podle parametru direction) kde má cesta začít není-li parametr zadán hledáme bod s minimálně cenou. Výstupní parametry: • cost [n × m] celková cena každého z pixelů. • prec [n × m] matice obsahující řádkovou či sloupcovou souřadnici (podle parametru direction) bodu v předchozí vrstvě, odkud vede do pixelu nejlevnější cesta. (2b.) 3. Naprogramujte funkci [track,tcost]=backtrackDP(cost,prec,direction,endp) která na základě celkové ceny a tabulky předchůdců vrátí nejlevnější cestu v grafu (obrázku). Vstupní parametry: • cost [n × m] celková cena každého z pixelů. • prec [n × m] matice předchůdců. • direction string ’lr’ pro hledání cesty z levého okraje obrázku do pravého (i opačně) nebo ’tb’ pro cestu shora dolů. • endp [1 × 1] číslo řádku či sloupce (podle parametru direction) kde má cesta končit není-li parametr zadán hledáme bod s minimální cenou.
8
100 200 300 400 500 600 700 100
200
300
400
500
600
700
Obrázek 3: Fotografie očního pozadí (ang. fundus image) s vyznačenou nejlevnější cestou podle bodu 6 zadání. Výstupní parametry: • track [n×1] vektor řádkových či sloupcových souřadnic (podle parametru direction) bodů ležících na nejkratší cestě (body seřaďte ve smyslu rostoucí osy x, resp. y, první bod cesty bude mít souřadnice [1,track(1)] resp. [track(1),1]. • tcost [1 × 1] celková cena nalezené cesty. (2b.) 4. Pomocí DP najděte cévu v obrázku 1. Nalezenou cestu zakreslete do vstupniho obrázku. Výsledek vložte do reportu. (0,5b.) 5. Pomocí DP najděte rozhraní v obrázku 2, počáteční bod [1,85], koncový bod [200,92]. Před výpočtem gradientu snižte úroveň šumu mediánovou filtrací im=medfilt2(im,[5,5],’symmetric’);. Nalezenou cestu zakreslete do vstupního obrázku. Výsledek vložte do reportu. (0,5b.) 6. BONUS Modifikujte funkci DPcost tak, aby využívala parametr mstep parametr α volte 1/10*mean(image(:)). Výsledek pro obrázek fundus, mstep=10 počátek v bodě [1,305] konec [768,220] je na obrázku 3(+1b.) 7. BONUS S využitím předchozích funkcí napište program, který dokáže pomocí DP najít nejkratší cestu mezi dvěma libovolnými body obrázku (jednu z variant). Počáteční a koncový bod nechte uživatele zadat kliknutím do obrazu. (+1b.) 8. BONUS Dijkstrovým algoritmem nalezněte nejlevnější cestu mezi dvěmi libovolnými body obrazu. Hrany grafu nechť vedou vždy mezi sousedícími pixely ve 4-okolí. Body nechte uživatele zadat myší. (+5b)
9
