NCS Demo Simulation (OpenCL).

In order for this to work you need the following in your Python path:

The reference NCS/python3-6 Python implementation.
The NCS/opencl Python module.



In [ ]:

    
import matplotlib
import matplotlib.pyplot as pyplot
import numpy
import os
import time

# python3-6 NCS.
import pyNCS
import pyNCS.denoisetools as ncs

# python3 and OpenCL 
import pyOpenCLNCS.ncs as ncsOC

# Generate the same random noise each time.
numpy.random.seed(1)



In [ ]:

    
py_ncs_path = os.path.dirname(os.path.abspath(pyNCS.__file__))
print(py_ncs_path)

pyNCS analysis

This is a basically a copy of NCS/python3-6/NCSdemo_simulation.py



In [ ]:

    
# create normalized ideal image
fpath1 = os.path.join(py_ncs_path, "../randwlcposition.mat")
imgsz = 128
zoom = 8
Pixelsize = 0.1
NA = 1.4
Lambda = 0.7
t = time.time()
res = ncs.genidealimage(imgsz,Pixelsize,zoom,NA,Lambda,fpath1)
elapsed = time.time()-t
print('Elapsed time for generating ideal image:', elapsed)
imso = res[0]
pyplot.imshow(imso,cmap="gray")
    
# select variance map from calibrated map data
fpath = os.path.join(py_ncs_path, "../gaincalibration_561_gain.mat")        
noisemap = ncs.gennoisemap(imgsz,fpath)
varsub = noisemap[0]*10 # increase the readout noise by 10 to demonstrate the effect of NCS algorithm
gainsub = noisemap[1]
    
# generate simulated data
I = 100
bg = 10
offset = 100
N = 50
dataimg = ncs.gendatastack(imso,varsub,gainsub,I,bg,offset,N)
imsd = dataimg[1]



In [ ]:

    
# generate noise corrected image
Rs = 8
iterationN = 15
alpha = 0.1    

out_name = os.path.join(py_ncs_path, "../../out.npy")

# This is useful for debugging as it takes a long time for this approach to
# to reduce the noise of an image. Once you've done this once you can just
# load the reference result.
if not os.path.exists(out_name):
    
    # This takes ~24 seconds on my laptop.
    out = ncs.reducenoise(Rs,imsd[0:1],varsub,gainsub,imgsz,Pixelsize,NA,Lambda,alpha,iterationN)
    numpy.save(out_name, out)

else:
    out = numpy.load(out_name)
    
print(out.shape)



In [ ]:

    
f,(ax1,ax2) = pyplot.subplots(1,2,sharey=False,figsize = (8,8))
ax1.imshow(imsd[0],aspect='equal',cmap="gray")
ax2.imshow(out[0],aspect ='equal',cmap="gray")
pyplot.show()

OpenCL analysis

Mixed OpenCL and Python NCS analysis.



In [ ]:

    
# Make an OTF mask that is the correct size.
rcfilter = ncs.genfilter(16,Pixelsize,NA,Lambda,'OTFweighted',1,0.7)

print(rcfilter.shape)

pyplot.imshow(rcfilter, cmap = "gray")
pyplot.show()



In [ ]:

    
# Calculate gamma and run Python/C NCS.
gamma = varsub/(gainsub*gainsub)



In [ ]:

    
# Create a list with all 50 images.
ims = []
for i in range(imsd.shape[0]):
    ims.append(numpy.copy(imsd[i]))



In [ ]:

    
# Analyze all 50 images at once.
#
# This takes about 1 second on my laptop with an Intel (integrated?) GPU, so
# a 24 x 50 = 1200X speed improvement.
#
out_oc = ncsOC.reduceNoise(ims, gamma, rcfilter, alpha, verbose = True)

Compare results to reference implementation.

The differences are larger here because we are using a different sub-region size (16x16 vs 10x10 for the reference). If we compare against the C version, also with 16x16 sub-region size the results are much closer.



In [ ]:

    
f,(ax1,ax2) = pyplot.subplots(1,2,sharey=False,figsize = (8,8))
ax1.imshow(out[0],aspect='equal',cmap="gray")
ax2.imshow(out_oc[0],aspect ='equal',cmap="gray")
pyplot.show()



In [ ]:

    
pyplot.figure(figsize = (6,6))
pyplot.imshow(out[0] - out_oc[0], cmap = "gray")
pyplot.show()



In [ ]:

    
print("Maximum pixel difference is {0:.3f}e-".format(numpy.max(numpy.abs(out[0] - out_oc[0]))))



In [ ]: