Matematické Fórum

Fra · 27. 01. 2009 21:08

Ahoj všem.

Ze školy jsem již nějaký ten pátek, takže hledám pomoc.
Řeším takový problém. Mým koníčkem je astronomie a astrofotografie (bilek.astronomy.cz). Momentálně pracuji na jednoduchém software, který má za úkol pomoci se zaostřením při astrofotografii. Princip je jednoduchý- nafotím sérii tzv. pokusných expozic (každá expozice s jinou polohou zaostření) a kontrolou jednotlivých expozic vybrat tu nejostřejší a podle ní zpětně nastavit zaostření. Neastrofotograf to asi nepochopí, takže to nebudu dále rozvádět...
Pro vyhodnocení zaostření se u profi software používá funkce FWHM (Full Width at Half Maximum). Předpokládám, že matematikům nemusím vysvětlovat, o co jde.
Na internetu jsem našel spoustu informací o FWHM, ale všechny se týkaly jen obecného popisu, definice, ale nenašel jsem nic o tom, jak tuto hodnotu stanovit z dat.
Mám takovou představu. Hvězda se na snímku většinou rozprostírá na několika pixelech. Udělal bych řez touto hvězdou, čímž bych získal hodnoty jednotlivých pixelů napříč hvězdou. Např. hvězda zabírá na snímku pět pixelů, hodnoty jednotlivých pixelů (jasnost, ADU, ...) např. 2, 8, 22,10,3.
Vynesením dat do grafu-křivky lze stanovit kvalitu zaostření- čím je křivka ostřejší a užší, tím je i hvězda ostřejší. Já to ale nechci kontrolovat pohledem na křivku, já to potřebuji stanovit matematicky. Na webu jsem našel vztah mezi FWHM a standardní deviací: FWHM= ~2.35*st_deviace. Výsledek ale někdy výrazně rozdílný od výsledku, který získám kontrolní "grafickou metodou"- na papíře vynesu data do grafu, od oka proložím křivkou a FWHM změřím pravítkem :). Je to sice hrubá metoda, ale pro kontrolu dostatečně přesná.
Máte někdo nápad, jak na to?
Franta

Kondr · 27. 01. 2009 22:08

Ten odhad, že FWHM= ~2.35*st_deviace je založen na tom, že ta křivka je podobná normálnímu rozdělení. Možná by šlo lepších výsledků dosáhnout, kdyby jsi ty hodnoty proložil nějakou funkcí: http://cs.wikipedia.org/wiki/Lagrangeova_interpolace a zkusil numericky najít body, v nichž tento polynom protíná polovinu maxima. Případně to proložit lomenou čarou -- pak by stačilo spočítat body, které jsou nad polovinou maxima a pro ty těsně okolo určit, z kolika procent jsou nad polovinou maxima, a to pomocí jednoduchého dělení a násobení. Ale on se určitě přidá někdo, kdo toho ví o numerických metodác daleko víc jak já.

Lishaak · 27. 01. 2009 23:13

Ja tedy te astrofotograficke terminologii vubec nerozumim, ale pokud prokladam krivku (takzvana interpolace) lagrangeovym polynomem, tak ten polynom sice prochazi nemarenymi hodnotami, ale obcas se stane, ze se v okrajovych bodech neprijemne rozkmita, takze potom to v intervelech mezi temi namerenymi hodnotami muze docela uletet a vubec se to nepodoba krivce, kterou by clovek nakreslil od oka rukou. Treba to zrovna v tomto pripade bude fungovat. Chce to si s tim trosku zaexperimentovat.
Kdyby to ale nefungovalo, da se pouzit interpolace napriklad pomoci takzvanych uniformnich kubickych B-spline krivek, ktere se nerozkmitaji a velmi se podobaji tomu, co by clovek nakreslil od oka rukou. Tyto interpolace (lagrange, mozna b-spline a pak jeste nejake zalozene na goniometrickych funkcich) umi krasne pocitat a zobrazovat program jmenem Graph (staci zadat do googlu). Takze je mozne si na konkretnich merenich nechat spocitat danou interpolaci a podivat se, jak to vypada.

Fra · 28. 01. 2009 17:52

Kondr napsal(a):
Ten odhad, že FWHM= ~2.35*st_deviace je založen na tom, že ta křivka je podobná normálnímu rozdělení. Možná by šlo lepších výsledků dosáhnout, kdyby jsi ty hodnoty proložil nějakou funkcí: http://cs.wikipedia.org/wiki/Lagrangeova_interpolace a zkusil numericky najít body, v nichž tento polynom protíná polovinu maxima. Případně to proložit lomenou čarou -- pak by stačilo spočítat body, které jsou nad polovinou maxima a pro ty těsně okolo určit, z kolika procent jsou nad polovinou maxima, a to pomocí jednoduchého dělení a násobení. Ale on se určitě přidá někdo, kdo toho ví o numerických metodác daleko víc jak já.

No, tak to jsi mi moc nepomohl :) . Mám jen maturitu, sice z matiky, ale ze školy jsem dvacet let, a od té doby jsem s matematikou do styku moc nepřišel...

Potřebuji poradit, JAK proložit hodnoty funkcí- polopatě.
Pokud vím, tak se v tomto případě (FWHM) používá Gaussova křivka, ale s tím si bez cizí pomoci neporadím. Když čtu tu tvou odpověď, tak si myslím, že by to takto jistě šlo, jen nevím, jak to prakticky udělat (dvacet let je proklatě dlouhá doba...)

Franta

Lishaak · 28. 01. 2009 18:43

Az ted jsem si poradne precetl to FWHM a je to vlastne docela jednoducha vec. Je treba ale rict, jak slozity vypocet jeste povazujete za unosny. Pokud to chcete pocitat rucne, tak se budem asi muset spokojit s jednodussimi metodami (nakonec by byla stejne nejlepsi ta vase metoda od ruky). Pokud to chete programovat, moznosti jsou takrka neomezene.

Vas problem spociva v tom nalezt takovou interpolaci, ktera si nejvic podoba vasi rucni interpolaci. Me osobne se nejvic libi interpolace pomoci B-spline krivek, protoze to jsou polynomy a s temi se hezky pocita. Nejsem si jist, jestli ten Graph umi tu B-spline interpolaci, nemuzu to ani overit, protoze to je program pro Windows, ktere ja napouzivam. Jsem ochoten vam poradit napriklad jak naprogramovat tu b-spline interpolaci a jak z toho spoitat to FWHM. Ale chtelo by to nejdriv najit nejaky software, ktery umi tu krivku ukazat. Mrknu se po nem.

Ted me napada, ze ani rucne by se to FWHM z tech B-spline nepocitalo az zas tak slozite.

Lishaak · 28. 01. 2009 18:49

Tak jsem to nasel. Jak vypada interpolace pomoci kubickych spline krivek si muzete vyzkouset tady. Dokonce si tam muzete vyzkouset i tu Kondrem navrhovanou Lagrangeovu interpolaci.

Fra · 28. 01. 2009 19:20

Lishaak napsal(a):
...Jsem ochoten vam poradit napriklad jak naprogramovat tu b-spline interpolaci a jak z toho spoitat to FWHM...

O radu samozřejmě stojím, jinak bych sem nepsal, takže uvítám jakoukoli radu. Pro mě je ideální rada typu "příklad", prostě ukázat postup řešení na ukázkovém příkladu.

Např. mám soustavu pěti bodů (hodnoty pixelů např. v ADU)
x-ová souřadnice s přírůstkem po jedné (pixely jsou vedle sebe), y-ová souřadnice jsou hodnoty pixelů v ADU):
(1,2), (2,7), (3,13), (4,8), (5,4) - pro zjednodušení jsem použil nízké a celé hodnoty

Tak a můžeme začít s výukou :)))

F.

P.S. Nejsem si jistý, zda je to možné řešit tou b-spline, pokud vím, tak se to vždy řeší jakýmsi Gaussovským modelem, ale zkusit to můžemě...

Lishaak · 28. 01. 2009 19:25

Fra napsal(a):
P.S. Nejsem si jistý, zda je to možné řešit tou b-spline, pokud vím, tak se to vždy řeší jakýmsi Gaussovským modelem, ale zkusit to můžemě...

Aha. A nemate nejaky odkaz na material, ze ktereho jste toto vycetl? Pokud na to primo existuje konkretni metoda, tak tim lip.

Pavel Brožek · 28. 01. 2009 19:32

↑ Fra:

Pokud se to má řešit jako Gaussova křivka, pak opravdu bude nutné použít ten vztah $\text{FWHM}=2\sqrt{2\ln2}\sigma\approx2,35\sigma$ . Bude tedy nutné určit $\sigma$ proložením Gaussovy křivky těmi body. Jestliže Vám ale vychází výsledek výrazně jiný než grafickou metodou, pak asi Gaussova křivka nebude ta správná závislost. Jestli přesto chcete prokládat Gaussovou křivkou, pak se můžu podívat, jak se $\sigma$ z dat spočte.

Fra · 28. 01. 2009 19:42

Například toto mi poradil jede kolega z oboru:

"Nejsem odborník, ale podle mě standardní procetura je fitování nějaké zvolené hladké funkce na průběh jasnosti hvězdy. Logický první typ je gaussovka, ale možná se používá i něco složitějšího. Nevím, jak moc ti říká pojem "fitovat", zkrátka to znamená najít takové parametry funkce (poloha středu, výška maxima, "roztaženost"), aby odchylka od hodnot v jednotlivých pixelech (typicky součet druhých mocnin odchylek v každém pixelu) byla co nejmenší. Stanovit pak FWHM je už snadné, pro funkci s danými parametry je známo, kde to je."

To je sice dobrá odpověď, ale mě nic neříkající, já to musím vidět na příkladu (jak jsem psal, nejsem matematik)...

To BrozekP: určitě se podívejte, sem s tím!

Franta

Lishaak · 28. 01. 2009 20:02

↑ Fra:
Jasne. Tomu, co vas kolega nazyva "fitovanim" se rika aproximace. Mezi aproximaci a interpolaci je rozdil v tom, ze kdyz aproximuji, hledam krivku urcitych vlastnosti, ktera se od tech namerenych hodnot lisis co nejmin, naopak interpolace je, ze krivka prochazi primo temi namerenymi hodnotami. Ted je otazka, co je ve vasem pripade vyhodnejsi. To se obecne neda rict.

Pokud budeme pouzivat Gaussovu krivku, tak pujde o aproximaci, protoze tezko muzeme prepokladat, ze se nejake namerene hodnoty strefi do te krivky uplne presne. Me by ale zajimalo, kde jste prisel ke slovum "Gaussuv model" a zda jste narazil na nejaky zdroj, kde se popisuje, jak se to vyuzije pri pocitani toho FWHM.

Pavel Brožek

↑ Fra:

Pokud ty body budu značit $(x_i,y_i),\,i=1,\ldots n$ , pak určím

$\mu=\frac{\sum_{i=1}^nx_iy_i}{\sum_{i=1}^ny_i}$

$\sigma^2={\frac{\sum_{i=1}^ny_i(x_i-\mu)^2}{\sum_{i=1}^ny_i}}$

Dá se to spočítat jednoduše i pomocí excelu.

Nejsem si těmi vzorci ale příliš jist, spíš jsem si to tak logicky odvodil, pokusím se ještě podívat po internetu, jestli by k tomu někde nebyl nějaký zdroj.

Pro (1,2), (2,7), (3,13), (4,8), (5,4) vyjde $\sigma=1,061$ ,
pro (1,2), (2,8), (3,22), (4,10), (5,3) z prvního příspěvku $\sigma=0,915$ .

Fra · 28. 01. 2009 21:05 — Editoval Fra (28. 01. 2009 21:07)

To Lishaak: Za posledních pár dní jsem prošel stovky stránek, kde bylo něco o FWHM, kde bylo něco o tom Gaussově modelu už opravdu nevím. Tím Gaussovým modelem určitě mysleli proložení bodů gassovou křivkou. Myslím si ale, že se jistě jedná o aproximaci a ne o interpolaci. Když si představím grafy FWHM, co jsem na webu viděl, tak určitě neprocházely přes naměřená data-body...

Zde je třeba vidět profil měřené hvězdy a něco málo o FWHM- používají "a two dimensions gauss function" (je to manuál k programu Iris):
http://www.astrosurf.com/buil/iris/tuto … c38_us.htm

To BrozekP: to již vypadá hezky, tomu již bych možná byl schopný porozumět, nicméně jak jsem psal, jsem ze školy dvacet let, a tyto symboly mi již mnoho neříkají, vidím jen, že tam jsou nějaké sumy, ... :) . Abych pochopil ty vztahy, musel bych prohrabat sklep, najít knížky, prostudovat, ...
Uvítal bych názorný postup, ne jen vzorec a výsledek. Vím, že toho chci asi moc, moc se za to omlouvám, ale budu rád, pokud mi to názorně předvedete.

Od onoho kolegy z oboru jsem dostal další odpověď:

"Pro jednoduchost uvažujme jen jeden rozměr (tedy jakože máš jenom jednu řádku obrázku), zobecnění na skutečný dvojrozměrný obrázek je myšlenkově snadné, technicky trochu pracnější.

Signál, který tvůj dalekohled spolu s aatmosférou vytvoří z bodové hvězdy, je nějaká funkce, terou pro svoje účely asi klidně můžeš za gausovku, jako v tom článku na wikipedii (http://en.wikipedia.org/wiki/Full_width_at_half_maximum). Taková gaussovka je charakterizována třemi parametry, ve vzorci v článku je to x-nula, sigma, a pak to celé ještě může být přenásobeno nějakou konstantou N, která udává, jak jasná hvězda vlastně byla.

Detektor (CCD, cokoliv) ti ovšem neřekne tu funkci, ale pouze její hodnoty v několika bodech - pixelech. Pro jednoduchost to můžeš brát tak, že hodnota každého pixelu odpovídá hodnotě funkce v jeho prostředku. Takže máž třeba 10 bodů na té funkci (kam padla hvězda, okolo je tma a šum) a teď hledáš, s jakou kombinacé x-nula, sigma a N získáš funkci, která co nejlépe projde těmito body. Když se ti to podaří, tak uděláš 2.35 krat sigma a máš vytouženou FWHM.

Jednoduchou metodou jak tohohole docílit je metoda nejmenších čtverců (least squares, najdi si o tom něco googlem). Naivně funguje tak, jak jsem o to už zavadil, pro každou kombinaci x-nula, sigma a N se spočítá součet druhých mocnin odchylek takové funkce od těch bodů, co si naměřil, a hledá se kombinaci, pro níž je tohle minimum.
"

Matematik toto jistě hravě pochopí, ale jak jsem již několikrát psal, potřebuji názorný příklad na konkrétních číslech...

Franta

Pavel Brožek · 28. 01. 2009 21:22

↑ Fra:

Tak abych nemusel tolik vypisovat, předpokládejme pouze 3 naměřené pixely (1,1) (2,13) (3,5).

$\mu=\frac{\sum_{i=1}^nx_iy_i}{\sum_{i=1}^ny_i}=\frac{x_1y_1+x_2y_2+x_3y_3}{y_1+y_2+y_3}=\frac{1\cdot1+2\cdot13+3\cdot5}{1+13+5}=\frac{1+26+15}{19}=\frac{42}{19}$

$\sigma^2=\frac{\sum_{i=1}^ny_i(x_i-\mu)^2}{\sum_{i=1}^ny_i}=\frac{1\cdot(1-\frac{42}{19})^2+13\cdot(2-\frac{42}{19})^2+5\cdot(3-\frac{42}{19})^2}{1+13+5}\approx0,271$

$\sigma=\sqrt{\sigma^2}=\sqrt{0,271}\approx0,521$

Fra · 28. 01. 2009 21:32

BrozekP napsal(a):
↑ Fra:

Tak abych nemusel tolik vypisovat, předpokládejme pouze 3 naměřené pixely (1,1) (2,13) (3,5).

$\mu=\frac{\sum_{i=1}^nx_iy_i}{\sum_{i=1}^ny_i}=\frac{x_1y_1+x_2y_2+x_3y_3}{y_1+y_2+y_3}=\frac{1\cdot1+2\cdot13+3\cdot5}{1+13+5}=\frac{1+26+15}{19}=\frac{42}{19}$

$\sigma^2=\frac{\sum_{i=1}^ny_i(x_i-\mu)^2}{\sum_{i=1}^ny_i}=\frac{1\cdot(1-\frac{42}{19})^2+13\cdot(2-\frac{42}{19})^2+5\cdot(3-\frac{42}{19})^2}{1+13+5}\approx0,271$

$\sigma=\sqrt{\sigma^2}=\sqrt{0,271}\approx0,521$

Tak toto je již naprosto jasné, dík za to!
Jen prosím o kontrolu tech vztahů (jak jste psal výše), ať mám jistotu

Dík moc.
Přivítám ještě další nápady (metody), jak to řešit...
F.

Fra · 28. 01. 2009 21:34

Všiml jsem si, že již nejsem "Zelenáč", ale "Matemetik" :)))))))))))))))))))))))))))))))))))))
(pod nickem)
F.

Lishaak · 28. 01. 2009 21:39

Pokud byste to chtel resit tou metodou nejmensich ctvercu, tak vzhledem k tomu, ze mame nelinearni funkci (gaussovu) tak to rucne neupocitate. Existuji ale programy (napriklad program R, zdarma dostupny pro Linux i Windows), ktery umi nejmensi ctverce pro jakoukoliv funkci, kterou mu zadate. Pokud byste si tedy stahl tento program, nemusite pocitat vubec nic. Staci zadat namerene hodnoty a program to spocita za vas. Mate-li zajem, muzu napsat skript pro tento program, takze staci jen do souboru napsat nemerene hodnoty a spustit tento skript.

Fra · 28. 01. 2009 21:48 — Editoval Fra (28. 01. 2009 21:55)

Lishaak napsal(a):
Pokud byste to chtel resit tou metodou nejmensich ctvercu, tak vzhledem k tomu, ze mame nelinearni funkci (gaussovu) tak to rucne neupocitate. Existuji ale programy (napriklad program R, zdarma dostupny pro Linux i Windows), ktery umi nejmensi ctverce pro jakoukoliv funkci, kterou mu zadate. Pokud byste si tedy stahl tento program, nemusite pocitat vubec nic. Staci zadat namerene hodnoty a program to spocita za vas. Mate-li zajem, muzu napsat skript pro tento program, takze staci jen do souboru napsat nemerene hodnoty a spustit tento skript.

Já ale potřebuji tu funkci implementovat do mé aplikace. Jedině, že by můj program zapsal naměřená data do souboru, na pozadí spustil onen program a nějak z něj nějak získat výsledek...
(prolétl jsem manuál k programu a nikde jsem nenašel informaci, že by program automaticky ukládal výsledek, tak nevím, zda bych to využil)

Jinak, nevím, jestli si rozumíme, já nechci nic ručně počítat, já to chci svěřit programu- aplikaci, kterou vyvíjím.

F.

Lishaak · 28. 01. 2009 23:24

Tak v tom pripade mate v podstate dve moznosti. Bud se spokojit s jednoduchym resenim ala BrozekP, ktere mozna bude vasim potrebam stacit. A kdyz ne, chce to vyzkouset sofistikovanejsi metody, kterym ale budete muset porozumet, nez se je budete snazit implementovat. Metoda nejmensich ctvercu pro gaussovu funkci urcite nebude zas tak tezke implementovat, ale da to asi vetci praci nez par jednoduchych sum.

Predpokladam, ze pred tim, nez se rozhodnete neco programovat, chcete prozkoumat ruzne moznosti a podivat se, jake davaji vysledky. Takze si muzete programem R nechat aproximovat nejaky objem pokusnych dat a kdyz to bude vychazet podle vasich predstav tak se snazit to sam naprogramovat. Stejne tak muzete vyzkouset metody od BrozekP, nebo nakonec i ty interpolacni polynomy jako je Lagrange nebo b-spline a vybrat si, co se vam libi nejvic.

Pavel Brožek · 28. 01. 2009 23:27

Zkoumal jsem ty mé vzorce. Fungují pouze přibližně a jen díky tomu, že rozložení x je rovnoměrné (tj. intervaly mezi pixely jsou vždy 1 pixel). Nespoléhal bych se na ně.

Pokud můžete spoléhat na cizí program, který umí fitovat, tak např. gnuplot by se myslím dal volat z jiného programu.

Matematické Fórum

#1 27. 01. 2009 21:08

fwhm

#2 27. 01. 2009 22:08

Re: fwhm

#3 27. 01. 2009 23:13

Re: fwhm

#4 28. 01. 2009 17:52

Re: fwhm

Kondr napsal(a):

#5 28. 01. 2009 18:43

Re: fwhm

#6 28. 01. 2009 18:49

Re: fwhm

#7 28. 01. 2009 19:20

Re: fwhm

Lishaak napsal(a):

#8 28. 01. 2009 19:25

Re: fwhm

Fra napsal(a):

#9 28. 01. 2009 19:32

Re: fwhm

#10 28. 01. 2009 19:42

Re: fwhm

#11 28. 01. 2009 20:02

Re: fwhm

#12 28. 01. 2009 20:30 — Editoval BrozekP (28. 01. 2009 20:33)

Re: fwhm

#13 28. 01. 2009 21:05 — Editoval Fra (28. 01. 2009 21:07)

Re: fwhm

#14 28. 01. 2009 21:22

Re: fwhm

#15 28. 01. 2009 21:32

Re: fwhm

BrozekP napsal(a):

#16 28. 01. 2009 21:34

Re: fwhm

#17 28. 01. 2009 21:39

Re: fwhm

#18 28. 01. 2009 21:48 — Editoval Fra (28. 01. 2009 21:55)

Re: fwhm

Lishaak napsal(a):

#19 28. 01. 2009 23:24

Re: fwhm

#20 28. 01. 2009 23:27

Re: fwhm

Zápatí