diff --git a/demos/demo.py b/demos/demo.py index 9136c16a55a02abed2b917c0d65516a755e180b9..f21fc5784fead894d5ebf5a30a8446011a97bb62 100644 --- a/demos/demo.py +++ b/demos/demo.py @@ -24,7 +24,7 @@ import starblade as sb if __name__ == '__main__': #specifying location of the input file: - path = 'hst_05195_01_wfpc2_f702w_pc_sci.fits' + path = 'data/hst_05195_01_wfpc2_f702w_pc_sci.fits' data = fits.open(path)[1].data data = data.clip(min=0.001) diff --git a/demos/multichannel_demo.py b/demos/multichannel_demo.py index dbda3a3f4e6bc0a626ce50fbc84a42ca6608a396..79a1f3f41605980058403f1c31527659c4062640 100644 --- a/demos/multichannel_demo.py +++ b/demos/multichannel_demo.py @@ -19,19 +19,22 @@ import numpy as np from matplotlib import pyplot as plt import starblade as sb +import nifty4 as ift if __name__ == '__main__': # data = plt.imread('10Keso1242a.tif') - data = plt.imread('data/eso1242a.jpg') + data = plt.imread('data/galaxy.jpg') data = data.astype(float) data = data.clip(0.0001) - alpha = 1.25 - MultiStarblade = sb.build_multi_starblade(data, alpha) + alpha = 1.3 + MultiStarblade = sb.build_multi_starblade(data, alpha, newton_iterations=1) + # power = MultiStarblade[0].power_spectrum / np.arange(len(MultiStarblade[0].power_spectrum.val))**0 + # MultiStarblade = sb.build_multi_starblade(data,alpha,newton_iterations=1, manual_power_spectrum=power) - for i in range(10): - MultiStarblade = sb.multi_starblade_iteration(MultiStarblade, multiprocessing=True) + 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) @@ -40,5 +43,13 @@ if __name__ == '__main__': 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.) + plt.figure() + + plt.plot(MultiStarblade[0].power_spectrum.val) + plt.yscale('log') + plt.xscale('log') + plt.savefig('power.jpg') + plt.close() diff --git a/starblade/starblade_energy.py b/starblade/starblade_energy.py index 85210d4665478d8b5625f57d788169acd2356330..da760702a88f25313a83bf5af99ada2f119b7310 100644 --- a/starblade/starblade_energy.py +++ b/starblade/starblade_energy.py @@ -16,7 +16,7 @@ # # Starblade is being developed at the Max-Planck-Institut fuer Astrophysik -from nifty4 import Energy, Field, log, exp, DiagonalOperator +from nifty4 import Energy, Field, log, exp, DiagonalOperator, create_power_operator from nifty4.library import WienerFilterCurvature from nifty4.library.nonlinearities import PositiveTanh @@ -38,9 +38,9 @@ class StarbladeEnergy(Energy): Slope parameter of the point-source prior q : Field Cutoff parameter of the point-source prior - correlation : Field - A field in the Fourier space which encodes the diagonal of the prior - correlation structure of the diffuse component + power_spectrum : callable or Field + An object that contains the power spectrum of the diffuse component + as a function of the harmonic mode. FFT : FFTOperator An operator performing the Fourier transform inverter : ConjugateGradient @@ -57,9 +57,12 @@ class StarbladeEnergy(Energy): self.inverter = parameters['inverter'] self.d = parameters['data'] self.FFT = parameters['FFT'] - self.correlation = parameters['correlation'] + self.power_spectrum = parameters['power_spectrum'] + self.correlation = create_power_operator(self.FFT.domain, self.power_spectrum) self.alpha = parameters['alpha'] self.q = parameters['q'] + self.update_power = parameters['update_power'] + self.newton_iterations = parameters['newton_iterations'] pos_tanh = PositiveTanh() self.S = self.FFT * self.correlation * self.FFT.adjoint self.a = pos_tanh(self.position) diff --git a/starblade/sugar.py b/starblade/sugar.py index 403e70455332a1024513b21eeb6968bc11322171..7a235c8a5f5362f14ece9157af9c417b03eff39f 100644 --- a/starblade/sugar.py +++ b/starblade/sugar.py @@ -23,7 +23,8 @@ import nifty4 as ift from .starblade_energy import StarbladeEnergy -def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500): +def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500, newton_iterations = 3, + manual_power_spectrum = None): """ Setting up the StarbladeEnergy for the given data and parameters Parameters ---------- @@ -35,6 +36,11 @@ def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500): The cutoff parameter of the point source prior (default: 1e-40). cg_iterations : int Maximum number of conjugate gradient iterations for numerical operator inversion (default: 500). + newton_iterations : int + Number of consecutive Newton steps within one algorithmic step.(default: 3) + manual_power_spectrum : None, Field or callable + Option to set a manual power spectrum which is kept constant during the separation. + If it is not specified, the algorithm will try to infer it via critical filtering. """ s_space = ift.RGSpace(data.shape, distances=len(data.shape) * [1]) @@ -43,9 +49,14 @@ def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500): FFT = ift.FFTOperator(h_space, target=s_space) binbounds = ift.PowerSpace.useful_binbounds(h_space, logarithmic = False) p_space = ift.PowerSpace(h_space, binbounds=binbounds) - initial_spectrum = ift.power_analyze(FFT.inverse_times(ift.log(data)), binbounds=p_space.binbounds) - initial_spectrum /= (p_space.k_lengths+1.)**2 - initial_correlation = ift.create_power_operator(h_space, initial_spectrum) + if manual_power_spectrum is None: + initial_spectrum = ift.power_analyze(FFT.inverse_times(ift.log(data)), + binbounds=p_space.binbounds) + initial_spectrum /= (p_space.k_lengths+1.)**2 + update_power = True + else: + initial_spectrum = manual_power_spectrum + update_power = False initial_x = ift.Field(s_space, val=-1.) alpha = ift.Field(s_space, val=alpha) q = ift.Field(s_space, val=q) @@ -53,38 +64,40 @@ def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500): tol_abs_gradnorm=1e-5) inverter = ift.ConjugateGradient(controller=ICI) - parameters = dict(data=data, correlation=initial_correlation, + parameters = dict(data=data, power_spectrum=initial_spectrum, alpha=alpha, q=q, - inverter=inverter, FFT=FFT) + inverter=inverter, FFT=FFT, + newton_iterations=newton_iterations, update_power=update_power) Starblade = StarbladeEnergy(position=initial_x, parameters=parameters) return Starblade -def starblade_iteration(starblade, iterations=3): +def starblade_iteration(starblade): """ Performing one Newton minimization step Parameters ---------- starblade : StarbladeEnergy An instance of an Starblade Energy - iterations : int - The number of steps with the Newton scheme (default: 3). """ - controller = ift.GradientNormController(name="Newton", tol_abs_gradnorm=1e-8, iteration_limit=iterations) + controller = ift.GradientNormController(name="Newton", + tol_abs_gradnorm=1e-8, + iteration_limit=starblade.newton_iterations) minimizer = ift.RelaxedNewton(controller=controller) energy, convergence = minimizer(starblade) new_position = energy.position - h_space = energy.correlation.domain[0] - FFT = energy.FFT - binbounds = ift.PowerSpace.useful_binbounds(h_space, logarithmic=False) - new_power = ift.power_analyze(FFT.inverse_times(energy.s), binbounds=binbounds) - # new_power /= (new_power.domain[0].k_lengths+1.)**2 - - new_correlation = ift.create_power_operator(h_space, new_power) new_parameters = energy.parameters - # new_parameters['correlation'] = new_correlation + if energy.update_power: + h_space = energy.correlation.domain[0] + FFT = energy.FFT + binbounds = ift.PowerSpace.useful_binbounds(h_space, logarithmic=False) + new_power = ift.power_analyze(FFT.inverse_times(energy.s), binbounds=binbounds) + # new_power /= (new_power.domain[0].k_lengths+1.)**2 + new_parameters['power_spectrum'] = new_power + NewStarblade = StarbladeEnergy(new_position, new_parameters) return NewStarblade -def build_multi_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500): +def build_multi_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500, + newton_iterations=3, manual_power_spectrum = None): """ Builds a list of StarbladeEnergies for the given multi-channel dataset Parameters ---------- @@ -96,24 +109,33 @@ def build_multi_starblade(data, alpha=1.5, q=1e-40, cg_iterations=500): The cutoff parameter of the point source prior (default: 1e-40). cg_iterations : int Maximum number of conjugate gradient iterations for numerical operator inversion (default: 500). + newton_iterations : int + Number of consecutive Newton steps within one algorithmic step.(default: 3) + manual_power_spectrum : None, Field or callable + Option to set a manual power spectrum which is kept constant during the separation. + If it is not specified, the algorithm will try to infer it via critical filtering. """ MultiStarblade = [] for i in range(data.shape[-1]): - starblade = build_starblade(data[...,i],alpha=alpha, q=q, cg_iterations=cg_iterations) + starblade = build_starblade(data[...,i],alpha=alpha, q=q, + cg_iterations=cg_iterations, + newton_iterations=newton_iterations, + manual_power_spectrum = manual_power_spectrum) MultiStarblade.append(starblade) return MultiStarblade -def multi_starblade_iteration(MultiStarblade, multiprocessing = False): +def multi_starblade_iteration(MultiStarblade, processes = 1): """ Performing one Newton minimization step for all entries of the MultiStarblade list. Parameters ---------- MultiStarblade : list of StarbladeEnergy A list of instances of an Starblade Energy - iterations : int - The number of steps with the Newton scheme (default: 3). + processes : int + Each channel can be computed independently. + This number specifies how many processes are set up.(default: 1) """ - if multiprocessing: - NewStarblades = list(Pool(processes=3).map(starblade_iteration, + if processes>1: + NewStarblades = list(Pool(processes=processes).map(starblade_iteration, MultiStarblade)) else: NewStarblades = []