aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xpython/msvst_starlet.py61
1 files changed, 40 insertions, 21 deletions
diff --git a/python/msvst_starlet.py b/python/msvst_starlet.py
index b90adf9..909cfac 100755
--- a/python/msvst_starlet.py
+++ b/python/msvst_starlet.py
@@ -11,9 +11,11 @@
#
# Aaron LI
# Created: 2016-03-17
-# Updated: 2016-04-20
+# Updated: 2016-04-22
#
# ChangeLog:
+# 2016-04-22:
+# * Show more verbose information/details
# 2016-04-20:
# * Add argparse and main() for scripting
#
@@ -25,8 +27,8 @@ And multi-scale variance stabling transform (MS-VST), which can be used
to effectively remove the Poisson noises.
"""
-__version__ = "0.2.0"
-__date__ = "2016-04-20"
+__version__ = "0.2.1"
+__date__ = "2016-04-22"
import sys
@@ -379,7 +381,7 @@ class IUWT_VST(IUWT): # {{{
data_ivst = data ** 2 + cb - self.vst_coef[scale]["c"]
return data_ivst
- def is_significant(self, scale, fdr=0.1, independent=False):
+ def is_significant(self, scale, fdr=0.1, independent=False, verbose=False):
"""
Multiple hypothesis testing with false discovery rate (FDR) control.
@@ -393,8 +395,8 @@ class IUWT_VST(IUWT): # {{{
https://en.wikipedia.org/wiki/False_discovery_rate
"""
coef = self.get_detail(scale)
- pvalues = 2 * (1 - sp.stats.norm.cdf(np.abs(coef) / \
- self.vst_coef[scale]["std"]))
+ std = self.vst_coef[scale]["std"]
+ pvalues = 2.0 * (1.0 - sp.stats.norm.cdf(np.abs(coef) / std))
p_sorted = pvalues.flatten()
p_sorted.sort()
N = len(p_sorted)
@@ -406,9 +408,13 @@ class IUWT_VST(IUWT): # {{{
comp = (p_sorted < p_comp)
# cutoff p-value after FDR control/correction
p_cutoff = np.max(p_sorted[comp])
+ if verbose:
+ print("std/sigma: %g, p_cutoff: %g" % (std, p_cutoff),
+ flush=True, file=sys.stderr)
return (pvalues <= p_cutoff, p_cutoff)
- def denoise(self, fdr=0.1, fdr_independent=False):
+ def denoise(self, fdr=0.1, fdr_independent=False, start_scale=1,
+ verbose=False):
"""
Denoise the wavelet coefficients by controlling FDR.
"""
@@ -418,17 +424,29 @@ class IUWT_VST(IUWT): # {{{
# supports of significant coefficients of each scale
self.sig_supports = [None] # make index match the scale
self.p_cutoff = [None]
+ if verbose:
+ print("MSVST decomposing ...", flush=True, file=sys.stderr)
for scale in range(1, self.level+1):
coef = self.get_detail(scale)
- sig, p_cutoff = self.is_significant(scale, fdr, fdr_independent)
- coef[np.logical_not(sig)] = 0.0
+ if verbose:
+ print("\tScale %d: " % scale, end="",
+ flush=True, file=sys.stderr)
+ if scale < start_scale:
+ if verbose:
+ print("skipped", flush=True, file=sys.stderr)
+ sig, p_cutoff = None, None
+ else:
+ sig, p_cutoff = self.is_significant(scale, fdr=fdr,
+ independent=fdr_independent, verbose=verbose)
+ coef[np.logical_not(sig)] = 0.0
+ #
self.denoised.append(coef)
self.sig_supports.append(sig)
self.p_cutoff.append(p_cutoff)
# append the last approximation
self.denoised.append(self.get_approx())
- def decompose(self, level, boundary="symm"):
+ def decompose(self, level=5, boundary="symm", verbose=False):
"""
2D IUWT decomposition with VST.
"""
@@ -439,7 +457,12 @@ class IUWT_VST(IUWT): # {{{
self.calc_vst_coef()
self.decomposition = []
approx = self.data
+ if verbose:
+ print("IUWT decomposing (%d levels): " % level,
+ end="", flush=True, file=sys.stderr)
for scale in range(1, level+1):
+ if verbose:
+ print("%d..." % scale, end="", flush=True, file=sys.stderr)
# approximation:
approx2 = signal.convolve2d(self.data, self.filters[scale],
mode="same", boundary=self.boundary)
@@ -449,6 +472,8 @@ class IUWT_VST(IUWT): # {{{
if scale == level:
self.decomposition.append(approx2)
approx = approx2
+ if verbose:
+ print("DONE!", flush=True, file=sys.stderr)
return self.decomposition
def reconstruct_ivst(self, denoised=True, positive_project=True):
@@ -476,7 +501,7 @@ class IUWT_VST(IUWT): # {{{
self.reconstruction = reconstruction
return reconstruction
- def reconstruct(self, denoised=True, niter=20, verbose=False):
+ def reconstruct(self, denoised=True, niter=10, verbose=False):
"""
Reconstruct the original image using iterative method with
L1 regularization, because the denoising violates the exact inverse
@@ -552,7 +577,7 @@ def main():
action="store_true", default=False,
help="whether the FDR null hypotheses are independent")
parser.add_argument("-n", "--niter", dest="niter",
- type=int, default=20,
+ type=int, default=10,
help="number of iterations for reconstruction")
parser.add_argument("-v", "--verbose", dest="verbose",
action="store_true", default=False,
@@ -576,15 +601,9 @@ def main():
img = imgfits[0].data
# Remove Poisson noises
msvst = IUWT_VST(data=img)
- if args.verbose:
- print("INFO: MSVST decomposing ...", file=sys.stderr)
- msvst.decompose(level=args.level)
- if args.verbose:
- print("INFO: MSVST denosing ...", file=sys.stderr)
- msvst.denoise(fdr=args.fdr, fdr_independent=args.fdr_independent)
- if args.verbose:
- print("INFO: MSVST reconstructing (this may take a while) ...",
- file=sys.stderr)
+ msvst.decompose(level=args.level, verbose=args.verbose)
+ msvst.denoise(fdr=args.fdr, fdr_independent=args.fdr_independent,
+ verbose=args.verbose)
msvst.reconstruct(denoised=True, niter=args.niter, verbose=args.verbose)
img_denoised = msvst.reconstruction
# Output