From b9b5a09f044c339b4de42095bfdedc79f2e4902b Mon Sep 17 00:00:00 2001 From: tjof2 Date: Mon, 16 May 2016 08:05:25 +0100 Subject: [PATCH] Updated example and added OMP options --- examples/PGURE-SVT HyperSpy Demo.ipynb | 111 +++++++------- examples/mixed_noise_figure.png | Bin 0 -> 5007 bytes pguresvt/pguresvt.py.in | 191 +++++++++++++++++-------- src/lib-PGURE-SVT.cpp | 151 ++++++++----------- 4 files changed, 247 insertions(+), 206 deletions(-) create mode 100644 examples/mixed_noise_figure.png diff --git a/examples/PGURE-SVT HyperSpy Demo.ipynb b/examples/PGURE-SVT HyperSpy Demo.ipynb index 81705e7..b26a0ce 100644 --- a/examples/PGURE-SVT HyperSpy Demo.ipynb +++ b/examples/PGURE-SVT HyperSpy Demo.ipynb @@ -4,12 +4,20 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# PGURE-SVT HyperSpy Demonstration\n", + "# PGURE-SVT Demonstration\n", "\n", "### Tom Furnival ([tjof2@cam.ac.uk](mailto:tjof2@cam.ac.uk))\n", "\n", "PGURE-SVT is an algorithm designed to denoise image sequences acquired in microscopy. It exploits the correlations between consecutive frames to form low-rank matrices, which are then recovered using a technique known as nuclear norm minimization. An unbiased risk estimator for mixed Poisson-Gaussian noise is used to automate the selection of the regularization parameter, while robust noise and motion estimation maintain broad applicability to many different types of microscopy.\n", "\n", + "You can read more about the algorithm and its applications in:\n", + "\n", + "> T Furnival, RK Leary, PA Midgley _\"Denoising time-resolved microscopy sequences with singular value thresholding.\"_ (Article in press). DOI:[10.1016/j.ultramic.2016.05.005](http://dx.doi.org/10.1016/j.ultramic.2016.05.005)\n", + "\n", + "The source code is available to download from [http://tjof2.github.io/pgure-svt/](http://tjof2.github.io/pgure-svt/).\n", + "\n", + "### Use with HyperSpy\n", + "\n", "This example notebook shows how PGURE-SVT can be combined with [HyperSpy](http://hyperspy.org), which is an open-source Python library that makes signal handling and processing straightforward in Python, with a friendly API. While you can use `pguresvt.pguresvt.SVT` to denoise a numpy array directly, `pguresvt.hspysvt.HSPYSVT` can directly denoise a HyperSpy signal." ] }, @@ -22,33 +30,15 @@ "outputs": [], "source": [ "# Configures the plotting backend\n", - "#%matplotlib inline\n", - "%matplotlib qt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ + "%matplotlib qt\n", + "\n", + "# Import NumPy and the HyperSpy API\n", "import numpy as np\n", - "# Import the HyperSpy API\n", - "import hyperspy.api as hs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ + "import hyperspy.api as hs\n", + "\n", "# Import the HyperSpy wrapper for PGURE-SVT\n", - "from pguresvt import hspysvt" + "from pguresvt import hspysvt\n", + "from pguresvt.pguresvt import PoissonGaussianNoiseGenerator" ] }, { @@ -80,7 +70,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we corrupt it with a mixture of Poisson and Gaussian noise." + "Now we corrupt the movie with using PGURE-SVT's built-in noise generator for mixed Poisson-Gaussian noise, according to the equation:\n", + "\n", + "where the parameters are defined as:\n", + "```\n", + "alpha = detector gain\n", + "mu = detector offset\n", + "sigma = detector noise\n", + "```" ] }, { @@ -91,25 +88,19 @@ }, "outputs": [], "source": [ - "# First extract the data and rescale to [0,1] range\n", - "clean = movie._data_aligned_with_axes\n", - "clean = clean / np.amax(clean)\n", - "\n", "# Detector gain\n", - "gain = 0.1\n", + "detector_gain = 0.1\n", "# Detector offset\n", - "offset = 0.1\n", - "# Detector variance\n", - "sigma = 0.1\n", - "\n", - "def addnoise(x):\n", - " return gain * np.random.poisson(x / gain) + offset + sigma * np.random.randn()\n", - "addnoise = np.vectorize(addnoise, otypes=[np.float])\n", - "\n", - "noisy = addnoise(clean)\n", - "\n", + "detector_offset = 0.1\n", + "# Detector noise\n", + "detector_sigma = 0.1\n", + "\n", + "noisy = PoissonGaussianNoiseGenerator(movie._data_aligned_with_axes,\n", + " alpha=detector_gain,\n", + " mu=detector_offset,\n", + " sigma=detector_sigma)\n", "noisy_movie = hs.signals.Image(noisy)\n", - "noisy_movie.plot()" + "noisy_movie.plot(navigator='slider')" ] }, { @@ -131,10 +122,13 @@ " arpssize=7, \n", " tol=1e-7,\n", " median=5,\n", - " hotpixelthreshold=10)\n", + " hotpixelthreshold=10,\n", + " numthreads=1)\n", "```\n", "\n", - "In this example we do not use the noise estimation procedure, and instead provide the known parameters to the algorithm directly. This information is used by the PGURE optimizer to calculate the threshold." + "In this example we do not use the noise estimation procedure, and instead provide the known parameters to the algorithm directly. This information is used by the PGURE optimizer to calculate the threshold.\n", + "\n", + "*Note:* If you have a multicore machine, you can set `numthreads > 1` to speed up the calculation." ] }, { @@ -149,10 +143,11 @@ "# Initialize with suggested parameters\n", "svt = hspysvt.HSPYSVT(patchsize=4,\n", " estimatenoise=False,\n", - " alpha=gain,\n", - " mu=offset,\n", - " sigma=sigma,\n", - " tol=1e-6)" + " alpha=detector_gain,\n", + " mu=detector_offset,\n", + " sigma=detector_sigma,\n", + " tol=1e-5,\n", + " numthreads=2)" ] }, { @@ -174,16 +169,18 @@ "denoised_movie = svt.denoise(noisy_movie)\n", "\n", "# Plot denoised data\n", - "denoised_movie.plot()" + "denoised_movie.plot(navigator='slider')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 2. Time-resolved ADF-STEM image sequence\n", + "## 2. Time-resolved ADF-STEM image sequence\n", "\n", - "In this example we apply PGURE-SVT to an experimental dataset of a nanoparticle acquired using ADF-STEM. Here the noise levels are not known, so a noise estimation procedure is used." + "In this example we apply PGURE-SVT to an experimental dataset of a nanoparticle acquired using ADF-STEM. For larger images, you can use the `patchoverlap` parameter to control the trade-off between speed and accuracy of the denoising procedure.\n", + "\n", + "Here only the detector offset (`mu`) was known beforehand, so a noise estimation procedure is used for the other values." ] }, { @@ -194,19 +191,21 @@ }, "outputs": [], "source": [ - "# Load example dataset\n", + "# Load example dataset and plot\n", "expt_movie = hs.load(\"experimentalsequence.tif\")\n", + "expt_movie.plot(navigator='slider')\n", "\n", "# Initialize with suggested parameters\n", "expt_svt = hspysvt.HSPYSVT(patchsize=4,\n", - " patchoverlap=2, \n", - " tol=1e-6)\n", + " patchoverlap=2,\n", + " mu=0.1,\n", + " numthreads=2)\n", "\n", "# Run the denoising\n", "denoised_movie = expt_svt.denoise(expt_movie)\n", "\n", "# Plot denoised data\n", - "denoised_movie.plot()" + "denoised_movie.plot(navigator='slider')" ] } ], diff --git a/examples/mixed_noise_figure.png b/examples/mixed_noise_figure.png new file mode 100644 index 0000000000000000000000000000000000000000..ec2fc06998a5add24c905215a5f50e6e4d8a1ca1 GIT binary patch literal 5007 zcmbtYivE+d&~${duZFw zaSq3m!y~+t?Jd?a5=Y~H-F{h0U&ZReTFjrj3_ShGj*6{JXV-ky&objB{-;oJJ`u%g zmh^6(4We&)ZaXMXk&?=k=L>VAx@YYhVDhBQhm*CjsoVEJ$uW!|zhQ0i+k13rXx;%a3~!Aki8Rvmk5qeSMSz9hjP3M~~)+cD}@Ro3+dPMFFnHT;PlbqOV)E5ZG5a zU?qHPjm{#9_N(eGj?GL6Y5ZW$_evgwB9}i9nU>n(lrFi84Tv&;_2E?bv}wCwR1cA~ zv`}N{?VI~=yK@^f+LL@F| zg=08AeiONY)_-X;HCYuYB`dHw%A|2~?uSN6#FVDzs8H`dLEqPWcVBwN>U;g#TauHU zl8<#?PZSS@xVle|;IvI^NF>J2tjI0~F0Fk+ZRu;fr}VvqJ9n#HhqkTAB8?JoAT%9 zq@{|SX{Eebg2P3#T8c+aOSRld2~Dm~jKV0V7R-^WeRF=#_BAJ#Pr-Dj+ERyHOu@qF zjSHVnKfU~qKReSC6HMZU&<3BgG%=wBe((xUK_9coF~qb<+R5WWU2WXw!BuyedgNNr zz<{=VA^CV#(r3g$RTK|DA5tfg>Mw5m%N0@LD%WqruJKL%d-X#kk2N}D@etebx7tdwR zRn8!>wg_3~GVVmX;0U2lro|Zw#Kdl6W7B3wrx+5_rm$zp{Y(ZA30h06f#0@{UVG_V zq6waDm&RI_qFi2cdw#0_s)UStH)D8jr+mk`>`s?JYv=}O0&$wE)$bddBdW2MSW>+4bs-(W1YP@PJ zbenCzif@Lkp`^j9cC*WQxrT(jSxg39kV;@6S*sqP)RRFHLEA3Hb@b^K0fZ_UfsZ?s zjisCb)}OYL5aWahTpcUkk$f)GXs|HXuKM;gNMeM3ELKgF?lDZ2CzKmsqs2%k3B}dk zEKq1VH<^0_fRc}qVw5a3V|Ep35?}iKsTO_xD6@cL{1Cx#b&i<#L?*R)^IP2#$)*Bs zp*aO(ukJfX?bqnaUCH;NKW3x9vxKt$ZVTH&2(g0LhZr`qAK}{vd(b)c17N~=RKS3epu~*9A zhS-@+sk@;@{)C7aB?wA6*yUHHY5FWoF{h~0;BU&H*}ld`0WauHRyRS^!QBF3Lz*NU z1?79I;6)KNtX2p;x*|%b-4rHli2_4 zFz1F~f{GS>Y?hYm`@zc(p5z!Fzsk707l|He)bUIs{uNy%+>kRiN?@YXq6kPMUi(3B zHtA!X@b~3x;hrucFx2u+pH@YpJ`fGd`k9GLp&t=m*%R3WRp>iF*$iMWCtuszO{t{k zN?p3QptsXG^%?zFoi+Fps{6zo_O~o!d!dO0%2o29Ac`J7gUu4$=l+$mjoOSNhJzq_ zz*PO*!Dv?6)c#SP0FVrY>$pmAEDNVvmU=LUq!r71V{&d33`xX|8MJ5gH~!G@deJG( zlILt>h3@;BLy2b5Xy+rO01j4b=Xqt0(IiSB(JLABcUgV33&tK19RTv+wBy!En9^|O`! z6E68W`tg@ZsX3e9vTQ!4=KL~%t=5%}P)TC09vr>qep>eat>w`%&2PoeOuOji_9=29 zf1l?<>(@4`O3i-{9G>JAexb@^wSip+o$H#^zds=?L%@oQ9>+pF1MO&6mX(&PHpjm= z+#M6M4(S!KTMUMSD=CJQ+II~J?JWDv4oM&KA&FErx}K2x=Ja*{>bS^=_>+>;J|a=q z>5rg0d<7vZY_&n9i3v1x?{vZEap?$>+tQkztEa2Y8D%Dwsa`lmQfFDaP_sz-xcBQK+X-exd84!i*__#gfIoa%A_d0 z&}hs`Ig)8q!0zgUxmW9(AdP-+G6x-0v*aGSB8eJ9sXs=T0p(!w*RgqFWsu5kRm9UN zXuv*|SmecYR%E0kZs@DvsD{9tg5o&mC+3X--{;l2y2Fs}{wyR2;Gv!9LH4Q1e?-R! zR;3xuD$^m;H&V51fOXJdZIPAAvCQ?rWYuaiG!7)&oC`o*_(khP4&7iE_5Jm&s&6aN zfJ=+{Q~8GWAkrLmR!`~P(d$sfyG+_tedBkG*&{tZS;S|$7|Id5j4{OZGQWRyo%bm~ z$Aj06nqYNl$&rYyY2~@z$PNqPCifzw1yJ71hsD&Y%-2VUL~!jNPtP9x8ijpYhXA@A zcZZxhnx|{4tu&q*<*7T_qH*u?aq8bCVEpYL;uN8M%#i)j{+Ra9=0AOs;V~2u-{0>x z(9P>7CqZF##%&+QI{H`%GdUSq1cAMuF-Tf}*!Xm4`tdr=yV5lBID@Xl`#dc5w|WZ4 z`x>YJj2DQ0(T=YV&rUo=zX&tuPjOxu2tu&Ow_wKUgAU=yVH38Et0Ew6?I?s6#|=U zqma(lpfRi)wt$j8%WZ6yXn%cTh`~N!EKjnAHLzO1y^R?TUGtr}!1BN0q30&gEsI@+ z#qlqd3r*MXdIYs|PyKZNa3=YibWmf!LkzYhjq;#3r&dNX1Jrjjjppjt24ao~>oatD zTETuwm1O+k;yc)CrTLrRZ+ij@1P^{ZAKTGMUrw)7_gm4Zy!^f%72z4Q0$kS|>|8jp zu*w4HL!XWrja$hGygVRXvmRc{0A>IF+-+Z6CiHSdFp5_4dg50|cAyeEl2Ht<9&o~Z zt9g@Ix+Nms9TGPzzFUOSVQt_hDw-%gKz!#Qz+DROoYahemu*S2x6t@<%CEB{lYWwP z%CK$`$U%0TTQv>m^?puaIz{!&&x<*=vF>$5z zVf!AvYx*Gu`AYwgkTcRcTPklgU1fbrRc?C8j2qa3nAw)u65P#uvzRCSZ6vLtCLnrE z^uY>7BX9C2^9Vz;42Kb8;*QRX*2aUdY zaI1uoJ4prBG@>pag=NFllS=gJI*merKZalFd zN9J3W{G%IM=CHaySDz&`dduvrmko|fkfo$om0K1mKixU$muHZ9tstGk@;#;ightN% z)M6}7&vuuys&Ti;?@%w0?6jl%zMnybkNiKljL8iNhERd#e`Iz24q~|AkCUzJ0Wo%d z@m1^&(T-%D-b-aRfgYL8uK8y6sEMRsWZsXk2GY5s&XPrHPn|u!3EZ%m&n^$|mU*>d z+CXjZPa+i({C*B{Wa~5LNao8M|9_q%PiwQhxas{Z!Osh_HWjSG zsRgd_p1&>EtbhHtQ91yZwYmZjF)gUODWG<*!bnGDZBN# zn4$KB!NG1+eey}V$a6faP_IbjMN=f}SR^&QePu8ti?Agsx@gl~bZF{L_!vV_LFo6H z^%vKwOtbhkr-}L9Ly>-yinle)|?b@J?6_l}?E>DoA#WfIN1t@_OK zZoi~*-m5~~l;wOZMZ``_!YSKiXpTqNdnq{kC(_$+)QkqktM9NXJ;HU!))v58134$K zw6Ba%GaBbKHTsHl`dJoiy(IOI$9cpV(1XmwV2?}yC&vhsK2|k-@^y8kBi4CRA=rpF zBfW`fU!Lvl@^!h1=`Jg2)!6ImBD(SS9Q~y0cM`v4yR|;&gB|kTev`?qkRa(=rx>5> z{1W5c+|MKk1E#kKb6F|B+r0C!aYci$6jQ6$`Qn9;*!fAl_1t9YSw5=6)57{yh8Ims z7{E$t=4+<(t(J=N!PQVDOt_+lm_N#+JTVPj<8xY&O{TC_rPAf^BVVxpK|7>}cp-<+ zF8fH8ox7a^N1%$g<1z)bca~^pHsW<*)m2PkW4NEtbE^~#O#c?zWed)ng^Qt-249KO zot%^I(#UC;`qZHLW^ezfWC~#eD6=gMEF`Nfx=xODi4*Xxv5!ZavP&J~QklfGqyi?9`Z55ucMNRvv*HTfw;r^%^dj%w zsSh=KtTqX7%&!6C=D<3|D4SQxc^kDY-i~5|PmutMA7K{sZ zGhitAna777L&>8Ib@%LFdVYeGJdPdPS<9Ab330Vyc*DalqpFPkYz4~j(lU@)A`9if zGxco5it^BQEQg|@PNELtTjS(DGekd{xVs3O+Ge_V8E)yZ9ic_e1Aa>A`X~hMWir&J z>CD0dCRR1XryI#5JZIt{212MJ?%eyYmG90DZpf>UTZ@u)xPAO58UUi9uU?~SANe26 C{+QAL literal 0 HcmV?d00001 diff --git a/pguresvt/pguresvt.py.in b/pguresvt/pguresvt.py.in index 321209d..5846aad 100644 --- a/pguresvt/pguresvt.py.in +++ b/pguresvt/pguresvt.py.in @@ -34,98 +34,159 @@ # You should have received a copy of the GNU General Public License # along with PGURE-SVT. If not, see . -import ctypes, os +import ctypes, os, multiprocessing import numpy as np from numpy.ctypeslib import ndpointer -def is_power_of_two(n): - n = n/2 - if n == 2: - return True - elif n > 2: - return is_power_of_two(n) - else: - return False +def _addnoise(x, alpha, mu, sigma): + """Add Poisson-Gaussian noise to the data x + + Parameters + ---------- + x : float + The original data + + alpha : float + Level of noise gain + + mu : float + Level of noise offset + + sigma : float + Level of Gaussian noise + + Returns + ------- + y : float + The corrupted data + + """ + y = alpha * np.random.poisson(x / alpha) + mu + sigma * np.random.randn() + return y + +def PoissonGaussianNoiseGenerator(X, alpha=0.1, mu=0.1, sigma=0.1): + """Add Poisson-Gaussian noise to the data X + + Parameters + ---------- + X : array + The data to be corrupted + + alpha : float + Level of noise gain + + mu : float + Level of noise offset + + sigma : float + Level of Gaussian noise + + Returns + ------- + Y : array + The corrupted data + + """ + # Do some error checking + if alpha < 0. or alpha > 1.: + raise ValueError("alpha should be in range [0,1]") + # Vectorize noise function + addnoise = np.vectorize(_addnoise, otypes=[np.float]) + + # Rescale to [0,1] range + Xmax = np.amax(X) + X = X / Xmax + + # Add noise + Y = addnoise(X, alpha, mu, sigma) + + # Rescale to [0,1] range + Y = Y + np.abs(np.amin(Y)) + + # Rescale to X range + Y = Xmax * Y / np.amax(Y) + return Y class SVT(object): - """ + """ Parameters ---------- patchsize : integer - The dimensions of the patch in pixels + The dimensions of the patch in pixels to form a Casorati matrix (default = 4) length : integer Length in frames of the block to form a Casorati matrix. Must be odd (default = 15) - + optimize : bool Whether to optimize PGURE or just denoise according to given threshold (default = True) - + threshold : float Threshold to use if not optimizing PGURE (default = 0.5) - + alpha : float Level of noise gain, if negative then estimated online (default = -1) - + mu : float Level of noise offset, if negative then estimated online (default = -1) - + sigma : float Level of Gaussian noise, if negative then estimated online (default = -1) - + arpssize : integer Size of neighbourhood for ARPS search Must be odd (default = 7 pixels) - + tol : float Tolerance of PGURE optimizers (default = 1E-7) - + median : integer Size of initial median filter Must be odd (default = 5 pixels) - + """ - - def __init__(self, + def __init__(self, patchsize=4, patchoverlap=1, length=15, optimize=True, threshold=0.5, estimatenoise=True, - alpha=-1., - mu=-1., - sigma=-1., - arpssize=7, + alpha=-1., + mu=-1., + sigma=-1., + arpssize=7, tol=1e-7, median=5, - hotpixelthreshold=10 + hotpixelthreshold=10, + numthreads=1 ): - - # Load up parameters + + # Load up parameters self.patchsize = patchsize self.overlap = patchoverlap self.length = length self.optimize = optimize self.threshold = threshold - self.estimation = estimatenoise + self.estimation = estimatenoise self.alpha = alpha self.mu = mu self.sigma = sigma self.arpssize = arpssize - self.tol = tol + self.tol = tol self.median = median self.hotpixelthreshold = hotpixelthreshold - + self.numthreads = numthreads + # Do some error checking if self.overlap > self.patchsize: raise ValueError("Patch overlap should not be greater than patch size") @@ -135,71 +196,77 @@ class SVT(object): raise ValueError("Threshold should be in range [0,1]") if self.median % 2 == 0: raise ValueError("Median filter size should be odd") - + + # Check number of available CPU cores + num_cpu_cores = multiprocessing.cpu_count() + if self.numthreads > num_cpu_cores: + raise ValueError("Number of threads should be less than or equal to %d" % num_cpu_cores) + # Setup ctypes function self._PGURESVT = ctypes.cdll.LoadLibrary( - '${PYTHONLIBRARYPATH}/libpguresvt.so').PGURESVT - self._PGURESVT.restype = ctypes.c_int + '${PYTHONLIBRARYPATH}/libpguresvt.so').PGURESVT + self._PGURESVT.restype = ctypes.c_int self._PGURESVT.argtypes = [ndpointer(ctypes.c_double, flags="F"), ndpointer(ctypes.c_double, flags="F"), ndpointer(ctypes.c_int), - ctypes.c_int, + ctypes.c_int, ctypes.c_int, ctypes.c_int, ctypes.c_bool, ctypes.c_double, - ctypes.c_double, - ctypes.c_double, + ctypes.c_double, + ctypes.c_double, ctypes.c_double, ctypes.c_int, ctypes.c_double, ctypes.c_int, - ctypes.c_double] + ctypes.c_double, + ctypes.c_int] self.Y = None - + def denoise(self, X): """Denoise the data X - + Parameters ---------- X : array [nx, ny, time] The image sequence to be denoised - + Returns ------- self : object Returns the instance itself - - """ - + + """ + self._denoise(X) return self - + def _denoise(self, X): """Denoise the data X - + Parameters ---------- X : array [nx, ny, time] The image sequence to be denoised - + Returns ------- Y : array [nx, ny, time] Returns the denoised sequence - + """ # Check X is Fortran-order - X = self._check_array(X) + X = self._check_array(X) # Check sequence dimensions dims = np.asarray(X.shape).astype(np.int32) if dims[0] != dims[1] and self.estimation: raise ValueError("Quadtree noise estimation requires square images") if is_power_of_two(dims[0]) is False and self.estimation: raise ValueError("Quadtree noise estimation requires image dimensions 2^N") - # Create Y in memory - Y = np.zeros(X.shape, dtype=np.double, order='F') + # Create Y in memory + Y = np.zeros(X.shape, dtype=np.double, order='F') result = self._PGURESVT(X, Y, dims, @@ -214,25 +281,33 @@ class SVT(object): self.arpssize, self.tol, self.median, - self.hotpixelthreshold) + self.hotpixelthreshold, + self.numthreads) self.Y = Y return Y def _check_array(self, X): """Sanity-checks the data and parameters. - + Parameters ---------- X : array [nx, ny, time] The data as an array - + Returns ------- x : array [nx, ny, time] Returns the array in Fortran-order (column-major) - - """ + + """ x = np.copy(X.astype(np.double), order='F') return x - + def _is_power_of_two(self, n): + n = n/2 + if n == 2: + return True + elif n > 2: + return self._is_power_of_two(n) + else: + return False diff --git a/src/lib-PGURE-SVT.cpp b/src/lib-PGURE-SVT.cpp index 84e8fe2..f80930b 100644 --- a/src/lib-PGURE-SVT.cpp +++ b/src/lib-PGURE-SVT.cpp @@ -76,8 +76,8 @@ bool strToBool(std::string const& s) {return s != "0";}; // Main program extern "C" int PGURESVT(double *X, - double *Y, - int *dims, + double *Y, + int *dims, int Bs, int Bo, int T, @@ -89,7 +89,8 @@ extern "C" int PGURESVT(double *X, int MotionP, double tol, int MedianSize, - double hotpixelthreshold) { + double hotpixelthreshold, + int numthreads) { // Overall program timer auto overallstart = std::chrono::steady_clock::now(); @@ -101,52 +102,50 @@ extern "C" int PGURESVT(double *X, std::cout<<"Email: tjof2@cam.ac.uk"<= 0.) ? userLambda : 0.; - - int Nx = dims[0]; - int Ny = dims[1]; - int num_images = dims[2]; - - // Copy the image sequence into the a cube - arma::cube noisysequence(X, Nx, Ny, num_images); - - // Generate the clean and filtered sequences - arma::cube cleansequence(Nx, Ny, num_images); - arma::cube filteredsequence(Nx, Ny, num_images); - - cleansequence.zeros(); - filteredsequence.zeros(); - - // Parameters for median filter - unsigned short *Buffer = new unsigned short[Nx*Ny]; - unsigned short *FilteredBuffer = new unsigned short[Nx*Ny]; - int memsize = 512 * 1024; // L2 cache size - int filtsize = MedianSize; // Median filter size in pixels - - // Perform the initial median filtering - for (int i = 0; i < num_images; i++) { - arma::Mat curslice = arma::conv_to>::from(noisysequence.slice(i).eval()); - inplace_trans(curslice); - Buffer = curslice.memptr(); - ConstantTimeMedianFilter(Buffer, - FilteredBuffer, - Nx, Ny, Nx, Ny, - filtsize, 1, memsize); + // Set up OMP + omp_set_dynamic(0); + omp_set_num_threads(numthreads); + + int NoiseMethod = 4; + double lambda = (userLambda >= 0.) ? userLambda : 0.; + + int Nx = dims[0]; + int Ny = dims[1]; + int num_images = dims[2]; + + // Copy the image sequence into the a cube + arma::cube noisysequence(X, Nx, Ny, num_images); + + // Generate the clean and filtered sequences + arma::cube cleansequence(Nx, Ny, num_images); + arma::cube filteredsequence(Nx, Ny, num_images); + + cleansequence.zeros(); + filteredsequence.zeros(); + + // Parameters for median filter + unsigned short *Buffer = new unsigned short[Nx*Ny]; + unsigned short *FilteredBuffer = new unsigned short[Nx*Ny]; + int memsize = 512 * 1024; // L2 cache size + int filtsize = MedianSize; // Median filter size in pixels + + // Perform the initial median filtering + for (int i = 0; i < num_images; i++) { + arma::Mat curslice = arma::conv_to>::from(noisysequence.slice(i).eval()); + inplace_trans(curslice); + Buffer = curslice.memptr(); + ConstantTimeMedianFilter(Buffer, + FilteredBuffer, + Nx, Ny, Nx, Ny, + filtsize, 1, memsize); arma::Mat filslice(FilteredBuffer, Nx, Ny); inplace_trans(filslice); - filteredsequence.slice(i) = arma::conv_to::from(filslice); - } - - // Initial outlier detection (for hot pixels) - // using median absolute deviation - HotPixelFilter(noisysequence, hotpixelthreshold); - - ///////////////////////////// - // // - // START THE LOOP // - // // - ///////////////////////////// + filteredsequence.slice(i) = arma::conv_to::from(filslice); + } + + // Initial outlier detection (for hot pixels) + // using median absolute deviation + HotPixelFilter(noisysequence, hotpixelthreshold); // Print table headings int ww = 10; @@ -176,36 +175,25 @@ extern "C" int PGURESVT(double *X, u = noisysequence.slices(timeiter - framewindow, timeiter + framewindow); ufilter = filteredsequence.slices(timeiter - framewindow, timeiter + framewindow); } - + // Basic sequence normalization double inputmax = u.max(); u /= inputmax; ufilter /= ufilter.max(); - ///////////////////////////// - // // - // NOISE ESTIMATION // - // // - ///////////////////////////// - // Perform noise estimation if(pgureOpt) { - NoiseEstimator *noise = new NoiseEstimator; - noise->Estimate(u, - alpha, - mu, - sigma, - 8, - NoiseMethod); - delete noise; + NoiseEstimator *noise = new NoiseEstimator; + noise->Estimate(u, + alpha, + mu, + sigma, + 4, + NoiseMethod); + delete noise; } - - ///////////////////////////// - // // - // MOTION ESTIMATION // - // // - ///////////////////////////// + // Perform motion estimation MotionEstimator *motion = new MotionEstimator; motion->Estimate(ufilter, timeiter, @@ -216,12 +204,7 @@ extern "C" int PGURESVT(double *X, arma::icube sequencePatches = motion->GetEstimate(); delete motion; - ///////////////////////////// - // // - // PGURE OPTIMIZATION // - // // - ///////////////////////////// - + // Perform PGURE optimization PGURE *optimizer = new PGURE; optimizer->Initialize(u, sequencePatches, @@ -241,12 +224,6 @@ extern "C" int PGURESVT(double *X, } delete optimizer; - ///////////////////////////// - // // - // SEQUENCE RECONSTRUCTION // - // // - ///////////////////////////// - // Rescale back to original range v *= inputmax; @@ -269,29 +246,19 @@ extern "C" int PGURESVT(double *X, // Output a table row with noise estimates, lambda and timing std::ostringstream framestring; framestring << timeiter+1; - std::cout<(overallend - overallstart); std::cout<<"Total time: "<