Commit 34357ff5 authored by Jakob Knollmueller's avatar Jakob Knollmueller

cleanup

parent 8c9c6556
......@@ -17,70 +17,66 @@
# Starblade is being developed at the Max-Planck-Institut fuer Astrophysik
import numpy as np
from astropy.io import fits
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
import nifty4 as ift
from scipy.stats import invgamma
import starblade as sb
if __name__ == '__main__':
#specifying location of the input file:
path = 'data/hst_05195_01_wfpc2_f702w_pc_sci.fits'
path = 'data/frame-i-004874-3-0692.fits'
path = 'data/frame-i-007812-6-0100.fits'
path = 'data/frame-g-002821-6-0141.fits'
path = 'data/frame-g-002821-6-0141.fits'
path = 'data/frame-i-000752-1-0432.fits'
# path = 'data/frame-i-004858-1-0480.fits'
# data = fits.open(path)[1].data
# data = fits.open(path)[0].data#[750:,1000:]
xx = 250
yy = int(xx /0.75)
data = fits.open(path)[0].data[500:500+xx,1200:1200+yy]
data_unmod = fits.open(path)[0].data[500:500+xx,1200:1200+yy]
sex = fits.open('data/check.fits')[0].data[500:500+xx,1200:1200+yy]
data -= data.min() - 0.001
# data = 1.-plt.imread('data/sdss.png').T[0]
# data = fits.open(path)[1].data
def generate_mock_data():
s_space = ift.RGSpace([128,128])
h_space = s_space.get_default_codomain()
FFT = ift.FFTOperator(h_space)
p_spec = lambda k: (1./(1+k)**4)
S = ift.create_power_operator(h_space, power_spectrum=p_spec)
sh = S.draw_sample()
s = FFT(sh)
q = 1e-3
alpha = 1.5
brightness = ift.log(ift.Field(s_space, val=invgamma.rvs(alpha-1., scale=q, size=s_space.shape)))
data = data.clip(min=0.0001)
u = brightness
d = ift.exp(s) + ift.exp(u)
return d.val
hdu = fits.PrimaryHDU(data)
hdul = fits.HDUList([hdu])
hdul.writeto('new1.fits')
data = np.ndarray.astype(data, float)
vmin = np.log(data.min()+0.2)
vmax = np.log(data.max())*0.3
plt.gray()
lin_max = 2.
lin_min=0.1
plt.imsave('log_data.png', np.log(data),vmin=vmin,vmax=vmax)
plt.imsave('data.png', (data), vmax = lin_max,vmin =lin_min)
plt.imsave('sex.png', sex)#, vmax = lin_max, vmin=lin_min)
plt.imsave('log_sex.png', np.log(sex),)#, vmax = lin_max, vmin=lin_min)
plt.imsave('point_sex.png', (data_unmod - sex), vmax = lin_max,vmin =lin_min)
alpha = 1.4
Starblade = sb.build_starblade(data, alpha=alpha, cg_steps=100, newton_steps=10)
for i in range(1000):
Starblade = sb.starblade_iteration(Starblade, samples=5)
#plotting on logarithmic scale
plt.imsave('log_diffuse_component.png', Starblade.s.val, vmin=vmin, vmax=vmax)
plt.imsave('diffuse_component.png', np.exp(Starblade.s.val), vmin=lin_min, vmax=lin_max)
if __name__ == '__main__':
np.random.seed(42)
data = generate_mock_data()
myStarblade = sb.build_starblade(data=data, alpha=1.5, newton_steps=200, cg_steps=5, q=1e-3)
plt.imsave('log_pointlike_component.png', Starblade.u.val, vmin=vmin, vmax=vmax)
plt.imsave('pointlike_component.png', np.exp(Starblade.u.val), vmin=lin_min, vmax=lin_max)
for i in range(3): # not fully converged after 3 steps.
myStarblade = sb.starblade_iteration(myStarblade, samples=5, cg_steps=10,
newton_steps=100, sampling_steps=1000)
plt.figure()
k_lenghts = Starblade.power_spectrum.domain[0].k_lengths
plt.plot(k_lenghts, Starblade.power_spectrum.val)
plt.title('power spectrum')
plt.yscale('log')
plt.xscale('log')
plt.ylabel('power')
plt.xlabel('harmonic mode')
plt.savefig('power_spectrum.png')
### PLOTS ###
vmin = data.min()
vmax = data.max()*0.1
size = 15
# plot data
plt.figure()
plt.imshow(data, norm=LogNorm(), vmin=vmin, vmax=vmax)
cbar = plt.colorbar()
cbar.set_label('intensity', size=size)
plt.axis('off')
plt.title('data', size=size)
plt.savefig('data.pdf')
plt.close()
# plot diffuse component
plt.figure()
plt.imshow(np.exp(myStarblade.s.val), norm=LogNorm(), vmin=vmin, vmax=vmax)
cbar = plt.colorbar()
cbar.set_label('intensity', size=size)
plt.axis('off')
plt.title('diffuse component', size=size)
plt.savefig('diffuse.pdf')
plt.close()
# plot point-like component
plt.figure()
plt.imshow(np.exp(myStarblade.u.val), norm=LogNorm(), vmin=vmin, vmax=vmax)
cbar = plt.colorbar()
cbar.set_label('intensity', size=size)
plt.axis('off')
plt.title('point-like component', size=size)
plt.savefig('points.pdf')
plt.close()
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2017-2018 Max-Planck-Society
# Author: Jakob Knollmueller
#
# Starblade is being developed at the Max-Planck-Institut fuer Astrophysik
import numpy as np
from matplotlib import pyplot as plt
import starblade as sb
if __name__ == '__main__':
# data = plt.imread('10Keso1242a.tif')
data = plt.imread('data/galaxy.jpg')
data = data.astype(float)
data = data.clip(0.0001)
alpha = 1.3
MultiStarblade = sb.build_multi_starblade(data, alpha=alpha)
for i in range(50):
MultiStarblade = sb.multi_starblade_iteration(MultiStarblade, processes=1)
#plotting a three channel RGB image in each iteration
diffuse = np.empty_like(data)
point = np.empty_like(data)
for i in range(len(MultiStarblade)):
diffuse[...,i] = np.exp(MultiStarblade[i].s.val)
point[...,i] = np.exp(MultiStarblade[i].u.val)
plt.imsave('rgb_diffuse.jpg',diffuse/255.)
plt.imsave('rgb_point.jpg',point/255.)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2017-2018 Max-Planck-Society
# Author: Jakob Knollmueller
#
# Starblade is being developed at the Max-Planck-Institut fuer Astrophysik
import numpy as np
from astropy.io import fits
from matplotlib import pyplot as plt
import scipy.cluster.vq as sp
import starblade as sb
if __name__ == '__main__':
#specifying location of the input file:
path = '../data/hst_05195_01_wfpc2_f702w_pc_sci.fits'
path = '../data/frame-i-004874-3-0692.fits'
# path = '../data/frame-i-007812-6-0100.fits'
# path = '../data/frame-g-002821-6-0141.fits'
path = '../data/frame-i-000752-1-0432.fits'
# path = '../data/frame-u-000752-1-0432.fits'
# path = '../data/frame-i-006174-2-0094.fits'
# path = 'data/frame-i-004858-1-0480.fits'
# data = fits.open(path)[1].data
# data = fits.open(path)[0].data#[750:,1000:]
##docker run --rm -i -t --name sex -v ~/Projects/starblade/demos/sextractor:/work chbrandt/sextractor my_data.fits -c default.se
xx = 450
yy = int(xx /0.75)
x0 = 400#1050#770#450
y0 = 1000#1550#200#1150
data = fits.open(path)[0].data[x0:x0+xx,y0:y0+yy]
data_unmod = data.copy()
# sex = fits.open('check.fits')[0].data[500:500+xx,1200:1200+yy]
data -= data.min() - 0.001
# data = 1.-plt.imread('data/sdss.png').T[0]
# data = fits.open(path)[1].data
data = data.clip(min=0.0001)
sex = fits.open('check.fits')[0].data
hdu = fits.PrimaryHDU(data)
hdul = fits.HDUList([hdu])
hdul.writeto('my_data.fits',overwrite=True)
data = np.ndarray.astype(data, float)
vmin = np.log(data.min()+0.2)
vmax = np.log(data.max())*0.4
plt.gray()
lin_max = data.max()*0.01
lin_min=0.01
plt.imsave('log_data.png', np.log(data),vmin=vmin,vmax=vmax)
plt.imsave('data.png', (data), vmax = lin_max,vmin =lin_min)
plt.imsave('sex.png', sex, vmax = lin_max, vmin=lin_min)
plt.imsave('log_sex.png', np.log(sex), vmax = vmax, vmin=vmin)
plt.imsave('point_sex.png', (data_unmod - sex), vmax = lin_max,vmin =lin_min)
plt.imsave('log_point_sex.png', np.log((data_unmod - sex).clip(min=0.0001)), vmax = vmax,vmin =vmin)
alpha = 1.4
Starblade = sb.build_starblade(data, alpha=alpha, cg_steps=100, newton_steps=3)
for i in range(10):
Starblade = sb.starblade_iteration(Starblade, samples=5)
#plotting on logarithmic scale
plt.imsave('log_diffuse_component.png', Starblade.s.val, vmin=vmin, vmax=vmax)
plt.imsave('diffuse_component.png', np.exp(Starblade.s.val), vmin=lin_min, vmax=lin_max)
plt.imsave('log_pointlike_component.png', Starblade.u.val, vmin=vmin, vmax=vmax)
plt.imsave('pointlike_component.png', np.exp(Starblade.u.val), vmin=lin_min, vmax=lin_max)
plt.figure()
k_lenghts = Starblade.power_spectrum.domain[0].k_lengths
plt.plot(k_lenghts, Starblade.power_spectrum.val)
plt.title('power spectrum')
plt.yscale('log')
plt.xscale('log')
plt.ylabel('power')
plt.xlabel('harmonic mode')
plt.savefig('power_spectrum.png')
plt.close('all')
hdu = fits.PrimaryHDU(np.exp(Starblade.u.val))
hdul = fits.HDUList([hdu])
hdul.writeto('my_points.fits',overwrite=True)
hdu = fits.PrimaryHDU(np.exp(Starblade.s.val))
hdul = fits.HDUList([hdu])
hdul.writeto('my_diffuse.fits',overwrite=True)
star_sex = np.log((data_unmod - sex).clip(min=0.0000000001))
star_blade = Starblade.u.val
diffuse_sex = np.log((sex).clip(min=0.0000000001))
diffuse_blade = Starblade.s.val
stars = np.empty((2,len(star_blade.flatten())))
# stars = np.concatenate((star_blade.flatten(),star_sex.flatten()))
stars[0] = star_blade.flatten()
stars[1] = star_sex.flatten()
# stars = sp.whiten(stars)
sp.kmeans2(stars.T,2,iter=30)
#:kivy 1.0.9
TextInput:
on_parent:self.focus = True
<MyWidget>:
image_widget: image_widget
menu_widget: menu_widget
data_path: 'hst_05195_01_wfpc2_f702w_pc_sci.fits'
result_path: ''
alpha: 1.5
spacing: 10
orientation: 'horizontal'
MenuWidget:
id: menu_widget
size_hint: 0.5,1
ImageWidget:
id: image_widget
size_hint: 1,1
<MyAlphaWidget>:
orientation: 'horizontal'
Label:
text: unichr(945)
FloatInput
text: '1.5'
on_text: app.set_alpha(self.text)
<IterationWidget>:
orientation: 'horizontal'
Label:
text: 'iterations'
IntInput
text: '3'
on_text: app.set_iterations(self.text)
<FloatInput>:
multiline: False
<IntInput>:
multiline: False
<MyImage>:
text: self.text
source: self.source
img: img
orientation: 'vertical'
spacing: 10
Label:
text: self.parent.text
size_hint: 1,0.1
Image:
id: img
source: self.parent.source
<AllImageWidget>:
data: data
points: points
diffuse: diffuse
power:power
orientation: 'vertical'
spacing: 10
BoxLayout
spacing: 10
orientation: 'horizontal'
MyImage:
id: data
text: 'data'
source: app.data_image
MyImage:
id: diffuse
source: app.diffuse_image
text: 'diffuse'
BoxLayout:
spacing: 10
orientation: 'horizontal'
MyImage:
id:points
source: app.points_image
text: 'points'
MyImage:
id: power
source: app.power_image
text: 'power'
<DisplayWidget>
all: all
points: points
diffuse: diffuse
data: data
power: power
AllImageScreen:
id:all
name: 'all'
on_pre_enter: self.reload()
SingleImageScreen:
id:data
name: 'data'
text: 'data'
source: app.data_image
on_pre_enter: self.reload()
SingleImageScreen:
id: diffuse
name: 'diffuse'
text: 'diffuse'
source: app.diffuse_image
on_pre_enter: self.reload()
SingleImageScreen:
id: points
name: 'points'
text: 'points'
source: app.points_image
on_pre_enter: self.reload()
SingleImageScreen:
id:power
name: 'power'
text: 'power'
source: app.power_image
on_pre_enter: self.reload()
<SingleImageScreen>
img: img
name: ''
source: self.source
text: ''
MyImage
id: img
text: self.parent.text
source: self.parent.source
<AllImageScreen>
name: ''
all : all
AllImageWidget
id:all
<MenuWidget>:
orientation: 'vertical'
MyAlphaWidget:
size_hint: 1,0.1
IterationWidget:
size_hint: 1,0.1
ActionWidget
<ActionWidget>:
orientation: 'horizontal'
Button:
text: 'load data'
on_press: app.root.current = 'file'
Button:
text: 'run separation'
on_press: app.run_separation()
Button:
text: 'save results'
on_press: app.save_results()
<DisplayChoiceWidget>:
orientation: 'horizontal'
Button:
text: 'all'
on_press: app.image_widget.current = 'all'
Button:
text: 'data'
on_press: app.image_widget.current = 'data'
Button:
text: 'diffuse'
on_press: app.image_widget.current = 'diffuse'
Button:
text: 'points'
on_press: app.image_widget.current = 'points'
Button:
text: 'power'
on_press: app.image_widget.current = 'power'
<ImageWidget>:
orientation: 'vertical'
image_widget : image_widget
DisplayChoiceWidget:
size_hint: 1, 0.1
DisplayWidget:
id: image_widget
<MyPathBrowser>:
id: _filebrowser
dirselect: True
<MyFileBrowser>:
<GlobalScreenManager>:
file: file
main: main
path: path
MainScreen:
id: main
name: 'main'
FileScreen:
id:file
name: 'file'
PathScreen:
id: path
name: 'path'
<FileScreen>:
MyFileBrowser:
on_success: app.load_data(self.selection)
on_canceled: app.root.current = 'main'
<MainScreen>:
image_widget: image_widget
MyWidget:
id:image_widget
<PathScreen>:
MyPathBrowser
filters: [root.is_dir]
on_success: app.select_path(self.selection)
on_canceled: app.root.current = 'main'
import matplotlib
matplotlib.use('agg')
# matplotlib.use('module://kivy.garden.matplotlib.backend_kivy')
from astropy.io import fits
from kivy.app import App
from kivy.uix.widget import Widget
from kivy.uix.button import Button
from kivy.uix.label import Label
from kivy.uix.image import Image
from kivy.properties import ObjectProperty, StringProperty, NumericProperty
from kivy.uix.boxlayout import BoxLayout
from kivy.clock import Clock, mainthread
from os.path import sep, expanduser, isdir, dirname, join
from kivy.garden.filebrowser import FileBrowser
from kivy.utils import platform
from kivy.uix.textinput import TextInput
from kivy.uix.screenmanager import ScreenManager, Screen, NoTransition
import numpy as np
import re
import threading
import matplotlib.pyplot as plt
from kivy.uix.progressbar import ProgressBar
import starblade as sb
def load_data(path):
if path[-5:] == '.fits':
data = fits.open(path)[1].data
else:
data = plt.imread(path)
data = data.clip(min=0.001)
data = np.ndarray.astype(data, float)
return data
class FloatInput(TextInput):
pat = re.compile('[^0-9]')
def insert_text(self, substring, from_undo=False):
print substring
pat = self.pat
if '.' in self.text:
s = re.sub(pat, '', substring)
else:
s = '.'.join([re.sub(pat, '', s) for s in substring.split('.', 1)])
return super(FloatInput, self).insert_text(s, from_undo=from_undo)
class IntInput(TextInput):
pat = re.compile('[^0-9]')
def insert_text(self, substring, from_undo=False):
pat = self.pat
s = re.sub(pat, '', substring)
return super(IntInput, self).insert_text(s, from_undo=from_undo)
class MyImage(BoxLayout):
text = StringProperty('')
source = StringProperty('')
def reload(self):
self.img.reload()
class MyAlphaWidget(BoxLayout):
alpha = NumericProperty(None)
pass
class IterationWidget(BoxLayout):
iteration = NumericProperty(None)
pass
class ResultsPathWidget(BoxLayout):
pass
class DataPathWidget(BoxLayout):
pass
class DisplayWidget(ScreenManager):
def reload(self):
for child in self.children:
child.reload()
class ImageWidget(BoxLayout):
def reload(self):
self.image_widget.reload()
class MenuWidget(BoxLayout):
pass
class SingleImageScreen(Screen):
source = StringProperty('')
def reload(self):
self.img.reload()
pass
class AllImageScreen(Screen):
def reload(self):
self.all.reload()
pass
class AllImageWidget(BoxLayout):
def reload(self):
self.data.reload()
self.points.reload()
self.diffuse.reload()
self.power.reload()
pass
class ActionWidget(BoxLayout):
pass
class DisplayChoiceWidget(BoxLayout):
pass
class GlobalScreenManager(ScreenManager):
pass
class MainScreen(Screen):
pass
class FileScreen(Screen):
pass
class PathScreen(Screen):
def is_dir(self, directory, filename):
return isdir(join(directory, filename))
class MyWidget(BoxLayout):
image_widget = ObjectProperty(None)
menu_widget = ObjectProperty(None)
data_path = StringProperty(None)
result_path = StringProperty(None)
alpha = NumericProperty(None)
class MyPathBrowser(FileBrowser):
pass
class MyFileBrowser(FileBrowser):
filters = ['*.fits', '*.png', '*.jpg']
pass
class StarbladeApp(App):
stop = threading.Event()
data_path = StringProperty(None)
result_path = StringProperty(None)
alpha = NumericProperty(None)
data_image = StringProperty(None)
diffuse_image = StringProperty(None)
points_image = StringProperty(None)
power_image = StringProperty(None)
vmin = None
vmax = None
myEnergy = None
iterations = 3
user_path = ''
reconstructing = False
data_loaded = False
def build(self):