scattersim/spectroscopy.py at master · STIM-Lab/scattersim · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
import numpy as np
import os
import subprocess
import glob
import skimage as ski
import skimage.io
import matplotlib.pyplot as plt
import matplotlib
import ffmpeg
from tqdm import tqdm
import shutil
from ScatterView import *
import math
import time

scattervol_path = "/home/ruijiao/A_research/build/scattersim/"
source_path = "/home/ruijiao/A_research/scattersim/"      # Where the sample is from
data_path = source_path + "spectroscopy_even/"
sample_name = "3.jpg"


# Configuration
arr = np.loadtxt(source_path+"polystyrene_all_data_standard_format.csv",
                 delimiter=",", dtype=str)      # refractive index form. Format: wavelength, real part of refractive index, imag part of refractive index
# lambdas = [2.038, 9.3254, 13.227]
# n = [1.576+0.0000369j, 1.5412+0.0306j, 1.52708+0.105j]
size = [500, 500, 5]                       # Size of the sample. 4 is the diamter of the cylinder(z)
coefs = [40, 40]                          # Pixel number: 40*40
resolution = 8
save_axis = 2
relative_slice = size[2]/2           # Use the bottom of the sample
sample_npy = data_path + "3_npy.npy"
ref_npy = data_path + "3_npy_ref.npy"
recon_A = []

for i in tqdm(range(len(arr[:, 0]))):
    # Define file name
    ref_cw = data_path + arr[i, 0] + "_ref.cw"
    result_cw = data_path + arr[i, 0] + ".cw"
    out_ref_npy = data_path + arr[i, 0] + "field_3d_ref.npy"
    out_result_npy = data_path + arr[i, 0] + "field_3d.npy"
    if i%2!=0:
        continue
    # Define wave
    lambdai = float(arr[i, 0])
    ni = complex(float(arr[i, 1]), float(arr[i, 2]))
    start = time.time()
    print("")
    print("----------"+ str(i+1) + "th loop-------------")
    print("wavelength:" + str(lambdai))
    print("n:" + str(ni))
    # result_cw = scattervol_path + str(lambdas[lambdai]) + ".cw"

    # Create sample
    img = ski.io.imread(source_path + sample_name)[:, :, 0].astype(np.complex128) / 255
    img = (2 - img)
    img[img < 1.05] = 1 + 0j
    img[img > 1] = img[img > 1] / np.max(img) * ni
    img[img < 1] = 1
    img = np.expand_dims(img, 0)
    np.save(sample_npy, img)  # Real sample
    img[:, :, :] = 1
    np.save(ref_npy, img)  # Reference

    # ----------------------------------Solve the reference first-------------------------------------
    subprocess.run([scattervol_path+"scattervolume", "--sample", ref_npy, "--size", str(size[0]), str(size[1]), str(size[2]),
                    "--coef", str(1), str(1), "--lambda", str(lambdai), "--output", ref_cw], shell=False, capture_output=True)

    subprocess.run([scattervol_path+"scatterview", "--input", ref_cw, "--size", str(size[0]),
                   "--nogui", "--resolution", str(resolution), "--output", out_ref_npy,  "--axis", str(save_axis),
                    "--center", str(size[0]/2), str(size[0]/2), str(0), "--slice", str(relative_slice)], shell=False, capture_output=True)

    xy = np.load(out_ref_npy)
    intensity_ref = np.real(
        xy[:, :, 0] * np.conj(xy[:, :, 0]) + xy[:, :, 1] * np.conj(xy[:, :, 1]) + xy[:, :, 2] * np.conj(xy[:, :, 2]))
    ref = time.time()
    print("---------------Profiling---------------")
    print("Reference intensity takes " + str(ref-start) +"s.")


    # ----------------------------Calculate the intensity for real sample------------------------------
    subprocess.run([scattervol_path+"scattervolume", "--sample", sample_npy, "--size", str(size[0]), str(size[1]), str(size[2]),
                    "--coef", str(coefs[0]), str(coefs[1]), "--lambda", str(lambdai), "--output", result_cw], shell=False, capture_output=True)

    sample_volume = time.time()
    print("Sample scattervolume takes " + str(sample_volume-ref) +"s.")
    subprocess.run([scattervol_path+"scatterview", "--input", result_cw, "--size", str(size[0]),
                   "--nogui", "--resolution", str(resolution), "--output", out_result_npy,  "--axis", str(save_axis),
                    "--center", str(size[0]/2), str(size[0]/2), str(0), "--slice", str(relative_slice)], shell=False, capture_output=True)

    xy = np.load(out_result_npy)
    intensity_sample = np.real(
        xy[:, :, 0] * np.conj(xy[:, :, 0]) + xy[:, :, 1] * np.conj(xy[:, :, 1]) + xy[:, :, 2] * np.conj(xy[:, :, 2]))
    sample_view = time.time()
    print("Sample scatterview takes " + str(sample_view-sample_volume) +"s.")
    # plt.imshow(np.real(intensity_sample))
    # plt.colorbar()
    # plt.show()
    A = np.log(intensity_ref[2**(resolution-1), 2**(resolution-1)] / intensity_sample[2**(resolution-1), 2**(resolution-1)])
    # Let's run scatterview instead. Might take less time
    # # Read the .cw file and evaluate the intensity on the lower boundary
    # layer = coupledwave()
    # layer.load(result_cw)
    print("A = " + str(A)+"\n")
    recon_A.append(A)
    # plt.plot(arr[0:i, 0], recon_A)
    # plt.savefig(scattervol_path+"spectroscopy\\"+"A_for_" + str(lambdai+1) + "_points.png")
    np.save(data_path + str(lambdai)+"recon_A_odd_coef40_40.npy", recon_A)
np.save(source_path + "recon_A_odd_coef40_40.npy", recon_A)
#plt.plot(recon_A)
#plt.show()