Skip to content

Commit bf1b067

Browse files
author
brandon.nelson
committed
Merge branch 'main' of https://github.com/DIDSR/LCD_CT into main
2 parents 436d87d + 81c8d2e commit bf1b067

File tree

5 files changed

+415
-1
lines changed

5 files changed

+415
-1
lines changed

‎README.rst‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ Low Contrast Detectability for CT Toolbox
99
:width: 800
1010
:align: center
1111

12-
- **Regulatory Science Tool:** Check the FDA website for a description of the LCD-CT toolbox in the `Regulatory Science Tool Catalog <https://www.fda.gov/medical-devices/science-and-research-medical-devices/lcd-ct-low-contrast-detectability-lcd-test-assessing-advanced-nonlinear-ct-image-reconstruction-and>`_
12+
- **Regulatory Science Tool:** Check the FDA website for a description of the LCD-CT toolbox in the `Regulatory Science Tool Catalog <https://cdrh-rst.fda.gov/lcd-ct-low-contrast-detectability-lcd-test-assessing-advanced-nonlinear-ct-image-reconstruction-and>`_
1313

1414
.. |zenodo| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.7996580.svg
1515
:alt: Zenodo Data Access

‎additional_demos/demo_FDK_Water_Phantom_PCD‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
% Inputs
2+
% signal present 3D
3+
% signal absent 3D
4+
% ground truth 2D
5+
6+
close all
7+
clear all
8+
clc
9+
restoredefaultpath
10+
11+
mfilename = '/home/sriharsha.marupudi/TIGRE-master/MATLAB';
12+
13+
load_coords_filename = '/LCD_CT/data/coordinates/coordinates_20241008.mat';
14+
15+
16+
%% Load your data
17+
load('data_Iodine_Water_Validation_3_20241015_Recon_Binned.mat', 'img_low_binned');
18+
img = (img_low_binned(:,:,200:210));
19+
20+
% Convert to HU
21+
% P_water_actual = 75;
22+
% P_air_actual = -10;
23+
%
24+
% HU_water = 0;
25+
% HU_air = -1000;
26+
%
27+
% m = (HU_air - HU_water) / (P_air_actual - P_water_actual);
28+
% b = HU_water - m * P_water_actual;
29+
%
30+
% convertToHU = @(P) m * P + b;
31+
%
32+
% img = convertToHU(img);
33+
34+
figure; imagesc(img(:,:,5)); colormap gray; axis off; axis tight; axis equal;
35+
36+
%% Mask Image to remove container
37+
[rows, cols, slices] = size(img);
38+
centerX = round(cols / 2);
39+
centerY = round(rows / 2);
40+
radius = min(rows, cols) / 2;
41+
42+
[X, Y] = meshgrid(1:cols, 1:rows);
43+
circularMask = (X - centerX).^2 + (Y - centerY).^2 <= radius^2;
44+
45+
maskedImage = zeros(size(img));
46+
47+
for k = 1:slices
48+
currentSlice = img(:,:,k);
49+
maskedSlice = zeros(size(currentSlice));
50+
maskedSlice(circularMask) = currentSlice(circularMask);
51+
maskedImage(:,:,k) = maskedSlice;
52+
end
53+
54+
figure;
55+
imagesc(maskedImage(:,:,5));
56+
colormap gray; axis off; axis tight; axis equal;
57+
title('Masked Image');
58+
59+
%% Threshold masked image to get inserts
60+
th = 100; %if not HU
61+
% th = 300; % if HU Threshold value
62+
diskSize = 5; % Disk size for morphological operations
63+
areaThreshold = 50; % Minimum area for keeping objects in mask
64+
65+
maskedImage_ground_truth = mean(maskedImage(:,:,1:10), 3);
66+
binaryMask = maskedImage_ground_truth > th;
67+
68+
% Remove small objects from the binary mask
69+
cleanedMask = bwareaopen(binaryMask, areaThreshold);
70+
71+
% Apply morphological closing to smooth the mask
72+
se = strel('disk', diskSize);
73+
cleanedMask = imclose(cleanedMask, se);
74+
75+
cleanedGrayscaleImage = zeros(size(maskedImage_ground_truth));
76+
cleanedGrayscaleImage(cleanedMask) = img(cleanedMask);
77+
78+
figure;
79+
imshow(cleanedGrayscaleImage, []);
80+
colormap(gray);
81+
title('Cleaned Grayscale Image');
82+
83+
ground_truth = cleanedGrayscaleImage;
84+
85+
%% Signal Free Image
86+
% Apply insert locations to original image and set inserts to 0
87+
88+
signalFreeMask = ~cleanedMask;
89+
signalFreeImage = zeros(size(img));
90+
91+
for k = 1:slices
92+
currentSlice = maskedImage(:,:,k);
93+
maskedSlice = zeros(size(currentSlice));
94+
maskedSlice(signalFreeMask) = currentSlice(signalFreeMask);
95+
signalFreeImage(:,:,k) = maskedSlice;
96+
end
97+
98+
figure;
99+
imagesc(signalFreeImage(:,:,5));
100+
colormap(gray); axis equal; axis tight; axis off;
101+
title('Signal Free Image');
102+
103+
104+
%% Specify observers to use
105+
observers = {LG_CHO_2D()};
106+
% observers = {LG_CHO_2D(), DOG_CHO_2D(), GABOR_CHO_2D(), ... };
107+
108+
%% Specify base directory and run the measurement
109+
base_directory = '/your/base/directory/path'; % Adjust to your actual base directory
110+
111+
offset = 0;
112+
n_reader = 10;
113+
n_inserts = 6;
114+
insert_r = 30;
115+
res_table = measure_LCD1(maskedImage, signalFreeImage, observers, ground_truth, offset,n_reader, n_inserts,insert_r,load_coords_filename);
116+
117+
118+
%% Plot and summarize results
119+
% Plot results
120+
custom_insert_names = {'5 Rod', '5 Solution','7.5 Rod', '7.5 Solution', '10 Rod', '10 Solution'}; % mg/mL of iodine inserts
121+
set_ylim = [0 1.2];
122+
plot_results1(res_table, [], custom_insert_names);
123+
124+
% Display results
125+
res_table
126+
127+
% Summarize results
128+
if ~is_octave
129+
groupsummary(res_table, ["observer", "insert_HU", "dose_level"],["mean", "std"])
130+
end
131+
%% or define a custom summary table by printing mean and standard deviation results
132+
nreader = max(res_table.reader);
133+
for i=1:6
134+
mean_AUC(i) = mean(res_table.auc([1:nreader]+(i-1)*nreader));
135+
std_AUC(i) = std(res_table.auc([1:nreader]+(i-1)*nreader));
136+
mean_snr(i) = mean(res_table.snr([1:nreader]+(i-1)*nreader));
137+
std_snr(i) = std(res_table.snr([1:nreader]+(i-1)*nreader));
138+
end
139+
140+
insert_HU = res_table.insert_HU(1:nreader:end);
141+
mean_AUC = mean_AUC(:);
142+
std_AUC = std_AUC(:)
143+
mean_snr = mean_snr(:);
144+
std_snr = std_snr(:)
145+
AUC_res = table(insert_HU, mean_AUC, std_AUC, mean_snr, std_snr);
146+
AUC_res
147+
148+

‎src/LCD/get_ROI_from_manual_selection.m‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
function roi = get_ROI_from_manual_selection(img, x_center, y_center, r)
2+
% Extract a square region from the image centered on (x_center, y_center) with radius r
3+
4+
[rows, cols,slices] = size(img);
5+
6+
% Define the bounding box for the region
7+
x_min = max(1, x_center - r);
8+
x_max = min(cols, x_center + r);
9+
y_min = max(1, y_center - r);
10+
y_max = min(rows, y_center + r);
11+
12+
% Extract the region of interest
13+
roi = img(y_min:y_max, x_min:x_max,:);
14+
end

‎src/LCD/measure_LCD1.m‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
function res_table = measure_LCD_single_dose(signal_present_img, signal_absent_img, observers, ground_truth, offset, n_reader, n_inserts, insert_r, load_coords_filename)
2+
% Calculate low contrast detectability using signal-present and signal-absent images.
3+
%
4+
% :param signal_present_img: The signal-present image
5+
% :param signal_absent_img: The signal-absent image (signal-free)
6+
% :param observers: List of observer objects or strings of the name of observers
7+
% :param ground_truth: Image or filename of the image with no noise of MITA LCD phantom
8+
% :param offset: An optional offset subtracted from each image (default 0)
9+
% :param n_reader: Number of readers (default is 10)
10+
% :param n_inserts: Number of inserts to be selected manually
11+
% :param insert_r: The radius of the inserts (same for all inserts)
12+
% :param load_coords_filename: Optional filename to load coordinates from (if it exists)
13+
%
14+
% :return res_table: Table ready for saving or plotting
15+
16+
%% Define observers
17+
if ~exist('observers', 'var')
18+
observers = {LG_CHO_2D()}; % Default observer
19+
end
20+
pct_split = 0.5;
21+
seed_split = randi(10000, n_reader, 1);
22+
23+
%% Load or select coordinates
24+
if exist(load_coords_filename, 'file')
25+
% Load previously saved coordinates
26+
coords_data = load(load_coords_filename);
27+
x_coords = coords_data.x_coords;
28+
y_coords = coords_data.y_coords;
29+
disp('Loaded coordinates from file.');
30+
else
31+
% Manually select the coordinates for each insert using ginput
32+
centerSliceIndex = round(size(signal_present_img, 3) / 2);
33+
figure; imagesc(signal_present_img(:, :, centerSliceIndex)); colormap gray; axis off; axis tight; axis equal;
34+
title('Select centers of the inserts');
35+
disp(['Please select ', num2str(n_inserts), ' insert locations']);
36+
[x_coords, y_coords] = ginput(n_inserts); % Get n_insert points manually
37+
38+
% Automatically save the coordinates to a file
39+
coords_data.x_coords = x_coords;
40+
coords_data.y_coords = y_coords;
41+
save(load_coords_filename, '-struct', 'coords_data');
42+
end
43+
44+
% Use the same radius for all inserts, defined by the input parameter `insert_r`
45+
insert_rs = ones(1, n_inserts) * insert_r; % Set the same radius for all inserts
46+
47+
%% Subtract offset from signal-present and signal-absent images
48+
sp_raw_array = signal_present_img - offset;
49+
sa_raw_array = signal_absent_img - offset;
50+
51+
% Initialize arrays for storing results
52+
observer = [];
53+
dose_level = [];
54+
snr = [];
55+
auc = [];
56+
reader = [];
57+
insert_HU = [];
58+
insert_diameter_pix = [];
59+
60+
%% Iterate through inserts and observers
61+
for i = 1:length(observers)
62+
model_observer = observers{i};
63+
64+
for insert_idx = 1:n_inserts
65+
% Get the manually selected insert coordinates and radius
66+
x_center = round(x_coords(insert_idx));
67+
y_center = round(y_coords(insert_idx));
68+
69+
% Update observer properties (if necessary)
70+
if isa(model_observer, "LG_CHO_2D")
71+
model_observer.channel_width = 2 / 3 * insert_r; % Update channel width for LG_CHO observer
72+
end
73+
74+
%% Extract regions of interest (ROIs) using the manually selected center and radius
75+
sp_imgs = get_ROI_from_manual_selection(sp_raw_array, x_center, y_center, insert_r);
76+
sa_imgs = get_ROI_from_manual_selection(sa_raw_array, x_center, y_center, insert_r);
77+
78+
% Remove DC component (mean intensity) from images
79+
for i_sp = 1:size(sp_imgs, 3)
80+
sp_imgs(:, :, i_sp) = sp_imgs(:, :, i_sp) - mean(sp_imgs(:, :, i_sp), 'all');
81+
end
82+
for i_sa = 1:size(sa_imgs, 3)
83+
sa_imgs(:, :, i_sa) = sa_imgs(:, :, i_sa) - mean(sa_imgs(:, :, i_sa), 'all');
84+
end
85+
86+
%% Perform observer study
87+
for n = 1:n_reader
88+
% Split the data into training and testing sets
89+
[sa_train, sa_test, sp_train, sp_test] = train_test_split(sa_imgs, sp_imgs, pct_split, seed_split(n));
90+
91+
% Perform the observer study
92+
res = model_observer.perform_study(sa_train, sp_train, sa_test, sp_test);
93+
94+
% Store the results
95+
observer = [observer; string(model_observer.type)];
96+
insert_HU = [insert_HU; mode(ground_truth(y_center, x_center))]; % Use the ground truth at the selected point
97+
insert_diameter_pix = [insert_diameter_pix; 2 * insert_r];
98+
dose_level = [dose_level; 1]; % Since you have only one dose level, set it to 1
99+
snr = [snr; res.snr];
100+
auc = [auc; res.auc];
101+
reader = [reader; n];
102+
end
103+
end
104+
end
105+
106+
% Return results as a table
107+
res_table = table(observer, insert_HU, dose_level, snr, auc, reader, insert_diameter_pix);
108+
end

0 commit comments

Comments
 (0)