Commit 5b85ca73 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Moved Utilities to proxtoolbox.Utilities, split up file

Adjusted Input, import paths in JWST_data_processor
parent 8538e813
import sys import sys
sys.path.append('../proxtoolbox/Problems/Phase') sys.path.append('../proxtoolbox/Problems/Phase')
sys.path.append('..')
import JWST_data_processor import JWST_data_processor
import JWST_in import JWST_in
from phase import Phase from phase import Phase
......
# Fourier transform processor. FFT's data in the physical domain and
# corrects for the peculiarities of numpy's fft algorithm.
# First we fft the data, then we take the negative sign of every other
# mode in each direction (x and y), then we
# divide the bad boy by (sqrt(res))^2 to get the correct intensity scaling,
# finally we fftshift everything.
from numpy.random import rand
from numpy import exp, sqrt, log, tan, pi, floor, zeros, ones
from scipy.special import gammaln
import numpy as np
def FFT(f):
shape = f.shape;
res=max(shape);
if shape[0] == shape[1]:
F= np.fft.fft2(f);
F=F/(res);
F[1:res:2,:]=-F[1:res:2,:];
F[:,1:res:2]=-F[:,1:res:2];
else:
F = np.fft.fft(f,axis=0);
F = F.ravel('F'); #create 1-d array, colum-major order (matlab style), not really nice
F[1:res:2] =-F[1:res:2];
F = F.reshape(shape[1],shape[0]).T; #back to original shape
F=F/np.sqrt(res);
return np.fft.fftshift(F);
# Inverse Fourier transform processor. Takes data in the fourier domain
# that has been properly normalized and arranged (i.e. negative intensities
# resulting from Python's fft algorithm flipped to the propper sign) by the
# companion routine FFT.
# First we fftshift the data, then we take the negative sign of every other
# mode in each direction (x and y), and finally we
# multiply the bad boy by (sqrt(res))^2 to get the correct intensity scaling.
def IFFT(F):
shape = F.shape;
res=max(shape);
F= np.fft.fftshift(F);
if shape[0] == shape[1]:
F[1:res:2,:]=-F[1:res:2,:]
F[:,1:res:2]=-F[:,1:res:2];
F=res*F;
f= np.fft.ifft2(F);
else:
F = F.ravel('F'); #create 1-d array, colum-major order (matlab style), not really nice
F[1:res:2] =-F[1:res:2];
F = F.reshape(shape[1],shape[0]).T; #back to original shape
F=np.sqrt(res)*F;
f= np.fft.ifft(F,axis=0);
return f;
# Poisson random number generation with mean of xm
# April 6, 2005 Ning Lei
# see numerical recipe C++ version
# ref: y(i+1) = xm^i*exp(-xm)/factorial(i); see Leo Breiman Probability
def PoissonRan(xm):
oldm = -1;
if xm<12:
if xm != oldm:
oldm = xm;
g = exp(-xm);
em = -1;
t = 1;
em = em+1;
t = t*rand(1);
while t>g :
em = em+1;
t = t*rand(1);
else:
if xm != oldm:
oldm = xm;
sq = sqrt(2.0*xm);
alxm = log(xm);
g = xm*alxm-gammaln(xm+1);
y = tan(pi*rand(1));
em = sq*y+xm;
while em < 0:
y = tan(pi*rand(1));
em = sq*y+xm;
em = floor(em);
t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g);
while rand(1) > t:
y = tan(pi*rand(1));
em = sq*y+xm;
while em < 0:
y = tan(pi*rand(1));
em = sq*y+xm;
em = floor(em);
t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g);
return em;
# resizes arrays. Takes an array A that is M by N and resizes it to Mnew by
# Nnew.
# usage: Anew = Resize(A, Mnew,Nnew)
def Resize(A, Mnew,Nnew):
M = A.shape[0];
N = A.shape[1];
Anew = zeros((Mnew,Nnew));
if Nnew <= N:
Nstep = int(N/Nnew);
Mstep = int(M/Mnew);
countN = -1;
for j in range(0,N,Nstep):
countN = countN+1;
countM = -1;
for i in range(0,M,Mstep):
countM = countM+1;
Anew[countM,countN] = A[i,j];
else:
Nstep = int(Nnew/N);
Mstep = int(Mnew/M);
temp=ones((Mstep,Nstep));
countN = -1;
for j in range(N):
countN = countN+1;
countM = -1;
for i in range(M):
countM = countM+1;
Anew[i*Mstep:(i+1)*Mstep,j*Nstep:(j+1)*Nstep] = A[countM,countN]*temp;
return Anew
...@@ -25,7 +25,7 @@ from numpy import fromfile, exp, nonzero, zeros, pi, resize ...@@ -25,7 +25,7 @@ from numpy import fromfile, exp, nonzero, zeros, pi, resize
from numpy.random import rand from numpy.random import rand
from numpy.linalg import norm from numpy.linalg import norm
from numpy.fft import fftshift from numpy.fft import fftshift
from Utilities import Utilities import proxtoolbox.Utilities as Utilities
def JWST_data_processor(config): def JWST_data_processor(config):
...@@ -36,19 +36,19 @@ def JWST_data_processor(config): ...@@ -36,19 +36,19 @@ def JWST_data_processor(config):
noise = config['noise']; noise = config['noise'];
snr = config['data_ball']; snr = config['data_ball'];
with open('InputData/phase_retrieval/pupil.pmod','r') as fid: with open('../InputData/phase_retrieval/pupil.pmod','r') as fid:
# read lower endian float <f # read lower endian float <f
Xi_A = np.fromfile(fid, dtype='<f'); Xi_A = np.fromfile(fid, dtype='<f');
Xi_A = Xi_A.reshape((512,512)).T; Xi_A = Xi_A.reshape((512,512)).T;
diversity=3; diversity=3;
with open('InputData/phase_retrieval/phase_p37.pmod','r') as fid: with open('../InputData/phase_retrieval/phase_p37.pmod','r') as fid:
# read lower endian float <f # read lower endian float <f
temp1 = np.fromfile(fid, dtype='<f'); temp1 = np.fromfile(fid, dtype='<f');
temp1 = temp1.reshape(512,512).T; temp1 = temp1.reshape(512,512).T;
with open('InputData/phase_retrieval/phase_m37.pmod','r') as fid: with open('../InputData/phase_retrieval/phase_m37.pmod','r') as fid:
# read lower endian float <f # read lower endian float <f
temp2 = np.fromfile(fid, dtype='<f'); temp2 = np.fromfile(fid, dtype='<f');
temp2 = temp2.reshape(512,512).T; temp2 = temp2.reshape(512,512).T;
......
# Fourier transform processor. FFT's data in the physical domain and
# corrects for the peculiarities of numpy's fft algorithm.
# First we fft the data, then we take the negative sign of every other
# mode in each direction (x and y), then we
# divide the bad boy by (sqrt(res))^2 to get the correct intensity scaling,
# finally we fftshift everything.
from numpy.random import rand
from numpy import exp, sqrt, log, tan, pi, floor, zeros, ones
from scipy.special import gammaln
import numpy as np
def FFT(f):
shape = f.shape;
res=max(shape);
if shape[0] == shape[1]:
F= np.fft.fft2(f);
F=F/(res);
F[1:res:2,:]=-F[1:res:2,:];
F[:,1:res:2]=-F[:,1:res:2];
else:
F = np.fft.fft(f,axis=0);
F = F.ravel('F'); #create 1-d array, colum-major order (matlab style), not really nice
F[1:res:2] =-F[1:res:2];
F = F.reshape(shape[1],shape[0]).T; #back to original shape
F=F/np.sqrt(res);
return np.fft.fftshift(F);
# Inverse Fourier transform processor. Takes data in the fourier domain
# that has been properly normalized and arranged (i.e. negative intensities
# resulting from Python's fft algorithm flipped to the propper sign) by the
# companion routine FFT.
# First we fftshift the data, then we take the negative sign of every other
# mode in each direction (x and y), and finally we
# multiply the bad boy by (sqrt(res))^2 to get the correct intensity scaling.
from numpy.random import rand
from numpy import exp, sqrt, log, tan, pi, floor, zeros, ones
from scipy.special import gammaln
import numpy as np
def IFFT(F):
shape = F.shape;
res=max(shape);
F= np.fft.fftshift(F);
if shape[0] == shape[1]:
F[1:res:2,:]=-F[1:res:2,:]
F[:,1:res:2]=-F[:,1:res:2];
F=res*F;
f= np.fft.ifft2(F);
else:
F = F.ravel('F'); #create 1-d array, colum-major order (matlab style), not really nice
F[1:res:2] =-F[1:res:2];
F = F.reshape(shape[1],shape[0]).T; #back to original shape
F=np.sqrt(res)*F;
f= np.fft.ifft(F,axis=0);
return f;
# Poisson random number generation with mean of xm
# April 6, 2005 Ning Lei
# see numerical recipe C++ version
# ref: y(i+1) = xm^i*exp(-xm)/factorial(i); see Leo Breiman Probability
from numpy.random import rand
from numpy import exp, sqrt, log, tan, pi, floor, zeros, ones
from scipy.special import gammaln
import numpy as np
def PoissonRan(xm):
oldm = -1;
if xm<12:
if xm != oldm:
oldm = xm;
g = exp(-xm);
em = -1;
t = 1;
em = em+1;
t = t*rand(1);
while t>g :
em = em+1;
t = t*rand(1);
else:
if xm != oldm:
oldm = xm;
sq = sqrt(2.0*xm);
alxm = log(xm);
g = xm*alxm-gammaln(xm+1);
y = tan(pi*rand(1));
em = sq*y+xm;
while em < 0:
y = tan(pi*rand(1));
em = sq*y+xm;
em = floor(em);
t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g);
while rand(1) > t:
y = tan(pi*rand(1));
em = sq*y+xm;
while em < 0:
y = tan(pi*rand(1));
em = sq*y+xm;
em = floor(em);
t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g);
return em;
# resizes arrays. Takes an array A that is M by N and resizes it to Mnew by
# Nnew.
# usage: Anew = Resize(A, Mnew,Nnew)
from numpy.random import rand
from numpy import exp, sqrt, log, tan, pi, floor, zeros, ones
from scipy.special import gammaln
import numpy as np
def Resize(A, Mnew,Nnew):
M = A.shape[0];
N = A.shape[1];
Anew = zeros((Mnew,Nnew));
if Nnew <= N:
Nstep = int(N/Nnew);
Mstep = int(M/Mnew);
countN = -1;
for j in range(0,N,Nstep):
countN = countN+1;
countM = -1;
for i in range(0,M,Mstep):
countM = countM+1;
Anew[countM,countN] = A[i,j];
else:
Nstep = int(Nnew/N);
Mstep = int(Mnew/M);
temp=ones((Mstep,Nstep));
countN = -1;
for j in range(N):
countN = countN+1;
countM = -1;
for i in range(M):
countM = countM+1;
Anew[i*Mstep:(i+1)*Mstep,j*Nstep:(j+1)*Nstep] = A[countM,countN]*temp;
return Anew
...@@ -9,5 +9,9 @@ The "Utilities"-module contains various ... ...@@ -9,5 +9,9 @@ The "Utilities"-module contains various ...
from .Laplace_matrix import * from .Laplace_matrix import *
from .Procrustes import * from .Procrustes import *
from .FFT import *
from .IFFT import *
from .Resize import *
from .PoissonRan import *
__all__ = ["Laplace_matrix","procrustes"] __all__ = ["Laplace_matrix","procrustes", "FFT", "IFFT", "Resize", "PoissonRan"]
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment