Commit ddddb453 authored by Alessandro Mirone's avatar Alessandro Mirone
Browse files

Merge branch 'chambolle_pock_tv_positive' into alessandro

Conflicts:
	XRStools/roiSelectionWidget.py
parent 69a91b34
......@@ -2366,7 +2366,7 @@ def superR_getVolume(mydata):
else:
debin = mydata['debin']
h5f = h5py.File(scalprods_filename)
h5f = h5py.File(scalprods_filename, "r")
h5 = h5f [scalprods_groupname]
scalDS = h5["scalDS"][:]
......
......@@ -81,7 +81,7 @@ class MaskImageWidget(sole_MaskImageWidget.MaskImageWidget):
if Ct:
print( str(e.mimeData().text()))
Cc = string.atoi( str(e.mimeData().text()))
Cc = int( str(e.mimeData().text()))
# bytearray = e.mimeData().data('application/x-qabstractitemmodeldatalist')
# data = decode_data(bytearray)
# print data
......
......@@ -81,7 +81,7 @@ class MaskImageWidget(sole_MaskImageWidget.MaskImageWidget):
if Ct:
print( str(e.mimeData().text()))
Cc = string.atoi( str(e.mimeData().text()))
Cc = int( str(e.mimeData().text()))
# bytearray = e.mimeData().data('application/x-qabstractitemmodeldatalist')
# data = decode_data(bytearray)
# print data
......
......@@ -1093,7 +1093,7 @@ class roi_finder:
self.roi_obj.masks = xrs_rois.convert_roi_matrix_to_masks(self.roi_obj.roi_matrix)
self.roi_obj.number_of_rois = int(np.amax(self.roi_obj.roi_matrix))
def refine_rois_PW(self, hydra_obj, scan_numbers):
def refine_rois_PW(self, hydra_obj, scan_numbers,save_dataset=None):
""" **refine_rois_PW**
Allows for manual refinement of ROIs by plotting the spectra column-wise.
......@@ -1121,6 +1121,9 @@ class roi_finder:
counter = 0
for data, key in zip(cw_data, sorted(self.roi_obj.red_rois)):
for ii in range(data.shape[1]):
if save_dataset:
the_shelve = shelve.open(save_dataset)
the_key = str('%s_PixelNo%d'%(key,ii) )
plt.cla()
title_txt = 'Click above to keep/below black line to discard column for ROI %02d'%(counter+1) + '.'
plt.title(title_txt)
......@@ -1134,12 +1137,18 @@ class roi_finder:
# let user click a point on figure
user_input = np.array(plt.ginput(1,timeout=-1)[0])
if user_input[1] >= offset:
if save_dataset:
the_shelve[the_key] = { 'spectrum': data[:,ii], 'decision': 1 }
pass
elif user_input[1] < offset:
index = np.unravel_index([ii],self.roi_obj.red_rois[key][1].shape )
self.roi_obj.red_rois[key][1][index[0],index[1]] = 0
if save_dataset:
the_shelve[the_key] = { 'spectrum': data[:,ii], 'decision': 0 }
else:
print('Something fishy happened!')
if save_dataset:
the_shelve.close()
counter += 1
# reassign ROI object
self.roi_obj.roi_matrix = xrs_rois.convert_redmatrix_to_matrix(self.roi_obj.red_rois, np.zeros(self.roi_obj.input_image.shape))
......
......@@ -11,6 +11,178 @@ try:
except:
print( " ATTENTION : SKIMAGE RESORATION NOT LOADED ")
def mdiv(grad):
res = np.zeros(grad.shape[1:])
res[ :-1 , :, : ] += this_grad[0, :-1 , :,:]
res[ 1:-1 , :, : ] -= this_grad[0, :-2 , :,:]
res[ -1 , :, : ] -= this_grad[0, -2 , :,:]
res[ :, :-1 , : ] += this_grad[1, :, :-1 ,:]
res[ :, 1:-1 , : ] -= this_grad[1, :, :-2 ,:]
res[ :, -1 , : ] -= this_grad[1, :, -2 ,:]
res[ :, : , :-1 ] += this_grad[1, : ,:, :-1 ]
res[ :, : , 1:-1 ] -= this_grad[1, : ,:, :-2 ]
res[ :, : , -1 ] -= this_grad[1, : ,:, -2 ]
return res
def mygradient(img):
shape = [3 ] + list(img.shape)
gradient = np.zeros(shape, dtype=img.dtype)
gradient[0,:,:,: ] = np.diff(img, axis=0)
gradient[1,:,:,: ] = np.diff(img, axis=1)
gradient[2,:,:,: ] = np.diff(img, axis=2)
return gradient
def v_project(v,weight ):
norms = np.minimum( weight, np.sqrt( v[0]*v[0] + v[1]*v[1] ))
return v/ norms
def my_denoise_tv_chambolle_positive(image, weight=0.1, n_iter_max=200):
ndim = image.ndim
g = np.zeros_like(p)
x = np.zeros_like(image)
tmpxa = np.zeros_like(image)
v = np.zeros((image.ndim, ) + image.shape, dtype=image.dtype)
i = 0
sigma = 1.0/math.sqrt(8.0)
tau = 1.0/math.sqrt(8.0)
while i < n_iter_max:
tmpxa[:] = x + sigma * ( ( image-x) + mydiv( v ) )
tmpxa[:] = np.maximum (tmpxa)
tmpxa[:] = tmpxa-x
x[:] = x + tmpxa
tmpxa[:] = x + tmpxa
v[:] = v + tau * mygrad(tmpxa)
v = v_project(v,weight )
return x
def _denoise_tv_chambolle_nd(image, weight=0.1, eps=2.e-4, n_iter_max=200,
positivity=False):
"""Perform total-variation denoising on n-dimensional images.
Parameters
----------
image : ndarray
n-D input data to be denoised.
weight : float, optional
Denoising weight. The greater `weight`, the more denoising (at
the expense of fidelity to `input`).
eps : float, optional
Relative difference of the value of the cost function that determines
the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max : int, optional
Maximal number of iterations used for the optimization.
positivity : bool, optional
Adds positivity constraint
Returns
-------
out : ndarray
Denoised array of floats.
Notes
-----
Rudin, Osher and Fatemi algorithm.
LICENCE
-------
Copyright (C) 2011, the scikit-image team
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the
distribution.
3. Neither the name of skimage nor the names of its contributors may be
used to endorse or promote products derived from this software without
specific prior written permission.
this software is provided by the author ``as is'' and any express or
implied warranties, including, but not limited to, the implied
warranties of merchantability and fitness for a particular purpose are
disclaimed. in no event shall the author be liable for any direct,
indirect, incidental, special, exemplary, or consequential damages
(including, but not limited to, procurement of substitute goods or
services; loss of use, data, or profits; or business interruption)
however caused and on any theory of liability, whether in contract,
strict liability, or tort (including negligence or otherwise) arising
in any way out of the use of this software, even if advised of the
possibility of such damage.
"""
ndim = image.ndim
p = np.zeros((image.ndim, ) + image.shape, dtype=image.dtype)
g = np.zeros_like(p)
d = np.zeros_like(image)
i = 0
while i < n_iter_max:
if i > 0:
# d will be the (negative) divergence of p
d = -p.sum(0)
slices_d = [slice(None), ] * ndim
slices_p = [slice(None), ] * (ndim + 1)
for ax in range(ndim):
slices_d[ax] = slice(1, None)
slices_p[ax+1] = slice(0, -1)
slices_p[0] = ax
d[tuple(slices_d)] += p[tuple(slices_p)]
slices_d[ax] = slice(None)
slices_p[ax+1] = slice(None)
out_nopos = image + d
else:
out_nopos = image
if not positivity:
out = out_nopos
else:
out = np.maximum(0, out_nopos)
removed = np.minimum(out_nopos, 0)
d = d-removed
E = (d ** 2).sum()
# g stores the gradients of out along each axis
# e.g. g[0] is the first order finite difference along axis 0
slices_g = [slice(None), ] * (ndim + 1)
for ax in range(ndim):
slices_g[ax+1] = slice(0, -1)
slices_g[0] = ax
g[tuple(slices_g)] = np.diff(out, axis=ax)
slices_g[ax+1] = slice(None)
norm = np.sqrt((g ** 2).sum(axis=0))[np.newaxis, ...]
E += weight * norm.sum()
tau = 1. / (2.*ndim)
norm *= tau / weight
norm += 1.
p -= tau * g
p /= norm
E /= float(image.size)
if i == 0:
E_init = E
E_previous = E
else:
if np.abs(E_previous - E) < eps * E_init:
break
else:
E_previous = E
i += 1
return out
def superr( scalDD, scalDS, scalSS, niter=15, beta=1.0e-8):
"""
- scalDS which is an array [ZDIM,YDIM,XDIM] , type "d" .
......@@ -68,7 +240,9 @@ def Fista( scalDD, scalDS , scalSS, solution , niter=500, beta=0.1 ):
err = calculate_grad(scalDD, scalDS , scalSS, y, grad)
solution[:] = y - grad/Lip
solution[:]=skimage.restoration.denoise_tv_chambolle(solution, weight=beta, eps=0.000002)
# solution[:]=skimage.restoration.denoise_tv_chambolle(solution, weight=beta, eps=0.000002)
solution[:]=_denoise_tv_chambolle_nd(solution, weight=beta, eps=0.000002, positivity=True)
## solution[:] = np.maximum(solution, 0)
......
......@@ -780,7 +780,9 @@ class xyzBox:
# get a test box
pbcMols = getPeriodicTestBox_molecules(self.xyzMolecules,self.boxLength,numbershells=1)
# cut clusters and write files
for o_atom, ii in zip(o_atoms,range(len(o_atoms))):
for ii, mol in enumerate(self.xyzMolecules):
o_atom = mol.get_atoms_by_name(o_name)[0]
#for o_atom, ii in zip(o_atoms,range(len(o_atoms))):
cluster = []
for molecule in pbcMols:
coor = molecule.getCoordinates_name(o_name)
......@@ -999,6 +1001,19 @@ class xyzBox:
else:
return np.arccos( np.clip( dotp, -1.0, 1.0 ) )
def get_angle( self, atom1, atom2, atom3, degrees=True):
""" **get_angle**
Return angle between the three given atoms (as seen from atom2).
"""
vec1 = getDistVector(atom1, atom2)
vec2 = getDistVector(atom3, atom2)
dotp = np.dot(vec1/np.linalg.norm(vec1), vec2/np.linalg.norm(vec2))
if degrees:
return np.degrees( np.arccos( np.clip( dotp, -1.0, 1.0 ) ) )
else:
return np.arccos( np.clip( dotp, -1.0, 1.0 ) )
def count_neighbors( self, name1, name2, cutoff_low=0.0, cutoff_high=2.0, counter_name='num_OO_shell' ):
""" **count_neighbors**
......@@ -1304,7 +1319,7 @@ def getPeriodicTestBox_molecules(Molecules,boxLength,numbershells=1):
cpAtom.translateSelf(vector*boxLength)
cpAtoms.append(cpAtom)
pbc_molecules.append(xyzMolecule(cpAtoms))
return pbc_moleculesfor
return pbc_molecules
def getPeriodicTestBox(xyzAtoms,boxLength,numbershells=1):
vectors = []
......@@ -2205,9 +2220,9 @@ def writeOCEANinput(fname,headerfile,xyzBox,exatom,edge,subshell):
if atom.name == 'C':
inputf.write('%d ' %4)
if atom.name == 'Cl':
inputf.write('%d ' %5)
if atom.name == 'Na':
inputf.write('%d ' %3)
if atom.name == 'Na':
inputf.write('%d ' %4)
inputf.write('\n')
inputf.write('} \n')
inputf.write('\n')
......@@ -2257,17 +2272,17 @@ def writeOCEANinput_new(fname,headerfile,xyzBox,exatom,edge,subshell):
inputf.write('typat { \n')
for atom in xyzBox.xyzAtoms:
if atom.name == 'H':
inputf.write('%d ' %4)
if atom.name == 'O':
inputf.write('%d ' %1)
if atom.name == 'N':
if atom.name == 'O':
inputf.write('%d ' %2)
if atom.name == 'N':
inputf.write('%d ' %6)
if atom.name == 'C':
inputf.write('%d ' %3)
if atom.name == 'Cl':
inputf.write('%d ' %5)
if atom.name == 'Na':
if atom.name == 'Cl':
inputf.write('%d ' %3)
if atom.name == 'Na':
inputf.write('%d ' %4)
inputf.write('\n')
inputf.write('} \n')
inputf.write('\n')
......@@ -2283,7 +2298,7 @@ def writeOCEANinput_new(fname,headerfile,xyzBox,exatom,edge,subshell):
#ind = 1
#for atom in xyzBox.xyzAtoms:
# if atom.name == exatom:
inputf.write('%d %d %d \n' %(-exatom,edge,subshell))
inputf.write('%d %d %d \n' %(-xrs_utilities.element(exatom),edge,subshell))
# ind += 1
inputf.write('} \n')
inputf.write('\n')
......
......@@ -1035,7 +1035,7 @@ class Hydra:
# make sure scan exists
for scannum in scannums:
if not scannum in self.scan_numbers:
self.load_scan(scannum, direct=False, method=method)
self.load_scan(scannum, direct=True, method=method)
# make sure raw_signals is existing
if not self.scans['Scan%03d' % scannums[0]].raw_signals:
......
......@@ -721,7 +721,7 @@ class Scan:
if not interp:
for key in self.raw_signals:
if method in [ 'sum', 'pixel', 'pixel2', 'row', 'column']:
if method in [ 'sum', 'row', 'column']:
# recover raw counts
signals1 = self.raw_signals[key]*self.monitor
......@@ -739,7 +739,27 @@ class Scan:
# assign and renormalize
self.raw_signals[key] = av_signals/av_monitor
self.raw_errors[key] = av_errors/av_monitor
self.raw_errors[key] = av_errors/av_monitor
if method in [ 'pixel', 'pixel2']:
# recover raw counts
signals1 = self.raw_signals[key]*self.monitor[:,None,None]
signals2 = scan.raw_signals[key]*scan.monitor[:,None,None]
# sum signals
av_signals = np.sum((signals1, signals2), axis=0)
# sum monitors
av_monitor = np.sum((self.monitor, scan.monitor), axis=0)
# errors of summed signals
av_errors = np.sqrt(av_signals)
# assign and renormalize
self.raw_signals[key] = av_signals/av_monitor[:,None,None]
self.raw_errors[key] = av_errors/av_monitor[:,None,None]
if interp:
#rbfi = Rbf( scan.energy, scan.monitor, function='linear' )
......
......@@ -331,7 +331,7 @@ def find_diag_angles(q, x0, U, B, Lab, beam_in, lambdai, lambdao, tol=1e-8, meth
# least square minimization between wanted and guessed q
fitfctn = lambda x: np.sum(( q - get_UB_Q(x[0], x[1], x[2], x[3], x[4], \
**kwargs)[0] )**2)
ans=minimize(fitfctn, x0, bounds=((20.,80.),(-10.,110.),(-7.,7.),(-7.,7.),(None,None)), tol=tol, method=method)
ans=minimize(fitfctn, x0, bounds=((0.,0.),(-10.,110.),(-7.,7.),(-7.,7.),(None,None)), tol=tol, method=method)
print( ans )
return ans.x
......@@ -3608,6 +3608,10 @@ def hlike_Rwfn(n,l,r,Z):
return factor1*factor2*factor3*lag#*np.sqrt(n+1.0)
def compute_matrix_elements(R1,R2,k,r):
# for ii=1:length(q);
# fun=y3d.^2.*besselj(4,q(ii)*r);
# int4(ii)=simpson(r,fun);
# end
from scipy import special
q = np.linspace(0,30,len(r))
r2RsphBR = np.linspace(0,30,len(r))
......
Supports Markdown
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