From 5b4a3f1b2299b69010e5c8736847cfc2ff577d63 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 15:21:19 +0200
Subject: [PATCH 01/15] add blurred model image conv1 to observer

add blurred model image conv1 to observer

add conv1 to observer for use in algorithm convergence analysis
add some descriptions in the body of the algorithm
disambiguate data and inputData objects.
---
 .../simpleFlowDecAlgoConvergenceTest.py       | 139 ++++++++++++++++++
 python/flowdec/restoration.py                 |  11 +-
 2 files changed, 147 insertions(+), 3 deletions(-)
 create mode 100644 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
new file mode 100644
index 0000000..6d71be5
--- /dev/null
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -0,0 +1,139 @@
+# Dan White March 2019-2020
+# simpleFlowDecTest.py 
+# because the example python script isnt working on my setup
+# try this as a simpler more direct test - which seems to work.
+
+# flowdec pip package install/update requirements for anaconda on win10 64 bit:
+
+# install the pip package for flowdec with GPU support as per the flowdec install instructions:
+# github.com/hammerlab/flowdec/blob/master/README.md   but......
+
+
+# 22 april 2020 - what works today :
+# flowdec uses tensorflow which on nvidda GPU uses CUDA so need to install the stuff here
+# https://www.tensorflow.org/install/gpu   but.....those instructionsd seem to install incompatible stuff... so try
+# things were installed in this order and the script works and used the GPU
+#cuda toolkit 10.0 
+#cudnn-10.0 7.6.34.38 installed into C:/tools/
+#pip install flowdec
+#     not pip install flowdec[tf_gpu  
+# ommitting the tf_gpu option, so it leaves tensorflow-gpu uninstalled, becasue by default
+# by now it installs v2.1 of tensorflow which doesnt seem to work for flowdec? 
+# pip install tensorflow-gpu==1.14.0 (2.0 might work...? maybe needs higher cuda version)
+# Need windows env variables pointing to cuda stuff:  CUDA and CUPTI and another related library cuDNN
+# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\bin;%PATH%
+# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\extras\CUPTI\libx64;%PATH%
+# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\include;%PATH%
+# SET PATH=C:\tools\cuda\bin;%PATH%
+# some of these don't seem to be sticky???
+
+# for time benchmarking
+import time
+# for logging to text file
+import sys
+import logging
+
+#send std output to a log file
+#sys.stdout = open('FlowDecLog.txt', 'w')
+logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
+
+startImports = time.process_time()   
+from skimage.external.tifffile import imsave, imread
+from skimage.metrics import structural_similarity as ssim
+from flowdec import data as fd_data
+from flowdec import restoration as fd_restoration
+importsTime = (time.process_time() - startImports)
+
+# this seems to return 0 seconds, so maybe TF math is already imported by flowde? 
+startTFMathImport = time.process_time()
+from tensorflow import math as tfmath
+importTFMathTime = (time.process_time() - startTFMathImport)
+
+# Load test image from same dir as we execute in
+rawBig = 'C1-YeastTNA1_1516_conv_RG_26oC_003sub100.tif'
+rawSmall = 'C1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tif'
+rawImg = rawSmall
+
+raw = imread(rawImg)
+
+# Load psf kernel image from same dir
+# A cropped 64x64 PSF to reduce memory use, 21 z slices, 0.125 nm spacing in z
+PSFsmall = 'gpsf_3D_1514_a3_001_WF-sub105crop64.tif'
+# The same PSF but not cropped so much, 128x128, probably uses much more memory.  
+PSFbigger = 'gpsf_3D_1514_a3_001_WF-sub105crop128.tif'
+# choose the psf to use from the options above. 
+PSF = PSFsmall
+print (PSF)
+kernel = imread(PSF)
+
+#base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
+base_iter = 100
+
+'''
+# Create an observer function to monitor converegence, 
+# where the first argument is the current state of the
+# image as it is being deconvolved, i is the current iteration number (1-based),
+# third is the padded current guess image, and fourth is the value of covergence metric R,
+# and the remaining arguments (at TOW) only include the uncropped image result which
+# is generally not useful for anything other than development or debugging
+# imgs = []
+def observer(decon_crop, i, decon_pad, conv1, *args):
+    #imgs.append(decon_crop)
+    if i % 1 == 0:
+        sumRaw = raw.sum()
+        sumBlurredModel = decon_crop.sum()
+	    #compute convergence residuals between raw image and blurred model image
+		# with currectn sted data these get bigger not smaller with iterations... weird..
+        convergenceResiduals = abs(sumRaw - sumBlurredModel)
+        convergenceR = convergenceResiduals / sumRaw
+		# lets try structural similarity (as it's supposed to be better then mean square error)
+		# which indeed seems to track convergence of this dataset in a way that seems to match the image results. 
+        structSim = ssim(raw, decon_crop, data_range=decon_crop.max() - decon_crop.min())
+        print('Iter={}, RawSum={:.3f}, DeconSum={:.3f}, DeconMax={:.3f}, DeconStDev={:.3f}, SumResiduals={:.3f}, ConvergeR={:.16f}), SSIM= {:.3f})'.format(
+		   i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
+'''
+# Run the deconvolution process and note that deconvolution initialization is best kept separate from 
+# execution since the "initialize" operation corresponds to creating a TensorFlow graph, which is a 
+# relatively expensive operation and should not be repeated across multiple executions
+
+# initialize the TF graph for the deconvolution settings in use for certain sized input and psf images
+# works for doing the same input data multiple times with different iteractions
+# should work for doing different input data with same sizes of image and psf, 
+# eg a time series split into tiff 1 file per time point???? 
+startAlgoinit = time.process_time()   
+# with observer function to track concvergence
+#algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim, observer_fn=observer).initialize()
+#without observer function - much faster obvs. 
+algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim).initialize()
+TFinitTime = (time.process_time() - startAlgoinit)
+
+# run the deconvolution itself
+# in a loop making different numbers of iterations, multiples of base value of n_iter
+multiRunFactor = 1
+timingListIter = []
+timingListTime = []
+for i in range(1, multiRunFactor+1):
+  niter = (base_iter*i)
+  # start measuring time
+  startDec = time.process_time()
+  res = algo.run(fd_data.Acquisition(data=raw, kernel=kernel), niter=(niter)).data
+  # measure time here includes only the deconvolution, no file saving
+  DecTime = (time.process_time() - startDec)
+  # save the result # using skimage.external.tifffile.imsave
+  resultFileName = ('result' + rawImg + PSF + str(niter) + 'iterations.tif')
+  imsave(resultFileName, res)
+  # measure time here includes file saving
+  #DecTime = (time.process_time() - startDec)
+  print('Saved result image TIFF file ' + resultFileName)
+  print(str(DecTime) + ' is how many sec ' + str(niter) + ' iterations took.')
+  timingListIter.append(niter)
+  timingListTime.append(DecTime)
+
+#benchmarking data output
+print (str(importsTime) + ' seconds to import flowdec, TF and CUDA libs etc.')
+print (str(importTFMathTime) + ' seconds to import TF Math')
+print (str(TFinitTime) + ' sec is the tensorFlow initialisation time')
+print ('Pairs of values of iterations done vs time in seconds')
+print (timingListIter)
+print (timingListTime)
+print('Done')
\ No newline at end of file
diff --git a/python/flowdec/restoration.py b/python/flowdec/restoration.py
index dc074bb..352a221 100755
--- a/python/flowdec/restoration.py
+++ b/python/flowdec/restoration.py
@@ -289,17 +289,19 @@ def _build_tf_graph(self):
         def cond(i, decon):
             return i <= niter
 
-        def conv(data, kernel_fft):
-            return tf.math.real(fft_rev(fft_fwd(tf.cast(data, self.fft_dtype)) * kernel_fft))
+        def conv(inputData, kernel_fft):
+            return tf.math.real(fft_rev(fft_fwd(tf.cast(inputData, self.fft_dtype)) * kernel_fft))
 
         def body(i, decon):
             # Richardson-Lucy Iteration - logic taken largely from a combination of
             # the scikit-image (real domain) and DeconvolutionLab2 implementations (complex domain)
+            # conv1 is the current model blurred with the PSF
             conv1 = conv(decon, kern_fft)
 
             # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)
             blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
 
+            # conv2 is the blurred model convolved with the flipped PSF
             conv2 = conv(blur1, kern_fft_conj)
 
             # Positivity constraint on result for iteration
@@ -309,7 +311,10 @@ def body(i, decon):
             if self.observer_fn is not None:
                 # Remove any cropping that may have been added as this is usually not desirable in observers
                 decon_crop = unpad_around_center(decon, tf.shape(datah))
-                _, i, decon = tf_observer([decon_crop, i, decon], self.observer_fn)
+                # we can use these capurured observed tensors to evaluate eg convergence
+                # in eg. the observer function used.
+                _, i, decon, conv1  = tf_observer(
+				    [decon_crop, i, decon, conv1], self.observer_fn)
 
             return i + 1, decon
 

From d571d6eb718bed2aa22850dc56f15c4f54637ba0 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 15:22:40 +0200
Subject: [PATCH 02/15] Revert "add blurred model image conv1 to observer"

This reverts commit 5b4a3f1b2299b69010e5c8736847cfc2ff577d63.
---
 .../simpleFlowDecAlgoConvergenceTest.py       | 139 ------------------
 python/flowdec/restoration.py                 |  11 +-
 2 files changed, 3 insertions(+), 147 deletions(-)
 delete mode 100644 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
deleted file mode 100644
index 6d71be5..0000000
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ /dev/null
@@ -1,139 +0,0 @@
-# Dan White March 2019-2020
-# simpleFlowDecTest.py 
-# because the example python script isnt working on my setup
-# try this as a simpler more direct test - which seems to work.
-
-# flowdec pip package install/update requirements for anaconda on win10 64 bit:
-
-# install the pip package for flowdec with GPU support as per the flowdec install instructions:
-# github.com/hammerlab/flowdec/blob/master/README.md   but......
-
-
-# 22 april 2020 - what works today :
-# flowdec uses tensorflow which on nvidda GPU uses CUDA so need to install the stuff here
-# https://www.tensorflow.org/install/gpu   but.....those instructionsd seem to install incompatible stuff... so try
-# things were installed in this order and the script works and used the GPU
-#cuda toolkit 10.0 
-#cudnn-10.0 7.6.34.38 installed into C:/tools/
-#pip install flowdec
-#     not pip install flowdec[tf_gpu  
-# ommitting the tf_gpu option, so it leaves tensorflow-gpu uninstalled, becasue by default
-# by now it installs v2.1 of tensorflow which doesnt seem to work for flowdec? 
-# pip install tensorflow-gpu==1.14.0 (2.0 might work...? maybe needs higher cuda version)
-# Need windows env variables pointing to cuda stuff:  CUDA and CUPTI and another related library cuDNN
-# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\bin;%PATH%
-# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\extras\CUPTI\libx64;%PATH%
-# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\include;%PATH%
-# SET PATH=C:\tools\cuda\bin;%PATH%
-# some of these don't seem to be sticky???
-
-# for time benchmarking
-import time
-# for logging to text file
-import sys
-import logging
-
-#send std output to a log file
-#sys.stdout = open('FlowDecLog.txt', 'w')
-logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
-
-startImports = time.process_time()   
-from skimage.external.tifffile import imsave, imread
-from skimage.metrics import structural_similarity as ssim
-from flowdec import data as fd_data
-from flowdec import restoration as fd_restoration
-importsTime = (time.process_time() - startImports)
-
-# this seems to return 0 seconds, so maybe TF math is already imported by flowde? 
-startTFMathImport = time.process_time()
-from tensorflow import math as tfmath
-importTFMathTime = (time.process_time() - startTFMathImport)
-
-# Load test image from same dir as we execute in
-rawBig = 'C1-YeastTNA1_1516_conv_RG_26oC_003sub100.tif'
-rawSmall = 'C1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tif'
-rawImg = rawSmall
-
-raw = imread(rawImg)
-
-# Load psf kernel image from same dir
-# A cropped 64x64 PSF to reduce memory use, 21 z slices, 0.125 nm spacing in z
-PSFsmall = 'gpsf_3D_1514_a3_001_WF-sub105crop64.tif'
-# The same PSF but not cropped so much, 128x128, probably uses much more memory.  
-PSFbigger = 'gpsf_3D_1514_a3_001_WF-sub105crop128.tif'
-# choose the psf to use from the options above. 
-PSF = PSFsmall
-print (PSF)
-kernel = imread(PSF)
-
-#base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
-base_iter = 100
-
-'''
-# Create an observer function to monitor converegence, 
-# where the first argument is the current state of the
-# image as it is being deconvolved, i is the current iteration number (1-based),
-# third is the padded current guess image, and fourth is the value of covergence metric R,
-# and the remaining arguments (at TOW) only include the uncropped image result which
-# is generally not useful for anything other than development or debugging
-# imgs = []
-def observer(decon_crop, i, decon_pad, conv1, *args):
-    #imgs.append(decon_crop)
-    if i % 1 == 0:
-        sumRaw = raw.sum()
-        sumBlurredModel = decon_crop.sum()
-	    #compute convergence residuals between raw image and blurred model image
-		# with currectn sted data these get bigger not smaller with iterations... weird..
-        convergenceResiduals = abs(sumRaw - sumBlurredModel)
-        convergenceR = convergenceResiduals / sumRaw
-		# lets try structural similarity (as it's supposed to be better then mean square error)
-		# which indeed seems to track convergence of this dataset in a way that seems to match the image results. 
-        structSim = ssim(raw, decon_crop, data_range=decon_crop.max() - decon_crop.min())
-        print('Iter={}, RawSum={:.3f}, DeconSum={:.3f}, DeconMax={:.3f}, DeconStDev={:.3f}, SumResiduals={:.3f}, ConvergeR={:.16f}), SSIM= {:.3f})'.format(
-		   i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
-'''
-# Run the deconvolution process and note that deconvolution initialization is best kept separate from 
-# execution since the "initialize" operation corresponds to creating a TensorFlow graph, which is a 
-# relatively expensive operation and should not be repeated across multiple executions
-
-# initialize the TF graph for the deconvolution settings in use for certain sized input and psf images
-# works for doing the same input data multiple times with different iteractions
-# should work for doing different input data with same sizes of image and psf, 
-# eg a time series split into tiff 1 file per time point???? 
-startAlgoinit = time.process_time()   
-# with observer function to track concvergence
-#algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim, observer_fn=observer).initialize()
-#without observer function - much faster obvs. 
-algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim).initialize()
-TFinitTime = (time.process_time() - startAlgoinit)
-
-# run the deconvolution itself
-# in a loop making different numbers of iterations, multiples of base value of n_iter
-multiRunFactor = 1
-timingListIter = []
-timingListTime = []
-for i in range(1, multiRunFactor+1):
-  niter = (base_iter*i)
-  # start measuring time
-  startDec = time.process_time()
-  res = algo.run(fd_data.Acquisition(data=raw, kernel=kernel), niter=(niter)).data
-  # measure time here includes only the deconvolution, no file saving
-  DecTime = (time.process_time() - startDec)
-  # save the result # using skimage.external.tifffile.imsave
-  resultFileName = ('result' + rawImg + PSF + str(niter) + 'iterations.tif')
-  imsave(resultFileName, res)
-  # measure time here includes file saving
-  #DecTime = (time.process_time() - startDec)
-  print('Saved result image TIFF file ' + resultFileName)
-  print(str(DecTime) + ' is how many sec ' + str(niter) + ' iterations took.')
-  timingListIter.append(niter)
-  timingListTime.append(DecTime)
-
-#benchmarking data output
-print (str(importsTime) + ' seconds to import flowdec, TF and CUDA libs etc.')
-print (str(importTFMathTime) + ' seconds to import TF Math')
-print (str(TFinitTime) + ' sec is the tensorFlow initialisation time')
-print ('Pairs of values of iterations done vs time in seconds')
-print (timingListIter)
-print (timingListTime)
-print('Done')
\ No newline at end of file
diff --git a/python/flowdec/restoration.py b/python/flowdec/restoration.py
index 352a221..dc074bb 100755
--- a/python/flowdec/restoration.py
+++ b/python/flowdec/restoration.py
@@ -289,19 +289,17 @@ def _build_tf_graph(self):
         def cond(i, decon):
             return i <= niter
 
-        def conv(inputData, kernel_fft):
-            return tf.math.real(fft_rev(fft_fwd(tf.cast(inputData, self.fft_dtype)) * kernel_fft))
+        def conv(data, kernel_fft):
+            return tf.math.real(fft_rev(fft_fwd(tf.cast(data, self.fft_dtype)) * kernel_fft))
 
         def body(i, decon):
             # Richardson-Lucy Iteration - logic taken largely from a combination of
             # the scikit-image (real domain) and DeconvolutionLab2 implementations (complex domain)
-            # conv1 is the current model blurred with the PSF
             conv1 = conv(decon, kern_fft)
 
             # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)
             blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
 
-            # conv2 is the blurred model convolved with the flipped PSF
             conv2 = conv(blur1, kern_fft_conj)
 
             # Positivity constraint on result for iteration
@@ -311,10 +309,7 @@ def body(i, decon):
             if self.observer_fn is not None:
                 # Remove any cropping that may have been added as this is usually not desirable in observers
                 decon_crop = unpad_around_center(decon, tf.shape(datah))
-                # we can use these capurured observed tensors to evaluate eg convergence
-                # in eg. the observer function used.
-                _, i, decon, conv1  = tf_observer(
-				    [decon_crop, i, decon, conv1], self.observer_fn)
+                _, i, decon = tf_observer([decon_crop, i, decon], self.observer_fn)
 
             return i + 1, decon
 

From 550a08e7e438f2f382e3f546ece6c2981f43021d Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 15:25:51 +0200
Subject: [PATCH 03/15] Revert "Revert "add blurred model image conv1 to
 observer""

This reverts commit d571d6eb718bed2aa22850dc56f15c4f54637ba0.
---
 .../simpleFlowDecAlgoConvergenceTest.py       | 139 ++++++++++++++++++
 python/flowdec/restoration.py                 |  11 +-
 2 files changed, 147 insertions(+), 3 deletions(-)
 create mode 100644 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
new file mode 100644
index 0000000..6d71be5
--- /dev/null
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -0,0 +1,139 @@
+# Dan White March 2019-2020
+# simpleFlowDecTest.py 
+# because the example python script isnt working on my setup
+# try this as a simpler more direct test - which seems to work.
+
+# flowdec pip package install/update requirements for anaconda on win10 64 bit:
+
+# install the pip package for flowdec with GPU support as per the flowdec install instructions:
+# github.com/hammerlab/flowdec/blob/master/README.md   but......
+
+
+# 22 april 2020 - what works today :
+# flowdec uses tensorflow which on nvidda GPU uses CUDA so need to install the stuff here
+# https://www.tensorflow.org/install/gpu   but.....those instructionsd seem to install incompatible stuff... so try
+# things were installed in this order and the script works and used the GPU
+#cuda toolkit 10.0 
+#cudnn-10.0 7.6.34.38 installed into C:/tools/
+#pip install flowdec
+#     not pip install flowdec[tf_gpu  
+# ommitting the tf_gpu option, so it leaves tensorflow-gpu uninstalled, becasue by default
+# by now it installs v2.1 of tensorflow which doesnt seem to work for flowdec? 
+# pip install tensorflow-gpu==1.14.0 (2.0 might work...? maybe needs higher cuda version)
+# Need windows env variables pointing to cuda stuff:  CUDA and CUPTI and another related library cuDNN
+# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\bin;%PATH%
+# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\extras\CUPTI\libx64;%PATH%
+# SET PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\include;%PATH%
+# SET PATH=C:\tools\cuda\bin;%PATH%
+# some of these don't seem to be sticky???
+
+# for time benchmarking
+import time
+# for logging to text file
+import sys
+import logging
+
+#send std output to a log file
+#sys.stdout = open('FlowDecLog.txt', 'w')
+logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
+
+startImports = time.process_time()   
+from skimage.external.tifffile import imsave, imread
+from skimage.metrics import structural_similarity as ssim
+from flowdec import data as fd_data
+from flowdec import restoration as fd_restoration
+importsTime = (time.process_time() - startImports)
+
+# this seems to return 0 seconds, so maybe TF math is already imported by flowde? 
+startTFMathImport = time.process_time()
+from tensorflow import math as tfmath
+importTFMathTime = (time.process_time() - startTFMathImport)
+
+# Load test image from same dir as we execute in
+rawBig = 'C1-YeastTNA1_1516_conv_RG_26oC_003sub100.tif'
+rawSmall = 'C1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tif'
+rawImg = rawSmall
+
+raw = imread(rawImg)
+
+# Load psf kernel image from same dir
+# A cropped 64x64 PSF to reduce memory use, 21 z slices, 0.125 nm spacing in z
+PSFsmall = 'gpsf_3D_1514_a3_001_WF-sub105crop64.tif'
+# The same PSF but not cropped so much, 128x128, probably uses much more memory.  
+PSFbigger = 'gpsf_3D_1514_a3_001_WF-sub105crop128.tif'
+# choose the psf to use from the options above. 
+PSF = PSFsmall
+print (PSF)
+kernel = imread(PSF)
+
+#base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
+base_iter = 100
+
+'''
+# Create an observer function to monitor converegence, 
+# where the first argument is the current state of the
+# image as it is being deconvolved, i is the current iteration number (1-based),
+# third is the padded current guess image, and fourth is the value of covergence metric R,
+# and the remaining arguments (at TOW) only include the uncropped image result which
+# is generally not useful for anything other than development or debugging
+# imgs = []
+def observer(decon_crop, i, decon_pad, conv1, *args):
+    #imgs.append(decon_crop)
+    if i % 1 == 0:
+        sumRaw = raw.sum()
+        sumBlurredModel = decon_crop.sum()
+	    #compute convergence residuals between raw image and blurred model image
+		# with currectn sted data these get bigger not smaller with iterations... weird..
+        convergenceResiduals = abs(sumRaw - sumBlurredModel)
+        convergenceR = convergenceResiduals / sumRaw
+		# lets try structural similarity (as it's supposed to be better then mean square error)
+		# which indeed seems to track convergence of this dataset in a way that seems to match the image results. 
+        structSim = ssim(raw, decon_crop, data_range=decon_crop.max() - decon_crop.min())
+        print('Iter={}, RawSum={:.3f}, DeconSum={:.3f}, DeconMax={:.3f}, DeconStDev={:.3f}, SumResiduals={:.3f}, ConvergeR={:.16f}), SSIM= {:.3f})'.format(
+		   i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
+'''
+# Run the deconvolution process and note that deconvolution initialization is best kept separate from 
+# execution since the "initialize" operation corresponds to creating a TensorFlow graph, which is a 
+# relatively expensive operation and should not be repeated across multiple executions
+
+# initialize the TF graph for the deconvolution settings in use for certain sized input and psf images
+# works for doing the same input data multiple times with different iteractions
+# should work for doing different input data with same sizes of image and psf, 
+# eg a time series split into tiff 1 file per time point???? 
+startAlgoinit = time.process_time()   
+# with observer function to track concvergence
+#algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim, observer_fn=observer).initialize()
+#without observer function - much faster obvs. 
+algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim).initialize()
+TFinitTime = (time.process_time() - startAlgoinit)
+
+# run the deconvolution itself
+# in a loop making different numbers of iterations, multiples of base value of n_iter
+multiRunFactor = 1
+timingListIter = []
+timingListTime = []
+for i in range(1, multiRunFactor+1):
+  niter = (base_iter*i)
+  # start measuring time
+  startDec = time.process_time()
+  res = algo.run(fd_data.Acquisition(data=raw, kernel=kernel), niter=(niter)).data
+  # measure time here includes only the deconvolution, no file saving
+  DecTime = (time.process_time() - startDec)
+  # save the result # using skimage.external.tifffile.imsave
+  resultFileName = ('result' + rawImg + PSF + str(niter) + 'iterations.tif')
+  imsave(resultFileName, res)
+  # measure time here includes file saving
+  #DecTime = (time.process_time() - startDec)
+  print('Saved result image TIFF file ' + resultFileName)
+  print(str(DecTime) + ' is how many sec ' + str(niter) + ' iterations took.')
+  timingListIter.append(niter)
+  timingListTime.append(DecTime)
+
+#benchmarking data output
+print (str(importsTime) + ' seconds to import flowdec, TF and CUDA libs etc.')
+print (str(importTFMathTime) + ' seconds to import TF Math')
+print (str(TFinitTime) + ' sec is the tensorFlow initialisation time')
+print ('Pairs of values of iterations done vs time in seconds')
+print (timingListIter)
+print (timingListTime)
+print('Done')
\ No newline at end of file
diff --git a/python/flowdec/restoration.py b/python/flowdec/restoration.py
index dc074bb..352a221 100755
--- a/python/flowdec/restoration.py
+++ b/python/flowdec/restoration.py
@@ -289,17 +289,19 @@ def _build_tf_graph(self):
         def cond(i, decon):
             return i <= niter
 
-        def conv(data, kernel_fft):
-            return tf.math.real(fft_rev(fft_fwd(tf.cast(data, self.fft_dtype)) * kernel_fft))
+        def conv(inputData, kernel_fft):
+            return tf.math.real(fft_rev(fft_fwd(tf.cast(inputData, self.fft_dtype)) * kernel_fft))
 
         def body(i, decon):
             # Richardson-Lucy Iteration - logic taken largely from a combination of
             # the scikit-image (real domain) and DeconvolutionLab2 implementations (complex domain)
+            # conv1 is the current model blurred with the PSF
             conv1 = conv(decon, kern_fft)
 
             # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)
             blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
 
+            # conv2 is the blurred model convolved with the flipped PSF
             conv2 = conv(blur1, kern_fft_conj)
 
             # Positivity constraint on result for iteration
@@ -309,7 +311,10 @@ def body(i, decon):
             if self.observer_fn is not None:
                 # Remove any cropping that may have been added as this is usually not desirable in observers
                 decon_crop = unpad_around_center(decon, tf.shape(datah))
-                _, i, decon = tf_observer([decon_crop, i, decon], self.observer_fn)
+                # we can use these capurured observed tensors to evaluate eg convergence
+                # in eg. the observer function used.
+                _, i, decon, conv1  = tf_observer(
+				    [decon_crop, i, decon, conv1], self.observer_fn)
 
             return i + 1, decon
 

From f947936e240208289067b6a4a7598443f997201e Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 15:31:35 +0200
Subject: [PATCH 04/15] reset file name and description

---
 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 6d71be5..9968264 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -1,7 +1,6 @@
 # Dan White March 2019-2020
-# simpleFlowDecTest.py 
-# because the example python script isnt working on my setup
-# try this as a simpler more direct test - which seems to work.
+# simpleFlowDecAlgoConvergenceTest.py 
+# A simple script usiing flowdec, testing convergence of the algorithm on test data
 
 # flowdec pip package install/update requirements for anaconda on win10 64 bit:
 

From 090ab8ace14768721c70c357dceae6e1f40ded18 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 15:59:05 +0200
Subject: [PATCH 05/15] fix typos

---
 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 9968264..998f2d9 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -33,7 +33,7 @@
 import logging
 
 #send std output to a log file
-#sys.stdout = open('FlowDecLog.txt', 'w')
+sys.stdout = open('FlowDecLog.txt', 'w')
 logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
 
 startImports = time.process_time()   

From b9777e4a1e3b4b4788bc8f7829de962022ebb8f4 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 15:59:46 +0200
Subject: [PATCH 06/15] typo

---
 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 998f2d9..428076d 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -1,6 +1,6 @@
 # Dan White March 2019-2020
 # simpleFlowDecAlgoConvergenceTest.py 
-# A simple script usiing flowdec, testing convergence of the algorithm on test data
+# A simple script using flowdec, testing convergence of the algorithm on test data
 
 # flowdec pip package install/update requirements for anaconda on win10 64 bit:
 

From bda68c2c504a2c890de188ffcfba20182856343c Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 16:00:22 +0200
Subject: [PATCH 07/15] turn on stdout redirection to text file for covergence
 results collection

---
 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 428076d..7ac2c53 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -10,7 +10,7 @@
 
 # 22 april 2020 - what works today :
 # flowdec uses tensorflow which on nvidda GPU uses CUDA so need to install the stuff here
-# https://www.tensorflow.org/install/gpu   but.....those instructionsd seem to install incompatible stuff... so try
+# https://www.tensorflow.org/install/gpu   but.....those instructions seem to install incompatible stuff... so try
 # things were installed in this order and the script works and used the GPU
 #cuda toolkit 10.0 
 #cudnn-10.0 7.6.34.38 installed into C:/tools/

From 61f1b49ef7937ef49527764f922c6a85b91ea96b Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 16:36:02 +0200
Subject: [PATCH 08/15] turn off high level logging, fix typos

---
 .gitignore                                           |  7 +++++++
 .../scripts/simpleFlowDecAlgoConvergenceTest.py      | 12 ++++++------
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/.gitignore b/.gitignore
index ed98de9..20be2da 100644
--- a/.gitignore
+++ b/.gitignore
@@ -86,3 +86,10 @@ tmp/
 *~.nib
 local.properties
 .settings/
+python/examples/scripts/C1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tif
+python/examples/scripts/C1-YeastTNA1_1516_conv_RG_26oC_003sub100.tif
+python/examples/scripts/gpsf_3D_1514_a3_001_WF-sub105crop128.tif
+python/examples/scripts/gpsf_3D_1514_a3_001_WF-sub105crop64.tif
+python/examples/scripts/FlowDecLog.txt
+python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif100iterations.tif
+python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif50iterations.tif
diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 7ac2c53..189aed3 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -30,11 +30,11 @@
 import time
 # for logging to text file
 import sys
-import logging
+#import logging
 
 #send std output to a log file
 sys.stdout = open('FlowDecLog.txt', 'w')
-logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
+#logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
 
 startImports = time.process_time()   
 from skimage.external.tifffile import imsave, imread
@@ -43,7 +43,7 @@
 from flowdec import restoration as fd_restoration
 importsTime = (time.process_time() - startImports)
 
-# this seems to return 0 seconds, so maybe TF math is already imported by flowde? 
+# this seems to return 0 seconds, so maybe TF math is already imported by flowdec? 
 startTFMathImport = time.process_time()
 from tensorflow import math as tfmath
 importTFMathTime = (time.process_time() - startTFMathImport)
@@ -81,11 +81,11 @@ def observer(decon_crop, i, decon_pad, conv1, *args):
     if i % 1 == 0:
         sumRaw = raw.sum()
         sumBlurredModel = decon_crop.sum()
-	    #compute convergence residuals between raw image and blurred model image
-		# with currectn sted data these get bigger not smaller with iterations... weird..
+	    # compute convergence residuals between raw image and blurred model image
+		# with current test data these get bigger not smaller with iterations... weird..
         convergenceResiduals = abs(sumRaw - sumBlurredModel)
         convergenceR = convergenceResiduals / sumRaw
-		# lets try structural similarity (as it's supposed to be better then mean square error)
+		# let's try structural similarity (as it's supposed to be better then mean square error)
 		# which indeed seems to track convergence of this dataset in a way that seems to match the image results. 
         structSim = ssim(raw, decon_crop, data_range=decon_crop.max() - decon_crop.min())
         print('Iter={}, RawSum={:.3f}, DeconSum={:.3f}, DeconMax={:.3f}, DeconStDev={:.3f}, SumResiduals={:.3f}, ConvergeR={:.16f}), SSIM= {:.3f})'.format(

From ff32b76aeadc60685f82330a29aa86755d3d7245 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 16:36:58 +0200
Subject: [PATCH 09/15] do algo convergence stats every 5th iteration

---
 .../simpleFlowDecAlgoConvergenceTest.py       | 20 +++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 189aed3..1e2ab5b 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -66,10 +66,10 @@
 kernel = imread(PSF)
 
 #base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
-base_iter = 100
+base_iter = 50
 
-'''
-# Create an observer function to monitor converegence, 
+
+# Create an observer function to monitor convergence, 
 # where the first argument is the current state of the
 # image as it is being deconvolved, i is the current iteration number (1-based),
 # third is the padded current guess image, and fourth is the value of covergence metric R,
@@ -78,7 +78,7 @@
 # imgs = []
 def observer(decon_crop, i, decon_pad, conv1, *args):
     #imgs.append(decon_crop)
-    if i % 1 == 0:
+    if i % 5 == 0:
         sumRaw = raw.sum()
         sumBlurredModel = decon_crop.sum()
 	    # compute convergence residuals between raw image and blurred model image
@@ -88,9 +88,9 @@ def observer(decon_crop, i, decon_pad, conv1, *args):
 		# let's try structural similarity (as it's supposed to be better then mean square error)
 		# which indeed seems to track convergence of this dataset in a way that seems to match the image results. 
         structSim = ssim(raw, decon_crop, data_range=decon_crop.max() - decon_crop.min())
-        print('Iter={}, RawSum={:.3f}, DeconSum={:.3f}, DeconMax={:.3f}, DeconStDev={:.3f}, SumResiduals={:.3f}, ConvergeR={:.16f}), SSIM= {:.3f})'.format(
+        print('Iter,{},RawSum,{:.3f},DeconSum,{:.3f},DeconMax,{:.3f},DeconStDev,{:.3f},SumResiduals,{:.3f},ConvergeR,{:.16f},SSIM,{:.3f}'.format(
 		   i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
-'''
+
 # Run the deconvolution process and note that deconvolution initialization is best kept separate from 
 # execution since the "initialize" operation corresponds to creating a TensorFlow graph, which is a 
 # relatively expensive operation and should not be repeated across multiple executions
@@ -100,10 +100,10 @@ def observer(decon_crop, i, decon_pad, conv1, *args):
 # should work for doing different input data with same sizes of image and psf, 
 # eg a time series split into tiff 1 file per time point???? 
 startAlgoinit = time.process_time()   
-# with observer function to track concvergence
-#algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim, observer_fn=observer).initialize()
-#without observer function - much faster obvs. 
-algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim).initialize()
+# Run algorithm with observer function to track concvergence
+algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim, observer_fn=observer).initialize()
+# Run algorithm without observer function - much faster obvs. 
+#algo = fd_restoration.RichardsonLucyDeconvolver(raw.ndim).initialize()
 TFinitTime = (time.process_time() - startAlgoinit)
 
 # run the deconvolution itself

From 11e6734acb2d894f9593b2a653de5316e952fd02 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 28 Apr 2020 18:30:14 +0200
Subject: [PATCH 10/15] 1st try at Gold (Ratio) algorithm

needed to normalise tensors to prevent overflows when multipling etc.
needed at add delta parameter tonumerator and denominator to prevent div by close to zero
---
 .gitignore                                    |  4 +++
 .../simpleFlowDecAlgoConvergenceTest.py       |  6 ++--
 python/flowdec/restoration.py                 | 30 +++++++++++++++++--
 3 files changed, 34 insertions(+), 6 deletions(-)

diff --git a/.gitignore b/.gitignore
index 20be2da..7e4e5ab 100644
--- a/.gitignore
+++ b/.gitignore
@@ -93,3 +93,7 @@ python/examples/scripts/gpsf_3D_1514_a3_001_WF-sub105crop64.tif
 python/examples/scripts/FlowDecLog.txt
 python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif100iterations.tif
 python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif50iterations.tif
+python/examples/scripts/FlowDecLogGold50.txt
+python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif250iterations.tif
+python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif4iterations.tif
+python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif5iterations.tif
diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 1e2ab5b..7c871e4 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -33,7 +33,7 @@
 #import logging
 
 #send std output to a log file
-sys.stdout = open('FlowDecLog.txt', 'w')
+sys.stdout = open('FlowDecLogGold50.txt', 'w')
 #logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
 
 startImports = time.process_time()   
@@ -66,7 +66,7 @@
 kernel = imread(PSF)
 
 #base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
-base_iter = 50
+base_iter = 250
 
 
 # Create an observer function to monitor convergence, 
@@ -78,7 +78,7 @@
 # imgs = []
 def observer(decon_crop, i, decon_pad, conv1, *args):
     #imgs.append(decon_crop)
-    if i % 5 == 0:
+    if i % 10 == 0:
         sumRaw = raw.sum()
         sumBlurredModel = decon_crop.sum()
 	    # compute convergence residuals between raw image and blurred model image
diff --git a/python/flowdec/restoration.py b/python/flowdec/restoration.py
index 352a221..099e453 100755
--- a/python/flowdec/restoration.py
+++ b/python/flowdec/restoration.py
@@ -293,19 +293,43 @@ def conv(inputData, kernel_fft):
             return tf.math.real(fft_rev(fft_fwd(tf.cast(inputData, self.fft_dtype)) * kernel_fft))
 
         def body(i, decon):
-            # Richardson-Lucy Iteration - logic taken largely from a combination of
+            '''# Richardson-Lucy Iteration - logic taken largely from a combination of
             # the scikit-image (real domain) and DeconvolutionLab2 implementations (complex domain)
             # conv1 is the current model blurred with the PSF
             conv1 = conv(decon, kern_fft)
 
             # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)
-            blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
+			blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
 
             # conv2 is the blurred model convolved with the flipped PSF
             conv2 = conv(blur1, kern_fft_conj)
 
             # Positivity constraint on result for iteration
-            decon = tf.maximum(decon * conv2, 0.)
+			decon = tf.maximum(decon * conv2, 0.)
+            '''
+
+            # Gold algorithm, ratio method, simpler then RL, doesnt use flipped OTF
+            # conv1 is the current model blurred with the PSF
+            conv1 = conv(decon, kern_fft)
+
+            # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)
+            # we wont do it here as we will use the delta values in denom and numerrator of division to get blur2
+            # as per Stephan Ludwig et al 2019
+            # should normalise blur2 and decon each time because numbers get big and we risk overflow when multiplying in next step
+            conv1norm = conv1 / (tf.math.reduce_sum(conv1))
+            datatNorm = datat / (tf.math.reduce_sum(datat))
+            deltaParam = 0.001
+            ratio = (datatNorm + deltaParam) / (conv1norm + deltaParam)
+            #blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
+            #ratioNorm = ratio / (tf.math.reduce_sum(ratio))
+            #deconNorm = decon / (tf.math.reduce_sum(decon))
+            # decon is the  normalised blurred model multiplied by the model
+            # Positivity constraint on result for iteration
+            decon = tf.maximum(decon * ratio, 0.)
+            #decon = tf.maximum(decon * conv2, 0.)
+            
+            # TODO - Smoothing every 5 iterations with gaussian or wiener. 
+            # TODO rescale back to input data sum intensity. , probably need to adjust deltaParam too. 
 
             # If given an "observer", pass the current image restoration and iteration counter to it
             if self.observer_fn is not None:

From 43f01b4101061d4888f2541087b2f56e55692b2b Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Wed, 29 Apr 2020 16:17:34 +0200
Subject: [PATCH 11/15] add normalisation or raw image to make SSIM convergence
 measurement work

---
 .gitignore                                        |  1 +
 .../scripts/simpleFlowDecAlgoConvergenceTest.py   | 15 +++++++++++----
 2 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/.gitignore b/.gitignore
index 7e4e5ab..801764c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -97,3 +97,4 @@ python/examples/scripts/FlowDecLogGold50.txt
 python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif250iterations.tif
 python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif4iterations.tif
 python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif5iterations.tif
+python/examples/scripts/resultC1-YeastTNA1_1516_conv_RG_26oC_003_256xcropSub100.tifgpsf_3D_1514_a3_001_WF-sub105crop64.tif10iterations.tif
diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 7c871e4..7e64abf 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -39,6 +39,7 @@
 startImports = time.process_time()   
 from skimage.external.tifffile import imsave, imread
 from skimage.metrics import structural_similarity as ssim
+#from tensorflow.image import ssim as ssim
 from flowdec import data as fd_data
 from flowdec import restoration as fd_restoration
 importsTime = (time.process_time() - startImports)
@@ -77,9 +78,12 @@
 # is generally not useful for anything other than development or debugging
 # imgs = []
 def observer(decon_crop, i, decon_pad, conv1, *args):
+    #normalise the raw data so its sum is 1
+    rawNorm = raw / raw.sum()
+    #the sum should be 1
+    sumRaw = rawNorm.sum()
     #imgs.append(decon_crop)
-    if i % 10 == 0:
-        sumRaw = raw.sum()
+    if i % 1 == 0:
         sumBlurredModel = decon_crop.sum()
 	    # compute convergence residuals between raw image and blurred model image
 		# with current test data these get bigger not smaller with iterations... weird..
@@ -87,9 +91,12 @@ def observer(decon_crop, i, decon_pad, conv1, *args):
         convergenceR = convergenceResiduals / sumRaw
 		# let's try structural similarity (as it's supposed to be better then mean square error)
 		# which indeed seems to track convergence of this dataset in a way that seems to match the image results. 
-        structSim = ssim(raw, decon_crop, data_range=decon_crop.max() - decon_crop.min())
+        #SSIM in skimage:
+        structSim = ssim(rawNorm, decon_crop, data_range=decon_crop.max() - decon_crop.min())
+        #SSIM in TensorFlow should be faster:  imported ft.image.ssim as ssim but this doesnt work... fix if too slow otherwize. 
+        # structSim = ssim(raw, decon_crop, max_val=1.0, filter_size=11, filter_sigma=1.0, k1=0.01, k2=0.03)
         print('Iter,{},RawSum,{:.3f},DeconSum,{:.3f},DeconMax,{:.3f},DeconStDev,{:.3f},SumResiduals,{:.3f},ConvergeR,{:.16f},SSIM,{:.3f}'.format(
-		   i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
+            i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
 
 # Run the deconvolution process and note that deconvolution initialization is best kept separate from 
 # execution since the "initialize" operation corresponds to creating a TensorFlow graph, which is a 

From 972647cc528b9b6d61feb531eec7e58b960d38c3 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Wed, 29 Apr 2020 16:19:10 +0200
Subject: [PATCH 12/15] back off iterations for testing

---
 python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 7e64abf..5215250 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -67,7 +67,7 @@
 kernel = imread(PSF)
 
 #base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
-base_iter = 250
+base_iter = 5
 
 
 # Create an observer function to monitor convergence, 

From 48f10d01d1f9bef0e2c95c05de7f0ff0e97c5df3 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Wed, 29 Apr 2020 16:21:14 +0200
Subject: [PATCH 13/15] WIP: add gaussian kernel generation for gaussian
 smoothing operations with tf.nn.conv3d

---
 python/flowdec/restoration.py | 43 ++++++++++++++++++++++++++++++-----
 1 file changed, 37 insertions(+), 6 deletions(-)

diff --git a/python/flowdec/restoration.py b/python/flowdec/restoration.py
index 099e453..7da1d9c 100755
--- a/python/flowdec/restoration.py
+++ b/python/flowdec/restoration.py
@@ -291,8 +291,29 @@ def cond(i, decon):
 
         def conv(inputData, kernel_fft):
             return tf.math.real(fft_rev(fft_fwd(tf.cast(inputData, self.fft_dtype)) * kernel_fft))
+			
+        def gaussian_kernel(size: int,
+                    mean: float,
+                    std: float,
+                    ):
+            """Makes 3D gaussian Kernel for convolution."""
 
-        def body(i, decon):
+            d = tf.distributions.Normal(mean, std)
+
+            vals = d.prob(tf.range(start = -size, limit = size + 1, dtype = tf.float32))
+
+            gauss_kernel = tf.einsum('i,j,k->ijk',
+                                        vals,
+                                        vals,
+                                        vals)
+            # return the kernel normalised to sum =1
+            return gauss_kernel / tf.reduce_sum(gauss_kernel)
+
+        gaussKernel = gaussian_kernel(9, 1.0, 7.0)
+        # Expand dimensions of `gauss_kernel` for `tf.nn.conv3d` signature.
+        gaussKernel = gaussKernel[:, :, :, tf.newaxis, tf.newaxis, tf.newaxis]
+
+        def body(i, decon,):
             '''# Richardson-Lucy Iteration - logic taken largely from a combination of
             # the scikit-image (real domain) and DeconvolutionLab2 implementations (complex domain)
             # conv1 is the current model blurred with the PSF
@@ -312,8 +333,8 @@ def body(i, decon):
             # conv1 is the current model blurred with the PSF
             conv1 = conv(decon, kern_fft)
 
-            # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)
-            # we wont do it here as we will use the delta values in denom and numerrator of division to get blur2
+            # High-pass filter to avoid division by very small numbers (see DeconvolutionLab2)?
+            # we wont do it here as we will use the delta parameter in denom and numerrator of division to get blur2
             # as per Stephan Ludwig et al 2019
             # should normalise blur2 and decon each time because numbers get big and we risk overflow when multiplying in next step
             conv1norm = conv1 / (tf.math.reduce_sum(conv1))
@@ -326,16 +347,26 @@ def body(i, decon):
             # decon is the  normalised blurred model multiplied by the model
             # Positivity constraint on result for iteration
             decon = tf.maximum(decon * ratio, 0.)
-            #decon = tf.maximum(decon * conv2, 0.)
+			# Smooth the intermediate result image with Gaussian of sigma 1 every 5th iteration
+            # to control noise buildup that Gold method is succeptible to.
+            # Use tf.nn.conv3d to convolve a Gaussian kernel with an image:
+            # Make Gaussian Kernel with desired specs using gaussian_kernel function defined above
+            if i % 5 == 0:
+                # Convolve decon with gauss kernel.
+                tf.nn.conv3d(decon, filter=gaussKernel, strides=[1, 1, 1, 1, 1], padding="SAME")
+            # normalise the result so the sum of the data is 1
+            decon = decon / (tf.math.reduce_sum(decon))
             
             # TODO - Smoothing every 5 iterations with gaussian or wiener. 
-            # TODO rescale back to input data sum intensity. , probably need to adjust deltaParam too. 
+            # TODO rescale back to input data sum intensity -  probably need to adjust deltaParam too. 
 
             # If given an "observer", pass the current image restoration and iteration counter to it
             if self.observer_fn is not None:
                 # Remove any cropping that may have been added as this is usually not desirable in observers
                 decon_crop = unpad_around_center(decon, tf.shape(datah))
-                # we can use these capurured observed tensors to evaluate eg convergence
+                # normalise the result so the sum of the data is 1
+                decon_crop = decon_crop / (tf.math.reduce_sum(decon_crop))
+				# we can use these captured observed tensors to evaluate eg convergence
                 # in eg. the observer function used.
                 _, i, decon, conv1  = tf_observer(
 				    [decon_crop, i, decon, conv1], self.observer_fn)

From 249dbb6e349c1f82e9c8cc6d424b4629acb3fd29 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Wed, 29 Apr 2020 16:21:58 +0200
Subject: [PATCH 14/15] update deltaParam to a value just large enough  that
 seem to work well

---
 python/flowdec/restoration.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/python/flowdec/restoration.py b/python/flowdec/restoration.py
index 7da1d9c..8a5b625 100755
--- a/python/flowdec/restoration.py
+++ b/python/flowdec/restoration.py
@@ -339,7 +339,8 @@ def body(i, decon,):
             # should normalise blur2 and decon each time because numbers get big and we risk overflow when multiplying in next step
             conv1norm = conv1 / (tf.math.reduce_sum(conv1))
             datatNorm = datat / (tf.math.reduce_sum(datat))
-            deltaParam = 0.001
+            # this value seems to work well fo rthe images that are normalised to sum of 1
+            deltaParam = 1e-4
             ratio = (datatNorm + deltaParam) / (conv1norm + deltaParam)
             #blur1 = tf.where(conv1 < self.epsilon, tf.zeros_like(datat), datat / conv1, name='blur1')
             #ratioNorm = ratio / (tf.math.reduce_sum(ratio))

From 832aa04957f9159d3915fe48f8d86e56562e6c88 Mon Sep 17 00:00:00 2001
From: dan white <dan@chalkie.org.uk>
Date: Tue, 23 Jun 2020 09:37:23 +0200
Subject: [PATCH 15/15] work in progress with Gauss / Wiener  filtering in
 flowdec for Acellerated gold algo.

---
 python/examples/scripts/gaussKernTest.py      | 47 +++++++++++++++++++
 .../simpleFlowDecAlgoConvergenceTest.py       | 19 ++++----
 2 files changed, 57 insertions(+), 9 deletions(-)
 create mode 100644 python/examples/scripts/gaussKernTest.py

diff --git a/python/examples/scripts/gaussKernTest.py b/python/examples/scripts/gaussKernTest.py
new file mode 100644
index 0000000..04ded17
--- /dev/null
+++ b/python/examples/scripts/gaussKernTest.py
@@ -0,0 +1,47 @@
+import tensorflow as tf
+import sys
+
+
+def gaussian_kernel(size: int, mean: float, std: float):
+	d = tf.distributions.Normal(mean, std)
+	vals = d.prob(tf.range(size, dtype = tf.float32))
+	gauss_kernel = tf.einsum('i,j,k->ijk',
+							vals,
+							vals,
+							vals)             
+	# normalise gauss kernel to sum = 1
+	kerngausnorm = gauss_kernel / tf.reduce_mean(gauss_kernel)
+	return kerngausnorm
+
+def squareIt(tensor):
+    tensorpow2 = tf.multiply(tensor, tensor)
+    return tensorpow2
+
+sess = tf.compat.v1.Session()
+with sess.as_default():
+    tensor = gaussian_kernel(5, 2.0, 1.0)
+    tensor2 = gaussian_kernel(5, 3.0, 1.0)
+    tensor = tf.dtypes.cast(tensor, tf.float32)
+    tensor2 = tf.dtypes.cast(tensor2, tf.float32)
+    tensorSquared = squareIt(tensor)
+    tensorPlusTen = tensor + 10.0
+    # expand its dimensionality to fit into conv3d, input and filter have different dimension orders, filter has no batch but in and out channels as 2 last dims
+    tensor_expand = tf.expand_dims(tensor, 0)
+    tensor_expand = tf.expand_dims(tensor_expand, -1)
+    tensor_filter = tf.expand_dims(tensor2, -1)
+    tensor_filter = tf.expand_dims(tensor_filter, -1)
+    # why does tf.nn.conv3d output a tensor containing only 1 value???eg  [[[[[]35234.325]]]] is it the sum of the real 3D image output? 
+    tensorConvolved = tf.compat.v1.nn.conv3d(tensor_expand, filter=tensor_filter, strides=[1,1,1,1,1], padding="VALID", data_format='NDHWC')
+    print_op = tf.print("3D Gauss Kernel? :\n", tensor, " Sum = ", tf.reduce_sum(tensor), " Max = ", tf.reduce_max(tensor), "\n",
+                        "3D Gauss Kernel squared? :\n", tensorSquared, " Sum = ", tf.reduce_sum(tensorSquared), " Max = ", tf.reduce_max(tensorSquared), "\n",
+                        "3D Gauss Kernel convolved? :\n", tensorConvolved, " Sum = ",
+                        tf.reduce_sum(tensorConvolved), " Max = ", tf.reduce_max(tensorConvolved), " DimsShape = ", tf.shape(tensorConvolved), "\n",
+                        " Plus10 :\n", tensorPlusTen,
+                        output_stream=sys.stdout)
+    with tf.control_dependencies([print_op]):
+        the_tensor = tensor * 1.0
+        the_tensor2 = tensor2 * 1.0
+        the_tensorpow2 = tensorSquared * 1.0
+        the_tensorConvolved = tensorConvolved * 1.0
+    sess.run(the_tensorConvolved)
+
diff --git a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
index 5215250..5446b86 100644
--- a/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
+++ b/python/examples/scripts/simpleFlowDecAlgoConvergenceTest.py
@@ -32,10 +32,6 @@
 import sys
 #import logging
 
-#send std output to a log file
-sys.stdout = open('FlowDecLogGold50.txt', 'w')
-#logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
-
 startImports = time.process_time()   
 from skimage.external.tifffile import imsave, imread
 from skimage.metrics import structural_similarity as ssim
@@ -67,7 +63,7 @@
 kernel = imread(PSF)
 
 #base number of iterations - RL converges slowly so need tens of iterations or maybe hundreds. 
-base_iter = 5
+base_iter = 15
 
 
 # Create an observer function to monitor convergence, 
@@ -77,13 +73,13 @@
 # and the remaining arguments (at TOW) only include the uncropped image result which
 # is generally not useful for anything other than development or debugging
 # imgs = []
-def observer(decon_crop, i, decon_pad, conv1, *args):
+def observer(decon_crop, i, decon, conv1, kerngaussh, *args):
     #normalise the raw data so its sum is 1
     rawNorm = raw / raw.sum()
     #the sum should be 1
     sumRaw = rawNorm.sum()
     #imgs.append(decon_crop)
-    if i % 1 == 0:
+    if i % 5 == 0:
         sumBlurredModel = decon_crop.sum()
 	    # compute convergence residuals between raw image and blurred model image
 		# with current test data these get bigger not smaller with iterations... weird..
@@ -95,8 +91,8 @@ def observer(decon_crop, i, decon_pad, conv1, *args):
         structSim = ssim(rawNorm, decon_crop, data_range=decon_crop.max() - decon_crop.min())
         #SSIM in TensorFlow should be faster:  imported ft.image.ssim as ssim but this doesnt work... fix if too slow otherwize. 
         # structSim = ssim(raw, decon_crop, max_val=1.0, filter_size=11, filter_sigma=1.0, k1=0.01, k2=0.03)
-        print('Iter,{},RawSum,{:.3f},DeconSum,{:.3f},DeconMax,{:.3f},DeconStDev,{:.3f},SumResiduals,{:.3f},ConvergeR,{:.16f},SSIM,{:.3f}'.format(
-            i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max()))
+        print('Iter,{},RawSum,{:.3f},DeconSum,{:.3f},DeconMax,{:.3f},DeconStDev,{:.3f},SumResiduals,{:.3f},ConvergeR,{:.16f},SSIM,{:.3f},KerGMax,{:.3f}'.format(
+            i, sumRaw, sumBlurredModel, decon_crop.max(), decon_crop.std(), convergenceResiduals.max(), convergenceR.max(), structSim.max(), kerngaussh.max()))
 
 # Run the deconvolution process and note that deconvolution initialization is best kept separate from 
 # execution since the "initialize" operation corresponds to creating a TensorFlow graph, which is a 
@@ -118,6 +114,11 @@ def observer(decon_crop, i, decon_pad, conv1, *args):
 multiRunFactor = 1
 timingListIter = []
 timingListTime = []
+
+#send std output to a log file
+sys.stdout = open('FlowDecLogGold' + str(base_iter) + 'multi' + str(multiRunFactor) + 'GAUSS.txt', 'w')
+#logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
+
 for i in range(1, multiRunFactor+1):
   niter = (base_iter*i)
   # start measuring time