Skip to content

Commit 6334625

Browse files
nmelladonmellado
authored
Merge pull request #2 from FrancoisGaits/main
Scus final version
2 parents 1f98952 + 9bc925c commit 6334625

File tree

113 files changed

+64984
-2
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

113 files changed

+64984
-2
lines changed

‎CMakeLists.txt

+77
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
cmake_minimum_required(VERSION 3.8)
2+
project(Scus)
3+
4+
set(CMAKE_CXX_STANDARD 17)
5+
set(CMAKE_CXX_STANDARD_REQUIRED ON)
6+
7+
# search for DCMTK library and header files
8+
find_package(DCMTK REQUIRED)
9+
10+
add_compile_options(-Wunknown-pragmas -fopenmp)
11+
12+
# specify DCMTK header include directories
13+
include_directories(${DCMTK_INCLUDE_DIRS})
14+
15+
set(Files src/main.cpp
16+
src/utility/samplingUtils.hpp
17+
18+
src/samplers/SamplerDartThrowing.hpp
19+
20+
src/samplers/SamplerRegularGrid.hpp
21+
22+
src/samplers/SamplerCascadedSobol.h
23+
src/samplers/cascadedSobol/OwenScrambling.h
24+
src/samplers/cascadedSobol/Random.h
25+
src/samplers/cascadedSobol/SobolGenerator1D.h
26+
27+
src/pointsets/Domain.hpp
28+
src/pointsets/Histogram.hpp
29+
src/pointsets/Point.hpp
30+
src/pointsets/Pointset.hpp
31+
src/pointsets/Vector.hpp
32+
33+
src/grid/grid.cpp src/grid/grid.h
34+
src/grid/cell.cpp src/grid/cell.h
35+
36+
src/geometry/acquisitionZone.cpp src/geometry/acquisitionZone.h
37+
38+
src/volume/weightingFunction.h
39+
src/volume/sphere.h
40+
src/volume/plane.h
41+
src/volume/dicom.h src/volume/dicom.cpp
42+
43+
src/utility/json.hpp
44+
src/utility/stb_image.h
45+
src/utility/stb_image_write.h
46+
src/utility/L2Discrepancy.hpp
47+
src/utility/GeneralizedL2Discrepancy.hpp
48+
src/utility/Diaphony.hpp
49+
50+
src/utility/pointWriter.cpp src/utility/pointWriter.h
51+
src/utility/TestUtils.cpp src/utility/TestUtils.h
52+
src/volume/weightingFactory.h src/volume/cube.h
53+
src/samplers/SamplerUniformRandom.h src/samplers/Sampler.h
54+
src/samplers/SamplerFactory.h src/volume/dicom.cpp src/volume/dicom.h
55+
src/geometry/scatMap.cpp src/geometry/scatMap.h src/geometry/basicUS.cpp
56+
src/geometry/basicUS.h src/Scus.h
57+
src/geometry/PSF.h src/geometry/PSF.cpp
58+
59+
src/utility/utils.h
60+
src/geometry/realScatMap.cpp src/geometry/realScatMap.h
61+
src/utility/TestUtils.cpp src/utility/TestUtils.h
62+
63+
)
64+
65+
add_library (ScusLib ${Files} )
66+
67+
add_executable(${PROJECT_NAME} ${Files})
68+
69+
70+
# link DCMTK library files
71+
target_link_libraries(${PROJECT_NAME} ${DCMTK_LIBRARIES} tbb)
72+
target_link_libraries(ScusLib ${DCMTK_LIBRARIES} tbb)
73+
find_package(OpenMP)
74+
if(OpenMP_CXX_FOUND)
75+
target_link_libraries(${PROJECT_NAME} OpenMP::OpenMP_CXX)
76+
target_compile_definitions(${PROJECT_NAME} PUBLIC OPENMP=1)
77+
endif()

‎FieldII/GenerateImage.py

+64
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
import matlab.engine
2+
import shutil
3+
import os
4+
import subprocess
5+
import time
6+
7+
import visualisationV2
8+
9+
NB_INSTANCES = 5
10+
11+
# Copy points to folder
12+
file = "../cmake-build-release/raw_res.dat"
13+
14+
if os.path.exists(file) :
15+
shutil.copy(file, "raw_res.dat")
16+
else:
17+
print("Unable to locate :",file, "use point generation tool to create it")
18+
exit(-1)
19+
20+
# Clear previous rf (with backup)
21+
folder = "rf_data"
22+
if not os.path.exists(folder):
23+
os.mkdir(folder)
24+
else :
25+
if os.path.exists(folder + "_old"):
26+
shutil.rmtree(folder + "_old")
27+
os.rename(folder, folder + "_old")
28+
os.mkdir(folder)
29+
30+
for filename in os.listdir(folder):
31+
file_path = os.path.join(folder, filename)
32+
try:
33+
if os.path.isfile(file_path) or os.path.islink(file_path):
34+
os.unlink(file_path)
35+
elif os.path.isdir(file_path):
36+
shutil.rmtree(file_path)
37+
except Exception as e:
38+
print("Failed to delete %s. Reason: %s" % (file_path, e))
39+
40+
# Create first matlab instance
41+
# Translate .dat to .mat
42+
eng = matlab.engine.start_matlab()
43+
eng.mk_pht(nargout=0)
44+
print("Phantom updated")
45+
46+
# Creating more instances
47+
processes = []
48+
for i in range(NB_INSTANCES):
49+
processes.append(subprocess.Popen(["Python", "ParallelGenerate.py"]))
50+
51+
print("Process launched")
52+
print("Waiting for processes...")
53+
54+
#visualisationV2.plot_file("../cmake-build-release/res.dat")
55+
56+
for proc in processes:
57+
proc.wait()
58+
print("A process finished")
59+
60+
eng.make_image(nargout=0)
61+
62+
63+
input("End simulation")
64+
eng.quit()

‎FieldII/Mat_field.mexw64

155 KB
Binary file not shown.

‎FieldII/ParallelGenerate.py

+9
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
import matlab.engine
2+
3+
print("Launching Matlab...")
4+
eng = matlab.engine.start_matlab()
5+
6+
print("Matlab launched")
7+
eng.field_init(0)
8+
9+
eng.sim_img(nargout=0)

‎FieldII/calc_scat.m

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
% Procedure for calculating the received signal from a collection of scatterers.
2+
%
3+
% Calling: [scat, start_time] = calc_scat(Th1, Th2, points, amplitudes);
4+
%
5+
% Parameters: Th1 - Pointer to the transmit aperture.
6+
% Th2 - Pointer to the receive aperture.
7+
% points - Scatterers. Vector with three columns (x,y,z)
8+
% and one row for each scatterer.
9+
% amplitudes - Scattering amplitudes. Row vector with one
10+
% entry for each scatterer.
11+
%
12+
% Return: scat - Received voltage trace.
13+
% start_time - The time for the first sample in scat.
14+
%
15+
% Version 1.0, November 28, 1995 by Joergen Arendt Jensen
16+
17+
function [scat, start_time] = calc_scat (Th1, Th2, points, amplitudes)
18+
19+
% Check the point array
20+
21+
[m1,n]=size(points);
22+
if (n ~= 3)
23+
error ('Points array must have three columns');
24+
end
25+
[m2,n]=size(amplitudes);
26+
if (m1 ~= m2)
27+
error ('There must be the same number of rows for points and amplitudes arrays');
28+
end
29+
30+
% Call the C-part of the program to show aperture
31+
32+
[scat, start_time] = Mat_field (4005,Th1,Th2,points, amplitudes);
33+
34+

‎FieldII/field.m

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
% Start up the Field II simulation system
2+
%
3+
% Version 1.1, August 13, 2007, JAJ
4+
5+
% Replace the second path with the name of the directory
6+
% for your Field II m-files and executable
7+
8+
path(path, '/home/jaj/programs/field_II/M_files')
9+
10+
field_init

‎FieldII/field_init.m

+30
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
% Procedure for initializing the Field II program system. Must be
2+
% the first routine that is called before using the system.
3+
%
4+
% Calling: field_init (suppress);
5+
%
6+
% Return: nothing.
7+
%
8+
% Input: suppress: An optional argument suppress with a value
9+
% of zero can be given to suppress the
10+
% display of the initial field screen.
11+
% No ACII ouput will be given, if the argument
12+
% is -1. Debug messages will be written if
13+
% enable by field_debug, and all error messages
14+
% will also be printed.
15+
%
16+
% Version 1.2, January 20, 1999 by Joergen Arendt Jensen
17+
18+
function res = field_init (suppress)
19+
20+
% Call the C-part of the program to initialize it
21+
22+
res = 1;
23+
24+
if (nargin==1)
25+
Mat_field (5001,suppress);
26+
else
27+
Mat_field (5001,1);
28+
end
29+
30+

‎FieldII/field_logo.m

+21
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
% Function to display the logo for field
2+
%
3+
% Version 1.3, August 10, 2007 by Joergen Arendt Jensen
4+
% Error in loading filr fixed
5+
6+
function res = field_logo
7+
8+
% Create a window and display the Field II logo
9+
10+
h=figure;
11+
axes('position',[0 0 1 1]);
12+
place=which ('logo_field.mat');
13+
eval(['load ',place])
14+
image(data1);
15+
axis off
16+
colormap(map);
17+
drawnow;
18+
pause(5)
19+
close(h);
20+
21+

‎FieldII/logo_field.mat

279 KB
Binary file not shown.

‎FieldII/make_image.m

+68
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
% Compress the data to show 60 dB of
2+
% dynamic range for the cyst phantom image
3+
%
4+
% version 1.3 by Joergen Arendt Jensen, April 1, 1998.
5+
% version 1.4 by Joergen Arendt Jensen, August 13, 2007.
6+
% Calibrated 60 dB display made
7+
8+
f0=3.0e6; % Transducer center frequency [Hz]
9+
fs=4*f0; % Sampling frequency [Hz]
10+
c=1540; % Speed of sound [m/s]
11+
no_lines=161; % Number of lines in image
12+
image_width=50/1000; % Size of image sector
13+
d_x=image_width/no_lines; % Increment for image
14+
15+
% Read the data and adjust it in time
16+
17+
18+
disp('Finding the envelope')
19+
min_sample=0;
20+
for i=1:no_lines
21+
22+
% Load the result
23+
24+
cmd=['load rf_data/rf_ln',num2str(i),'.mat'];
25+
disp(cmd)
26+
eval(cmd)
27+
28+
% Find the envelope
29+
30+
rf_env=abs(hilbert([zeros(round(tstart*fs-min_sample),1); rf_data]));
31+
env( 1:max(size(rf_env)) ,i)=rf_env;
32+
end
33+
34+
save("envtst.mat","env");
35+
36+
% Do logarithmic compression
37+
38+
D=3; % Sampling frequency decimation factor
39+
40+
disp('Log Compression')
41+
log_env=env(1:D:max(size(env)),:)/max(max(env));
42+
log_env=20*log10(log_env);
43+
log_env=127/60*(log_env+60);
44+
45+
% Make an interpolated image
46+
47+
disp('Doing interpolation')
48+
ID=20;
49+
[n,m]=size(log_env);
50+
new_env=zeros(n,m*ID);
51+
52+
for i=1:n
53+
new_env(i,:)=abs(interp(log_env(i,:),ID));
54+
end
55+
56+
[n,m]=size(new_env);
57+
58+
fn=fs/D;
59+
clf
60+
a=((1:(ID*no_lines-1))*d_x/ID-no_lines*d_x/2)*1000;
61+
b=((1:n)/fn+min_sample/fs)*c/2*1000;
62+
image(a,b,new_env)
63+
xlabel('Lateral distance [mm]')
64+
ylabel('Axial distance [mm]')
65+
colormap(gray(127))
66+
axis('image')
67+
axis([-image_width*500 image_width*500 35 90])
68+

‎FieldII/mk_pht.m

+9
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
% Make the scatteres for a simulation and store
2+
% it in a file for later simulation use
3+
4+
% Joergen Arendt Jensen, Feb. 26, 1997
5+
6+
function mk_pht()
7+
[phantom_positions, phantom_amplitudes] = scus_phantom();
8+
9+
save pht_data.mat phantom_positions phantom_amplitudes

‎FieldII/saveEnv.m

+36
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
MSEs = [];
2+
NMSEs = [];
3+
4+
X = [1 2 3 4 5 6 7 8 9 10 15 20 27 64 125 343];
5+
6+
for i = X
7+
load("enveloppes/uniform/" + int2str(i));
8+
9+
envcrop = env(1100:3000,2:48);
10+
envcropf = reshape(envcrop.',1,[]);
11+
envcropf = (envcropf-min(envcropf))/(max(envcropf)-min(envcropf));
12+
13+
h = histfit(envcropf,64,"rayleigh");
14+
histdata = get(h(1),'YData');
15+
redcurve = get(h(2),'YData');
16+
17+
xgrid64 = linspace(0,1,64);
18+
xgrid100 = linspace(0,1,100);
19+
20+
21+
%redcurve2 = redcurve;
22+
%redcurve2 = interp1(xgrid100,redcurve,xgrid64);
23+
redcurve2 = [mean(reshape(redcurve, 2, [])) zeros(1,14)]; % Mean Of Bin Edges Fit
24+
fit_err = histdata - redcurve2; % Error
25+
26+
mse=norm(redcurve2-histdata,2)^2/length(redcurve2);
27+
sigEner=norm(redcurve2)^2;
28+
nmse=(mse/sigEner);
29+
30+
NMSEs = [NMSEs nmse];
31+
MSEs = [MSEs mse];
32+
33+
end
34+
35+
36+

‎FieldII/scus_phantom.m

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
function [positions, amp] = scus_phantom()
2+
load raw_res.dat
3+
4+
z_start = 30/1000;
5+
6+
x = raw_res(:,1)/1000;
7+
y = raw_res(:,2)/1000;
8+
z = raw_res(:,3)/1000 + z_start;
9+
amp = raw_res(:,4);
10+
11+
positions = [x y z];

0 commit comments

Comments
 (0)