Updated Siemens amplitude, still needs to be tested

proxtoolbox/Algorithms/RAAR.py

@@ -117,7 +117,7 @@ class RAAR(Algorithm):
               if hasattr(self, 'truth'):
                   if self.truth_dim[0] == 1:
                       z=tmp3[0,:]
-                  elif slef.truth_dim[1] == 1:
+                  elif self.truth_dim[1] == 1:
                       z=tmp3[:,0]
                   else:
                       z=tmp3;

proxtoolbox/Problems/Phase/Siemens_amplitude_in.py

@@ -61,7 +61,7 @@ new_config = {
     ## able to control (without too much damage)
 
     ## Algorithm:
-    'algortihm' : 'RAAR',#'RAAR', 'AP',  
+    'algorithm' : 'RAAR',#'RAAR', 'AP',  
     'numruns':1, # the only time this parameter will
     # be different than 1 is when we are
     # benchmarking...not something a normal user

proxtoolbox/Problems/Phase/Siemens_processor.py

@@ -5,6 +5,10 @@ from numpy.linalg import norm
 from numpy.fft import fft2, ifft2
 from numpy.random import rand
 import proxtoolbox.Utilities as Utilities
+from pathlib import Path
+import os
+import proxtoolbox.Utilities.GetData as GetData
+import urllib.request
 
 def Siemens_processor(config):
 
@@ -13,6 +17,20 @@ def Siemens_processor(config):
     noise = config['noise']
     snr = config['data_ball']
 
+    my_file = Path("../InputData/phase_retrieval/Siemens_star_200px.mat")
+    if not(my_file.is_file()):
+       print('******************************************************************')
+       print('* Input data missing.  Please download the phase input data from *') 
+       print('*    http://num.math.uni-goettingen.de/data/Phase.tar.gz         *')
+       print('* Save and unpack the Phase.tar.gz datafile in the               *')
+       print('*    ProxMatlab/InputData subdirectory                           *')
+       print('******************************************************************')
+       '''
+            Automatic download
+            if GetData.query_yes_no("Phase retrieval input data could not be found. Do you want to download the Phase retrieval input data?"):
+            os.makedirs('../InputData/phase_retrieval')
+            urllib.request.urlretrieve("http://num.math.uni-goettingen.de/~r.luke/proxtoolbox-matlab.tar.gz","../InputData/phase_retrieval/proxtoolbox-matlab.tar.gz", reporthook=GetData.dlProgress)'''
+
     print('Loading data file: Siemens_star_200px.mat')
     S = np.loadtxt('../InputData/phase_retrieval/Siemens_star_200px.mat')
 
@@ -36,12 +54,12 @@ def Siemens_processor(config):
         # Stmp=zeros(size(S))
         # Stmp((m/2-m/4):(m/2+m/4),(n/2-n/4):(n/2+n/4))=S((m/2-m/4):(m/2+m/4),(n/2-n/4):(n/2+n/4))
         # S=Stmp
-        config['supp_ampl'] = npnonzero(S)
+        config['support_idx'] = np.nonzero(S)
         
         # use the abs_illumination field to represent the
         # support constraint.
-        config['abs_illumination'] = config['norm_rt_data']*S/(norm(S,'fro'))
-        config['abs_illumination'] = config['abs_illumination']/norm(config['abs_illumination'],'fro')*config['norm_rt_data'][0]
+        config['amplitude'] = config['norm_rt_data']*S/(norm(S,'fro'))
+        config['amplitude'] = config['amplitude']/norm(config['amplitude'],'fro')*config['norm_rt_data'][0]
         config['supp_phase'] = find(ones(m,n))
         config['illumination_phase'] = find(ones(m,n))
     elif config['object'] == 'complex':
@@ -49,10 +67,10 @@ def Siemens_processor(config):
         points = config['Nx']
         config['norm_rt_data']=norm(S,'fro')
         config['norm_data']=config['norm_rt_data']**2
-        # use the abs_illumination field to represent the
+        # use the amplitude field to represent the
         # support constraint.
-        config['abs_illumination'] = config['norm_rt_data']*S/(norm(S,'fro'))
-        # input.abs_illumination = input.abs_illumination/norm(input.abs_illumination,'fro')*input.norm_rt_data(1)
+        config['amplitude'] = config['norm_rt_data']*S/(norm(S,'fro'))
+        # input.amplitude = input.amplitude/norm(input.amplitude,'fro')*input.norm_rt_data(1)
         G=np.zeros(S.shape)# Gaussian(points,10,[points/2+1,points/2+1])
         W=S*np.exp(1j*2*np.pi*G)
         M=abs(fft2(W))
@@ -63,7 +81,7 @@ def Siemens_processor(config):
         # in Luke.m this is changed to the magnitude.
         config['data']=M**M
         config['data_zeros'] = np.where(M==0)
-        config['supp_ampl'] = np.nonzero(S)    
+        config['support_idx'] = np.nonzero(S)    
     elif config['object'] == 'phase':
         # put some phase across S
         points = config['Nx']
@@ -82,7 +100,7 @@ def Siemens_processor(config):
         # in Luke.m this is changed to the magnitude.
         config['data']=M**M
         config['data_zeros'] = np.where(M==0)
-        config['supp_ampl'] = np.nonzero(W)    
+        config['support_idx'] = np.nonzero(W)    
         config['supp_phase'] = S
         config['illumination_phase'] = S    
     
@@ -93,6 +111,8 @@ def Siemens_processor(config):
 
     config['product_space_dimension'] = 1
     config['truth'] = S
+    config['truth_dim'] = np.shape(S)
+    config['norm_truth']=norm(S,'fro')
 
     
 

proxtoolbox/ProxOperators/proxoperators.py

@@ -10,7 +10,7 @@ The "ProxOperators"-module contains various specific operators that do the actua
 import numpy as np
 from numpy import conj, dot, empty, ones, sqrt, sum, zeros, exp, nonzero, log, tile, shape, real
 #from pyfftw.interfaces.scipy_fftpack import fft2, ifft2
-from numpy.fft import fft2, ifft2
+from numpy.fft import fft2, ifft2, ifft
 
 __all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp", "P_SP", "Approx_PM_Gaussian", "Approx_PM_Poisson", "P_S"]
 

proxtoolbox/Utilities/ZeroPad.py

@@ -54,16 +54,16 @@ def ZeroPad(m):
       tmp = np.zeros((rightpad,minor))
       padded_m = np.concatenate((padded_m,tmp),axis=0)
 
-    major = 2**n
+    major = int(2**n)
 
     if minor!=1:
       tmp = np.log2(minor)
       tmp = np.round(2**n-2**tmp)
       if np.mod(tmp,2)==1:
-        leftpad =  (tmp+1)/2
-        rightpad = (tmp-1)/2
+        leftpad =  int((tmp+1)/2)
+        rightpad = int((tmp-1)/2)
       else:
-        leftpad =  tmp/2
+        leftpad =  int(tmp/2)
         rightpad = leftpad
       if orient==2:
         tmp = np.zeros((leftpad,major))