Phantom_generate/synthetic_data.py

@@ -16,8 +16,8 @@ the blurred images. N is the run number.
 output_path = "synthetic-data/"
 output_height = 2192
 output_width = 3840
-generation_number = 100  # number of files to output
-
+generation_number = 20  # number of files to output
+dtype = np.float16
 line_thickness = (4,15)*255
 
 line_color = (0.6,0.98)
@@ -88,10 +88,10 @@ if blur_kernel == "gaussian":
 
 if blur_kernel == "custom":
     kernel = np.load(kernel_path)
-
+    print("kernel size: {}".format(kernel.shape))
 storage = None
 if output_single_npz:
-    storage = np.zeros((2,generation_number,output_height,output_width),dtype=np.float64)
+    storage = np.zeros((2,generation_number,output_height,output_width),dtype=dtype)
 
 for i in range(generation_number):
     gt = gen_image(output_height, output_width, background_color,line_thickness, line_color, line_count,square_count)
@@ -100,6 +100,9 @@ for i in range(generation_number):
     blured = blured + np.random.rand(blured.shape[0],blured.shape[1]) * random_noise
     blured = np.clip(blured,a_min=0,a_max=1)
 
+    #convert array to desired floating point precision
+    gt = gt.astype(dtype)
+    blured = blured.astype(dtype)
     if interactive:
         show_image(gt, "gt")
         show_image(blured, "blured")

README.md

@@ -27,7 +27,7 @@ See documentation for more details. There is one container currently provided fo
 ## Usage
 
 ### Synthetic Data
-There is a python [script](Phantom_generate/synthetic_data.py) to generate synthetic training data and blur it using the PSF function of the camera. See the python script for more details. There are also helper scripts available to determine the PSF of the camera from images of pinhole light sources. See [this](helper%20scripts/PSF/README.md) for more details. 
+There is a python [script](Phantom_generate/synthetic_data.py) to generate synthetic training data and blur it using the PSF function of the camera. See the python script for more details. There are also helper scripts available to determine the PSF of the camera from images of pinhole light sources. See [this](helper scripts/PSF/VC MIPI/README.md) for more details. 
 
 ### Training
 

RLN_code/RLN_single_Train.py

@@ -218,6 +218,11 @@ for epoch in range(0, epochs):
                 estimatation1 = estimatation1*std + mean
                 update1 = update1*std + mean
 
+                prediction = tf.clip_by_value(prediction,0,1,name="final_activation_clip")
+                estimatation1 = tf.clip_by_value(estimatation1,0,1,name="final_activation_clip")
+                update1 = tf.clip_by_value(update1,0,1,name="final_activation_clip")
+
+
                 tf.summary.image("test output", prediction, i,max_outputs=1)
                 tf.summary.image("test estimation", estimatation1, i,max_outputs=1)
 

RLN_code/RLN_single_inference.py

@@ -56,6 +56,7 @@ for path_i in a:
     x = (x - mean) / std
     x_i = model(x, training=False)
     x_i = x_i[0] * std + mean
+    x_i = tf.clip_by_value(x_i, 0, 1, name="final_activation_clip")
 
     et = time.time()
     dt = np.array(et - st)

helper scripts/PSF/Basler/PSF.py

+import pickle
+
+from helpers.image import show_image,read_y10
+from helpers.output import save_array,load_array
+import os
+import cv2
+import numpy as np
+import matplotlib.pyplot as plt
+
+"""
+this script is for actually determining the PSF for each measurement. For each measurement it will do:
+
+1) It finds the brightest point in the image and cuts a square with (2*delta+1) side length around the point.
+2) It will show the whole image and teh cut portion of the image. You can accept by [ENTER] or press c to cancel
+3) It will show the normed version of the cropped image (L1 norm). You can accept by [ENTER] or press c to cancel
+4) It will save the image
+"""
+
+delta = 3# box size
+detail_display = False # display the position of the brightest point on the whole image
+root = ("/home/guney/Project/Richardson-Lucy-Net/helper scripts/data/psf2")
+dirs = os.listdir(root)
+for i,path in enumerate(dirs): # iterate over the measurements
+    path = root + "/" + path+"/"
+
+    print(path)
+    try:
+        image = load_array(path+"avg_image.npz")
+    except FileNotFoundError:
+        continue
+
+    max_i = np.unravel_index(np.argmax(image),image.shape) # find the brightest point
+
+    print("Maximal Brightness at {}".format(max_i))
+
+    y_min, y_max = max_i[0]-delta,max_i[0]+delta+1 #determine the box around the point
+    x_min, x_max = max_i[1]-delta,max_i[1]+delta+1
+
+    top_left_point = (x_min,y_min)
+    bottom_right_point = (x_max,y_max)
+    if detail_display:
+        im_2 = np.copy(image)
+        cv2.rectangle(im_2,top_left_point,bottom_right_point,1.0,4)
+        show_image(im_2,"avg",norm=(0,1))
+    cropped = image[y_min:y_max, x_min:x_max]
+    #print(np.unravel_index(np.argmax(cropped),cropped.shape))
+    show_image(cropped,"avg_cropped",norm=(0,4))
+    #print(cropped.shape)
+    if cv2.waitKey() == ord('c'):
+        print("exiting")
+        cv2.destroyAllWindows()
+        exit(0)
+
+    avg_image = cropped
+
+
+    normed = np.zeros(np.shape(avg_image), dtype=np.float64)
+    cv2.normalize(avg_image,normed,1.0,0,cv2.NORM_L1) # normalize to L! norm
+
+    show_image(normed,"avg_L1",norm=(0,3))
+
+    plt.imshow(normed)
+    plt.savefig(path+"psf.png")
+    plt.show()
+
+    print("L1 norm is {}".format(np.sum(normed)))
+    if cv2.waitKey() == ord('c'):
+        print("exiting")
+        cv2.destroyAllWindows()
+        exit(0)
+
+    plt.figure()
+    plt.imshow(normed)
+    plt.savefig(path+"psf.png")
+    plt.show()
+
+    crop_limits = {"y_max" : y_max, "y_min" : y_min, "x_max" : x_max, "x_min" :x_min}
+    with open(path+'crop_limits.pickle', 'wb') as handle:
+        pickle.dump(crop_limits, handle, protocol=pickle.HIGHEST_PROTOCOL)
+    #cv2.imwrite(path+"psf.png",normed)
+    np.save(path+"cropped.npy",cropped)
+    np.save(path+"psf.npy",normed)
+
+    cv2.destroyAllWindows()
\ No newline at end of file

helper scripts/PSF/Basler/merge_PSF.py

+import os
+
+import numpy as np
+from matplotlib import pyplot as plt
+import cv2
+
+"""
+This script is used to merge the outputs of PSF.py. It merges the outputs by averaging and plots the individual PSFs.  
+"""
+root = ("/home/guney/Project/Richardson-Lucy-Net/helper scripts/data/psf2/")
+dirs = os.listdir(root)
+b = []
+# fig, axs = plt.subplots(2, 4)
+for i,path in enumerate(dirs):
+    # plot_index = np.unravel_index(i-1,axs.shape)
+    path = root + path+"/"
+    try:
+        a = np.load(path+"cropped.npy")
+    except FileNotFoundError:
+        continue
+    print(a.shape)
+    b.append(a)
+
+    # axs[plot_index].set_title("PSF of {}".format(i))
+    # axs[plot_index].imshow(a)
+
+# fig.show()
+
+b = np.array(b)
+b = np.mean(b,axis=0)
+normed = np.zeros(np.shape(b), dtype=np.float64)
+cv2.normalize(b, normed, 1.0, 0, cv2.NORM_L1)
+
+fig,ax = plt.subplots(1)
+ax.set_title("Point Spread Function of the camera")
+heatmap = ax.pcolor(normed)
+fig.colorbar(heatmap,ax=ax)
+fig.show()
+
+print("L1 Norm is {}".format(np.sum(normed)))
+np.save("avg_psf",normed)
+print(b.shape)

helper scripts/PSF/Basler/preprocess_PSF.py

+from helpers.image import show_image,read_y10
+
+import os
+import cv2
+import numpy as np
+import matplotlib.pyplot as plt
+
+from helpers.output import save_array
+"""
+This script reads PSF measurements in Y10 format and averages the measurements for each position.  I.E: It averages over all y10 images in a folder and outputs the average image as npz.
+"""
+path = ("/home/guney/Project/Richardson-Lucy-Net/helper scripts/data/psf2")
+print(os.getcwd())
+b = os.listdir(path)
+for i in b:
+
+    pathi = path + "/" + i+"/"
+    print(pathi)
+    a = os.listdir(pathi)
+
+    a = [pathi +i for i in a if i.split('.')[-1] == 'bmp']
+    if len(a) == 0:
+        continue
+    image = cv2.imread(a[0], cv2.IMREAD_GRAYSCALE)
+
+    # show_image(image,i)
+    # cv2.waitKey(0)
+
+    datas = np.zeros((len(a),image.shape[0],image.shape[1]),dtype=np.float64)
+
+    for i ,path_i in enumerate(a):
+        loaded = cv2.imread(path_i, cv2.IMREAD_GRAYSCALE)
+        datas[i,:,:] = loaded / ((2**8) - 1)
+        #datas[i,:,:] = loaded
+    print(datas.shape)
+    avg_image = np.mean(datas,axis=0,dtype=np.float64)
+    print(avg_image.dtype,avg_image.shape)
+    print(np.max(avg_image),np.mean(avg_image))
+    # show_image(avg_image,"avg")
+    # cv2.waitKey(0)
+
+
+    save_array(pathi+"avg_image.npz",avg_image)

helper scripts/PSF/VC MIPI/README.md

+## How to Determine the PSF
+The PSF is used to generate synthetic data to be used during the training process. Using the PSF we can model the effects of the camera on real world objects. Thus we can generate realistic input images from ground truths.
+
+### Measurement
+We did the measurements using a laser and A 1.4 um pinhole. This was barely smaller than the pixel size of the CCD. Thus we used the lens in 1:1 magnification mode at ca. 50 cm distance from the pinhole. We took 50 images per measurement and did 8 measurements with the light shining on different parts of the lens.
+
+### Analysis
+The code is divided into three parts:
+1. preprocess_PSF.py: This combines the 50 pictures into a single averaged picture for each measurement. Done only once during the process.
+2. PSF.py: This script actually calculates the PSF from the averaged images created in the last step. This is done by cropping the image to a (2\*delta+1 x 2\*delta+1) square around the brightest point and normalizing the output according to the L1 norm.
+3. merge_PSF.py: This combines the output of the PSF.py for each measurement into a combined PSF and generates graphs. 
+
+Read the relevant scripts for more details
\ No newline at end of file

helper scripts/PSF/merge_PSF.py → helper scripts/PSF/VC MIPI/merge_PSF.py

@@ -24,10 +24,12 @@ b = np.array(b)
 b = np.mean(b,axis=0)
 normed = np.zeros(np.shape(b), dtype=np.float64)
 cv2.normalize(b, normed, 1.0, 0, cv2.NORM_L1)
-plt.figure()
-plt.title("Average PSF")
-plt.imshow(normed)
-plt.show()
+
+fig,ax = plt.subplots(1)
+ax.set_title("Point Spread Function of the camera")
+heatmap = ax.pcolor(normed)
+ax.colorbar(heatmap)
+fig.show()
 
 print("L1 Norm is {}".format(np.sum(normed)))
 np.save("avg_psf",normed)

helper scripts/image analysis/read_inference_output_from_Server.py

@@ -10,30 +10,42 @@ This script is used to display the output of the inference that runs on the serv
 RLN_single_inference.py outputs .npz files with the arrays "x" for input, "dt" for calculation time, and "f_x" for 
 inference output. This display x and f_x for all .npz s in the path.
 """
-path = "/home/guney/Project/Richardson-Lucy-Net/data from server/traind model inference/"
+path = "/home/guney/Project/Richardson-Lucy-Net/data from server/traind model 2d/inference_output_2/"
 a = os.listdir(path)
-a = [path+i for i in a if i.split('.')[-1] == 'npz']
-print(a)
+a = [i for i in a if i.split('.')[-1] == 'npz']
+a = sorted(a,key = lambda x : int(x.split(".")[0]))
+a = [path + i for i in a]
 display = True
+delta = 800
 
 
 
 
-
-for data in a:
+for i,data in enumerate(a):
     print("==========================================================")
     print(data)
     print("==========================================================")
 
     data = np.load(data,allow_pickle=True)
-    inference = data["f_x"][0]
-    input = data["x"]
-    print(inference.shape)
+    inference = data["f_x"][0][:,:,0]
+    input = data["x"][0,:,:,0]
+    print(np.max(inference))
+    print(inference.dtype)
+    shape = inference.shape
+    print("Read: max: {} , min: {}, mean: {}".format(np.max(inference),np.min(inference),np.mean(inference)))
+    inference = np.clip(inference,0,1)
+    inference = inference[shape[0]//2-delta:shape[0]//2+delta,shape[1]//2-delta:shape[1]//2+delta]
+    # cv2.normalize(inference, inference, 0, 1.2, cv2.NORM_MINMAX)
+    input = input[shape[0]//2-delta:shape[0]//2+delta,shape[1]//2-delta:shape[1]//2+delta]
+    # cv2.normalize(input, input, 0, 1,cv2.NORM_MINMAX)
+    # cv2.imwrite("{}_inference.jpg".format(i),inference)
+    # cv2.imwrite("{}_input.jpg".format(i),input)
+    print(inference.dtype)
 
     if display:
-        show_image(inference,"inference")
+        show_image(inference,"inference",500,500)
         cv2.moveWindow("inference",900,000)
-        show_image(input,"input")
+        show_image(input,"input",500,500)
         print(np.max(inference),np.min(inference))
 
 
@@ -42,7 +54,6 @@ for data in a:
             cv2.destroyAllWindows()
             exit(0)
 
-
     print("==========================================================")
 
 

helper scripts/image analysis/view_y10.py

@@ -5,7 +5,7 @@ import cv2
 boilerplate script to view y10 files
 """
 from helpers.image import display_y10
-path = "/media/guney/9E79-C5A1/grid/3/"
+path = ("/media/guney/9E79-C5A1/photos/ind2/")
 
 a = os.listdir(path)
 a = [i for i in a if i.split('.')[-1] == "y10"]
@@ -16,5 +16,6 @@ for i in a:
     cv2.moveWindow(i,900,000)
     display_y10(path+i,"{} (normalized)" .format(i))
 
-    cv2.waitKey(0)
-    cv2.destroyAllWindows()
+    if cv2.waitKey(0) == ord('c'):
+        cv2.destroyAllWindows()
+        exit(0)

helper scripts/image analysis/y10_to_npz.py

@@ -10,22 +10,28 @@ This script is used to prepare the y10 files from the camera for the neural netw
 converts it into float32 representation, normalizes it to (0,1) and inverts the image. The network works best on 
 white features on black background.
 """
-
-path = ("/media/guney/9E79-C5A1/grid/3/")
-output= "test/"
+i = 1
+path = ("/media/guney/9E79-C5A1/photos/ind2/")
+output= "test/photosind2/"
 if not os.path.exists(output):
     os.makedirs(output)
 print(os.getcwd())
 a = os.listdir(path)
-a = [i for i in a if i.split('.')[-1] == 'y10']
+a = [i for i in a if i.split('.')[-1] == 'y10'][:]
 a = sorted(a,key = lambda x : int(x.split(".")[0]))
 
 for i,path_i in enumerate(a):
     print(path_i)
     datat = (read_y10(path+path_i)/ ((2**10) - 1))
+    print("Read: max: {} , min: {}, mean: {}".format(np.max(datat),np.min(datat),np.mean(datat)))
+
+
     cv2.normalize(datat,datat,0,1,cv2.NORM_MINMAX)
+    show_image(datat,"image_normed")
+
     datat = 1- datat
-    print(np.max(datat),np.min(datat))
+    print("Out: max: {} , min: {}, mean: {}".format(np.max(datat),np.min(datat),np.mean(datat)))
     np.savez_compressed(output+str(i) + ".npz",datat)
-    #show_image(datat,"image")
-    #cv2.waitKey(0)
+    show_image(datat,"image")
+    if cv2.waitKey(0) == ord('c'):
+        exit(0)
\ No newline at end of file