Merge branch 'master' of gitlab.gwdg.de:nam/ProxPython

demos/Elser_RRR.py

+import pdb
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.Elser_Experiment import Elser_Experiment
+
+Elser = Elser_Experiment(Atoms=400, category='H', algorithm='RRR', lambda_0=0.95, lambda_max=0.95, anim=True)
+Elser.run()
+Elser.show()

docs/Tutorial/addingExperiment.rst

@@ -69,7 +69,7 @@ A concrete experiment class needs to override the
       GetData.getData('Phase')
 
       print('Loading data file CDI_intensity')
-      f = loadmat('../InputData/Phase/CDI_intensity.mat')
+      f = loadmat(datadir/'Phase'/'CDI_intensity.mat')
       # diffraction pattern
       dp = f['intensity']['img'][0, 0]
       orig_res = max(dp.shape[0], dp.shape[1])  # actual data size

proxtoolbox/experiments/CT/ART_Experiment.py

@@ -6,6 +6,7 @@ from proxtoolbox.utils.graphics import addColorbar
 
 #for downloading data
 import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
 
 import numpy as np
 from scipy.io import loadmat
@@ -74,7 +75,7 @@ class ART_Experiment(Experiment):
         # load data
         if not self.silent:
             print('Loading data file ART_SheppLogan.mat ')
-        data_Shepp = loadmat('../InputData/CT/ART_SheppLogan.mat')
+        data_Shepp = loadmat(datadir/'CT'/'ART_SheppLogan.mat')
         N = data_Shepp['N'].item()
         p = data_Shepp['p'].item()
         theta = data_Shepp['theta']

proxtoolbox/experiments/orbitaltomography/degenerate_orbits.py

@@ -5,14 +5,16 @@ from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMole
 from proxtoolbox.utils.orbitaltomog import shifted_fft, fourier_interpolate
 from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
 
 class DegenerateOrbital(PlanarMolecule):
     @staticmethod
     def getDefaultParameters():
         defaultParams = {
             'experiment_name': '2D ARPES',
-            'data_filename': '..\\Inputdata\\OrbitalTomog'
-                             + '\\2020_10_27_coronene_Homo1+2_ARPES_2e6counts_corrected_80x80.tif',
+            'data_filename': datadir/'OrbitalTomog'/'2020_10_27_coronene_Homo1+2_ARPES_2e6counts_corrected_80x80.tif',
             'from_intensity_data': True,
             'object': 'real',
             'degeneracy': 2,  # Number of degenerate states to reconstruct
@@ -59,6 +61,10 @@ class DegenerateOrbital(PlanarMolecule):
             - threshold_for_support: float, in range [0,1], fraction of the maximum at which to threshold when
                 determining support or support for sparsity
         """
+
+        # make sure input data can be found, otherwise download it
+        GetData.getData('OrbitalTomog')
+        
         super(DegenerateOrbital, self).loadData()
 
     def createRandomGuess(self):

proxtoolbox/experiments/orbitaltomography/molecule_3d.py

@@ -7,6 +7,9 @@ from proxtoolbox.experiments.orbitaltomography.orbitExperiment import OrbitalTom
 from proxtoolbox.utils.visualization.stack_viewer import XYZStackViewer
 from proxtoolbox.utils.orbitaltomog import shifted_fft
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
 
 class Molecule3D(OrbitalTomographyExperiment):
     @staticmethod
@@ -56,6 +59,9 @@ class Molecule3D(OrbitalTomographyExperiment):
         """
         Load data and set in the correct format for reconstruction
         """
+        # make sure input data can be found, otherwise download it
+        GetData.getData('OrbitalTomog')
+        
         # TODO: copy from data processor
         raise NotImplementedError
 

proxtoolbox/experiments/orbitaltomography/orbitExperiment.py

@@ -6,6 +6,9 @@ from glob import glob
 from proxtoolbox.utils.h5 import write_dict_to_h5
 import numpy as np
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
 
 class OrbitalTomographyExperiment(Experiment):
     @staticmethod
@@ -30,6 +33,10 @@ class OrbitalTomographyExperiment(Experiment):
         self.output_ignore_keys = ['u1', 'u2', 'plots']
 
     def loadData(self):
+        
+        # make sure input data can be found, otherwise download it
+        GetData.getData('OrbitalTomog')
+
         """
         Load or generate the dataset that will be used for
         this experiment. Create the initial iterate.

proxtoolbox/experiments/orbitaltomography/orbit_kpeem.py

@@ -7,6 +7,10 @@ from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMole
 from proxtoolbox.utils.orbitaltomog import shifted_fft, fourier_interpolate, bin_array, shifted_ifft
 from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
+
 
 class OrbitalMomentumMicroscope(PlanarMolecule):
     """
@@ -17,8 +21,7 @@ class OrbitalMomentumMicroscope(PlanarMolecule):
     def getDefaultParameters():
         defaultParams = {
             'experiment_name': 'Momentum microscopy',
-            'data_filename': '..\\Inputdata\\OrbitalTomog'
-                             + '\\2020_11_05_Coronene_HOMO_near_degenerate_noisefree.tif',
+            'data_filename': datadir/'OrbitalTomog'/'2020_11_05_Coronene_HOMO_near_degenerate_noisefree.tif',
             'from_intensity_data': False,
             'object': 'real',
             'constraint': 'sparse real',
@@ -54,6 +57,10 @@ class OrbitalMomentumMicroscope(PlanarMolecule):
         self.momentum_axes = (-1, -2)
 
     def loadData(self):
+
+        # make sure input data can be found, otherwise download it
+        GetData.getData('OrbitalTomog')
+        
         """
         Load data and set in the correct format for reconstruction
         Parameters are taken from experiment class (self) properties, which must include::

proxtoolbox/experiments/orbitaltomography/orthogonal_orbits.py

@@ -7,13 +7,18 @@ from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMole
 from proxtoolbox.utils.orbitaltomog import shifted_fft, fourier_interpolate, bin_array, shifted_ifft
 from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
+
 
 class OrthogonalOrbitals(PlanarMolecule):
     @staticmethod
     def getDefaultParameters():
         defaultParams = {
             'experiment_name': '2D ARPES',
-            'data_filename': '..\\Inputdata\\OrbitalTomog'
+            'data_filename':  datadir
+                             + '\OrbitalTomog'
                              + '\\2020_10_27_coronene_Homo_stack_ARPES_2e6counts_corrected_80x80.tif',
             'from_intensity_data': False,
             'object': 'real',
@@ -52,6 +57,10 @@ class OrthogonalOrbitals(PlanarMolecule):
             - threshold_for_support: float, in range [0,1], fraction of the maximum at which to threshold when
                 determining support or support for sparsity
         """
+        
+        # make sure input data can be found, otherwise download it
+        GetData.getData('OrbitalTomog')
+
         # load data
         if self.data_filename is None:
             self.data_filename = input('Please enter the path to the datafile: ')

proxtoolbox/experiments/orbitaltomography/planar_molecule.py

@@ -7,6 +7,9 @@ from proxtoolbox.experiments.orbitaltomography.orbitExperiment import OrbitalTom
 from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
 from proxtoolbox.utils.orbitaltomog import bin_array, shifted_fft, shifted_ifft, fourier_interpolate, roll_to_pos
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
 
 class PlanarMolecule(OrbitalTomographyExperiment):
     @staticmethod
@@ -14,7 +17,7 @@ class PlanarMolecule(OrbitalTomographyExperiment):
         # TODO: optimize parameters and proxoperators to get good & consistent phase retrieval using the demo
         defaultParams = {
             'experiment_name': 'noisy 2D ARPES', # '2D ARPES', #
-            'data_filename': '../InputData/OrbitalTomog/coronene_homo1.tif',
+            'data_filename': datadir/'OrbitalTomog'/'coronene_homo1.tif',
             'from_intensity_data': False,
             'object': 'real',
             'constraint': 'sparse real',
@@ -74,6 +77,9 @@ class PlanarMolecule(OrbitalTomographyExperiment):
             - threshold_for_support: float, in range [0,1], fraction of the maximum at which to threshold when
                 determining support or support for sparsity
         """
+        # make sure input data can be found, otherwise download it
+        GetData.getData('OrbitalTomog')
+        
         # load data
         if self.data_filename is None:
             self.data_filename = input('Please enter the path to the datafile: ')

proxtoolbox/experiments/phase/CDI_Experiment.py

@@ -7,6 +7,7 @@ from proxtoolbox.utils.graphics import addColorbar
 
 #for downloading data
 import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
 
 import numpy as np
 from numpy import exp, sqrt, log2, log10, ceil, floor, unravel_index, argmax, zeros
@@ -77,7 +78,7 @@ class CDI_Experiment(PhaseExperiment):
         GetData.getData('Phase')
         if not self.silent:
             print('Loading data file CDI_intensity')
-        f = loadmat('../InputData/Phase/CDI_intensity.mat')
+        f = loadmat(datadir/'Phase'/'CDI_intensity.mat')
         # diffraction pattern
         dp = f['intensity']['img'][0, 0]
         orig_res = max(dp.shape[0], dp.shape[1])  # actual data size

proxtoolbox/experiments/phase/CDP_Experiment.py

@@ -13,6 +13,11 @@ from numpy.linalg import norm
 from numpy.fft import fft2, ifft2
 import time
 
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
+print(datadir)
+
 
 class CDP_Experiment(PhaseExperiment):
     """
@@ -69,6 +74,9 @@ class CDP_Experiment(PhaseExperiment):
         """
         Load CDP dataset. Create the initial iterate.
         """
+        #make sure input data can be found, otherwise download it
+        GetData.getData('Phase')
+        
         # Implementation of the Wirtinger Flow (WF) algorithm presented in the paper 
         # "Phase Retrieval via Wirtinger Flow: Theory and algorithms"
         # by E. J. Candes, X. Li, and M. Soltanolkotabi
@@ -91,18 +99,18 @@ class CDP_Experiment(PhaseExperiment):
         # make image
         if debug:
             if n2 == 1:
-                x_dict = loadMatFile('../InputData/Phase/CDP_1D_x.mat')
+                x_dict = loadMatFile(datadir/'Phase'/'CDP_1D_x.mat')
                 debug_image = x_dict['x']
-                masks_dict = loadMatFile('../InputData/Phase/CDP_1D_Masks.mat')
+                masks_dict = loadMatFile(datadir/'Phase'/'CDP_1D_Masks.mat')
                 debug_masks = masks_dict['Masks']
-                z0_dict = loadMatFile('../InputData/Phase/CDP_1D_z0.mat')
+                z0_dict = loadMatFile(datadir/'Phase'/'CDP_1D_z0.mat')
                 debug_z0 = z0_dict['z0']
             elif n2 == 256:
-                x_dict = loadMatFile('../InputData/Phase/CDP_2D_x.mat')
+                x_dict = loadMatFile(datadir/'Phase'/'CDP_2D_x.mat')
                 debug_image = x_dict['x']
-                masks_dict = loadMatFile('../InputData/Phase/CDP_2D_Masks.mat')
+                masks_dict = loadMatFile(datadir/'Phase'/'CDP_2D_Masks.mat')
                 debug_masks = masks_dict['Masks']
-                z0_dict = loadMatFile('../InputData/Phase/CDP_2D_z0.mat')
+                z0_dict = loadMatFile(datadir/'Phase'/'CDP_2D_z0.mat')
                 debug_z0 = z0_dict['z0']
             x = debug_image
         else:

proxtoolbox/experiments/phase/Elser_Experiment.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.phaseExperiment import PhaseExperiment
+from proxtoolbox import proxoperators
+import numpy as np
+from numpy import fromfile, exp, nonzero, zeros, pi, resize, real, angle
+from numpy.random import rand
+from numpy.linalg import norm
+from numpy.fft import fftshift, fft2, ifft2
+
+import proxtoolbox.utils as utils
+from proxtoolbox.utils.cell import Cell, isCell
+
+#for downloading data
+import proxtoolbox.utils.GetData as GetData
+from proxtoolbox.utils.GetData import datadir
+
+import matplotlib
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import subplots, show, figure
+
+
+class Elser_Experiment(PhaseExperiment):
+    '''
+    Elser's Crystallography experiment class
+    '''
+
+    @staticmethod
+    def getDefaultParameters():
+        defaultParams = {
+            'experiment_name' : 'Elser',
+            'object': 'complex',
+            'constraint': 'atom count',
+            'distance': 'far field',
+            'Nx': 128,
+            'Ny': 128,
+            'Nz': 1,
+            'Atoms':100,
+            'category': 'E',
+            'product_space_dimension': 2,
+            'MAXIT': 10000,
+            'TOL': 5e-5,
+            'lambda_0': 0.85,
+            'lambda_max': 0.85,
+            'lambda_switch': 20,
+            'data_ball': 1e-30,
+            'diagnostic': True,
+	    'iterate_monitor_name': 'FeasibilityIterateMonitor',
+            'verbose': 0,
+            'graphics': 1
+        }
+        return defaultParams
+
+
+    def __init__(self, Atoms=100, category='E', **kwargs):
+        super(Elser_Experiment, self).__init__(**kwargs)
+        self.fresnel_nr = None
+        self.use_farfield_formula = None
+        self.data_zeros = None
+        self.supp_ampl = None
+        self.indicator_ampl = None
+        self.abs_illumination = None
+        self.illumination_phase = None
+        self.FT_conv_kernel = Cell(1)
+        self.Atoms=Atoms
+        self.category=category
+    def loadData(self):
+    # def loadElserData(filename, M):
+        """
+        Load Elser dataset. Create the initial iterate.
+        """
+        #make sure input data can be found, otherwise download it
+        GetData.getData('Elser')
+	#defaultParams = getDefaultParameters()
+        # So far only Nx to set the desired
+        # resolution
+        newres = self.Nx
+
+        # The following value for snr does not work well in Python
+        # as this creates numerical issues with Poisson noise
+        # (data_ball is far too small)
+        # snr = self.data_ball
+        # Instead, we choose the following bigger value which
+        # works well.
+        snr = 1e-8
+
+        print('Loading data.')
+        fn="".join([str(datadir/"Elser"/"data"), str(self.Atoms), self.category])
+        # read text data file as a matrix of numbers
+        data = np.loadtxt(fn, delimiter="\t")
+        data = np.sqrt(data)/self.Nx
+        data = np.c_[data, np.zeros(self.Nx)] # add extra column
+        print(data)
+        data = data.flatten() # make it a one-dimensional array
+        M = fftshift(data)
+        # have to look at what Elser is actually doing with his data...
+        self.magn = 1  # magnification factor
+        self.farfield = True
+        print('Using farfield formula.')
+
+        self.rt_data = Cell(1)
+        self.rt_data[0] = M
+        # standard for the main program is that
+        # the data field is the magnitude SQUARED
+        # in Luke.m this is changed to the magnitude.
+        self.data = Cell(1)
+        self.data[0] = M**2
+        self.norm_rt_data = np.linalg.norm(ifft2(M), 'fro')
+        self.data_zeros = np.where(M == 0)
+        self.support_idx = np.nonzero(S)
+        self.sets = 2
+
+        # The support constraint is sparsity determined by the number of atoms indicated in the
+        # data file name...this needs to be installed
+
+        # initial guess: follow Elser's initialization...
+
+    def setupProxOperators(self):
+        """
+        Determine the prox operators to be used for this experiment
+        """
+        super(Elser_Experiment, self).setupProxOperators()  # call parent's method
+
+        self.propagator = 'Propagator_FreFra'
+        self.inverse_propagator = 'InvPropagator_FreFra'
+
+        # remark: self.farfield is always true for these problems
+        # self.proxOperators already contains a prox operator at slot 0.
+        # Here, we add the second one.
+        if self.constraint == 'phaselift':
+            self.proxOperators.append('P_Rank1')
+        elif self.constraint == 'phaselift2':
+            self.proxOperators.append('P_rank1_SR')
+        else:
+            self.proxOperators.append('Approx_Pphase_FreFra_Poisson')
+        self.nProx = self.sets
+
+        # Now we adjust data structures and prox operators according to the algorithm
+        # What follows wa set up for JWST and works in that context.  Not
+        # sure how much of this is trasferrable to Elsers problems (DRL 14.08.2020)
+        if self.algorithm_name == 'ADMM':
+            # set-up variables in cells for block-wise algorithms
+            # ADMM has three blocks of variables, primal, auxilliary (primal
+            # variables in the image space of some linear mapping) and dual.
+            # we sort these into a 3D cell.
+            u0 = self.u0
+            self.u0 = Cell(3)
+            self.proxOperators = []
+            self.productProxOperators = []
+            K = self.product_space_dimension - 1
+            self.product_space_dimension = K
+            self.u0[0] = u0[K]
+            prox_FreFra_Poisson_class = getattr(proxoperators, 'Approx_Pphase_FreFra_Poisson')
+            prox_FreFra_Poisson = prox_FreFra_Poisson_class(self)
+            self.u0[1] = Cell(K)
+            self.u0[2] = Cell(K)
+            for j in range(K):
+                tmp_u0 = prox_FreFra_Poisson.eval(u0[j], j)
+                self.u0[1][j] = fft2(self.FT_conv_kernel[j]*tmp_u0) / (self.Nx*self.Ny)
+                self.u0[2][j] = self.u0[1][j] / self.lambda_0
+            self.proxOperators.append('Prox_primal_ADMM_indicator')
+            self.proxOperators.append('Approx_Pphase_FreFra_ADMM_Poisson')
+        elif self.algorithm_name == 'AvP2':
+            # u_0 should be a cell...we change it into a cell of cells
+            u0 = self.u0
+            self.u0 = Cell(3)
+            prox2_class = getattr(proxoperators, self.proxOperators[1])
+            prox2 = prox2_class(self)
+            self.u0[2] = prox2.eval(u0)
+            P_Diag_class = getattr(proxoperators, 'P_diag')
+            P_Diag_prox = P_Diag_class(self)
+            tmp_u = P_Diag_prox.eval(self.u0[2])
+            self.u0[1] = tmp_u[0]
+            self.u0[0] = u0[self.product_space_dimension-1]
+        elif self.algorithm_name == 'PHeBIE':
+            # set up variables in cells for block-wise, implicit-explicit algorithms
+            u0 = self.u0
+            self.u0 = Cell(2)
+            self.product_space_dimension = self.product_space_dimension - 1
+            self.u0[0] = u0[self.product_space_dimension].copy()  # copy last element of prev u0 into first slot
+            self.u0[1] = Cell(self.product_space_dimension)
+            for j in range(self.product_space_dimension):
+                self.u0[1][j] = u0[j]
+            self.proxOperators = []
+            self.nProx = 2
+            self.proxOperators.append('P_amp_support')  # project onto the support/amplitude constraints
+            self.proxOperators.append('Prox_product_space')
+            self.productProxOperators = []
+            self.n_product_Prox = self.product_space_dimension
+            for _j in range(self.n_product_Prox):
+                self.productProxOperators.append('Approx_Pphase_FreFra_Poisson')
+        elif self.algorithm_name == 'Wirtinger':
+            # Contrary to the Matlab code, we use cells instead of 3D arrays
+            L = len(self.FT_conv_kernel)
+            self.product_space_dimension = L
+            self.masks = Cell(L)
+            for j in range(L):
+                self.masks[j] = self.indicator_ampl * self.FT_conv_kernel[j]
+            self.data_zeros = None
+            self.FT_conv_kernel = self.masks
+            self.proxOperators[1] = 'Approx_Pphase_JWST_Wirt'
+            self.use_farfield_formula = True
+            self.fresnel_nr = 0
+
+    def show(self):
+        """
+        Generate graphical output from the solution
+        """
+
+        # display plot of far field data
+        # figure(123)
+        f, (ax1) = subplots(1, 1,
+                                 figsize=(self.figure_width, self.figure_height),
+                                 dpi=self.figure_dpi)
+        im = ax1.imshow(log10(self.dp + 1e-15))
+        f.colorbar(im, ax=ax1, fraction=0.046, pad=0.04)
+        ax1.set_title('Far field data')
+        plt.subplots_adjust(wspace=0.3)  # adjust horizontal space (width)
+        # between subplots (default = 0.2)
+        f.suptitle('Elser Data')
+
+        # call parent to display the other plots
+        super(Elser_Experiment, self).show()
+
+
+        u_m = self.output['u_monitor']
+        if isCell(u_m):
+            u = u_m[0]
+            if isCell(u):
+                u = u[0]
+            u2 = u_m[len(u_m)-1]
+            if isCell(u2):
+                u2 = u2[len(u2)-1]
+        else:
+            u2 = u_m
+            if u2.ndim > 2:
+                u2 = u2[:,:,0]
+            u = self.output['u']
+            if isCell(u):
+                u = u[0]
+            elif u.ndim > 2:
+                u = u[:,:,0]
+
+        algo_desc = self.algorithm.getDescription()
+        title = "Algorithm " + algo_desc
+
+        # figure 904
+        titles = ["Best approximation amplitude - physical constraint satisfied",
+                  "Best approximation phase - physical constraint satisfied",
+                  "Best approximation amplitude - Fourier constraint satisfied",
+                  "Best approximation phase - Fourier constraint satisfied"]
+        f = self.createFourImageFigure(u, u2, titles)
+        f.suptitle(title)
+        plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height) between subplots (default = 0.2)
+
+        # figure 900
+        changes = self.output['stats']['changes']
+        time = self.output['stats']['time']
+        time_str = "{:.{}f}".format(time, 5) # 5 is precision
+        label = "Iterations (time = " + time_str + " s)"
+
+        f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, \
+                    figsize = (self.figure_width, self.figure_height),
+                    dpi = self.figure_dpi)
+        self.createImageSubFigure(f, ax1, abs(u), titles[0])
+        self.createImageSubFigure(f, ax2, real(u), titles[1])
+
+        ax3.semilogy(changes)
+        ax3.set_xlabel(label)
+        ax3.set_ylabel('Log of change in iterates')
+
+        if 'gaps' in self.output['stats']:
+            gaps = self.output['stats']['gaps']
+            ax4.semilogy(gaps)
+            ax4.set_xlabel(label)
+            ax4.set_ylabel('Log of the gap distance')
+
+        f.suptitle(title)
+        plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height) between subplots (default = 0.2)
+        plt.subplots_adjust(wspace = 0.3) # adjust horizontal space (width) between subplots (default = 0.2)
+
+        show()
+
+if __name__ == "__main__":
+	Elser = Elser_Experiment(Atoms=200, category='E', algorithm='RRR', lambda_0=0.95, lambda_max=0.95, anim=True)
+	Elser.run()
+	Elser.show()

proxtoolbox/experiments/phase/Elser_processor.m

+load('../../../InputData/Phase/Elser/data100E')
+x=data100E;
+datapow=0;
+M=max(size(x));
+fmag=zeros(M, M/2);
+Fmag=zeros(M,M);
+for(i=1:M/2)
+	j=1;
+	datapow=x(i,j)+datapow;
+	fmag(i,j)=sqrt(x(i,j))/M;
+    	Fmag(i,j) = fmag(i,j);
+	for(j=2:M/2)
+		datapow=2*x(i,j)+datapow;
+		fmag(i,j)=sqrt(x(i,j))/M;
+        	Fmag(i,j) = fmag(i,j);
+        	if(i>1)
+            		Fmag(M-(j-2),M-(i-2)) = fmag(i,j);
+        	end
+    	end
+end
+for(i=M/2+1:M)
+    	j=1;
+	datapow=x(i,j)+datapow;
+	fmag(i,j)=sqrt(x(i,j))/M;
+    	Fmag(i,j) = fmag(i,j);
+	for(j=2:M/2)
+		datapow=2*x(i,j)+datapow;
+		fmag(i,j)=sqrt(x(i,j))/M;
+        	Fmag(i,j) = fmag(i,j);
+        	Fmag(j,i) = fmag(i,j);
+    	end
+end