parfor with temporary variable worning

I need expert's help.

I attached a code below:

The worning signal indicated "temp_info".

I want to stack CT images as a 3D volume.

But, position information of each slice is on the temp_info (loaded by dicominfo fn.)

I will handle few hundreds patient data.

Seriously, i need your help...

Best regards

Wonjoong Cheon

function [outputArg1] = dicomread3D(path_)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

ct_list = dir(fullfile(path_, '*.dcm'));
num = length(ct_list);
temp_img_size = dicomread(fullfile(ct_list(1).folder, ct_list(1).name));
dcm_stack = zeros(temp_img_size(1), temp_img_size(2), num);

f = waitbar(0,'Please wait...');
parfor iter2 = 1: num
    temp_img = dicomread(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    temp_info = dicominfo(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    ct_location = temp_info.InstanceNumber;
    dcm_stack(:,:,ct_location) = temp_img;
    waitbar(iter2/num, f, 'Load 3D medical images...');
end
close(f)

try
    outputArg1 = (dcm_stack.*temp_info.RescaleSlope)+temp_info.RescaleIntercept;
catch
    outputArg1 = dcm_stack;
    disp('Rescale information is empty')
end

end

Answer 1

[Edric Ellis]

I think there are a couple of problems here. Firstly, I think you need to change the way you're assigning into dcm_stack because as written, you are not satisfying the "sliced variable" requirement for a parfor loop output. I think you can fix this by storing an additional sliced output containing the the ct_location, I think you should apply the slope and intercept inside the parfor loop. Here's a sketch:

ct_list = dir(fullfile(path_, '*.dcm'));
num = length(ct_list);
temp_img_size = dicomread(fullfile(ct_list(1).folder, ct_list(1).name));
dcm_stack = zeros(temp_img_size(1), temp_img_size(2), num);

% Additional output for the parfor loop:
ct_locations = zeros(1, num);

parfor iter2 = 1: num
    temp_img = dicomread(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    temp_info = dicominfo(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    % Store information from temp_info
    ct_locations(iter2) = temp_info.InstanceNumber;
    % Note we must store in stack indexed using "iter2":
    dcm_stack(:,:,iter2) = (temp_img .* temp_info.RescaleSlope) + temp_info.RescaleIntercept;
end

% You can now re-order the "pages" in dcm_stack like this:
dcm_stack = dcm_stack(:,:,ct_locations);

[Walter Roberson]

Each different file could potentially have its own slope and rescale (or be missing them) so I do recommend doing the transform inside the loop. Maybe

if ~isfield(temp_info, 'RescaleSlope'); temp_info.RescaleSlope = 1; end
if ~isfield(temp_info, 'RescaleIntercept'); temp_info.RescaleIntercept = 0; end

just before the assignment to dcm_stack(:,:,iter2)

[wonjoong cheon]

Dear Edric Ellis

I agree with your idea.

But, "dcm_stack = dcm_stack(:,:,ct_locations);"

is not working...

wht..happen?

Best regards

[Walter Roberson]

I need to block access from other threads such as gate keeper

If you are using parfor() then the only other thread that exists is the client, unless you used java to create threads, or you used system() or System.Diagnostics.Process to create additional processes

The parfor communications mechanisms ensure that different workers do not overwrite each other when setting variables that are indexed by the parfor index variable.

[wonjoong cheon]

It is not perfect

But It is working: I don't like that for statement for re-oredering

Thank you for your help.

function [outputArg1] = dicomread3D(path_)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

ct_list = dir(fullfile(path_, '*.dcm'));
num = length(ct_list);
temp_img_size = size(dicomread(fullfile(ct_list(1).folder, ct_list(1).name)));
temp_infomation = dicominfo(fullfile(ct_list(1).folder, ct_list(1).name));
dcm_stack = zeros(temp_img_size(1), temp_img_size(2), num);

sliceInfo = zeros(1, num);
WaitMessage = parfor_wait(num, 'Waitbar', true);
parfor iter2 = 1: num
    temp_img = dicomread(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    temp_info = dicominfo(fullfile(ct_list(iter2).folder, ct_list(iter2).name));

    sliceInfo(iter2) = temp_info.InstanceNumber;
    dcm_stack(:,:,iter2) = temp_img;
    
    WaitMessage.Send;
    labBarrier
end
WaitMessage.Destroy

% Re-ordering 
for iter1= 1: num
    [~, c] = find(iter1== sliceInfo);
    dcm_stackN(:,:,iter1) = dcm_stack (:,:,c);
end

try
    outputArg1 = (dcm_stackN.*temp_infomation .RescaleSlope)+temp_infomation .RescaleIntercept;
catch
    outputArg1 = dcm_stackN;
    disp('Rescale information is empty')
end

end

[Walter Roberson]

Why isn't the reodering

[~, sortorder] = sort(sliceInfo);
dcm_stackN = dcm_stack(:,:,sortorder);

[Walter Roberson]

outputArg1 = (dcm_stackN.*temp_infomation .RescaleSlope)+temp_infomation .RescaleIntercept;

No, every file could have a different rescale slope and intercept (or not have the information.) You should be rescaling inside the parfor.

    temp_img = dicomread(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    temp_info = dicominfo(fullfile(ct_list(iter2).folder, ct_list(iter2).name));
    sliceInfo(iter2) = temp_info.InstanceNumber;
    
    if ~isfield(temp_info, 'RescaleSlope'); temp_info.RescaleSlope = 1; end
    if ~isfield(temp_info, 'RescaleIntercept'); temp_info.RescaleIntercept = 0; end
    
    dcm_stack(:,:,iter2) = (temp_img .* temp_info.RescaleSlope) + temp_info.RescaleIntercept;

[wonjoong cheon]

Dear Walter Roberson

IT IS DONE.

Great skill from the level 10 MVP.

Best regards.

Answer 2

[Walter Roberson]

temp_info is assigned to inside the parfor in a way that is not indexed by the loop variable. parfor automatically makes it a temporary variable, and that means it cannot be used after the loop.

[wonjoong cheon]

Do you have any tips for this situnation ?