extract rescale slope value

Question

mohd akmal masud 2018-2-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/383405-extract-rescale-slope-value

评论： mohd akmal masud 2019-3-24

Hi all,

I want to extract the RescaleSlope value from dicominfo for each slice. But I have 135 slice images. This is my code to extract RescaleSlope. But still failed.

P = zeros(256, 256, 135);
for K = 1 : 135
    petname = sprintf('PET_I1001_PT%03d.dcm', K);
    P(:,:,K) = dicominfo(petname);
end
info=dicominfo(P(:,:,K));
[r,c,slice] = findND (info.RescaleSlope);

Anyone who can help me solve this problem???

2 个评论
显示无隐藏无

Image Analyst 2018-2-19

I fixed your formatting for you. http://www.mathworks.com/matlabcentral/answers/13205#answer_18099

mohd akmal masud 2018-2-19

Thank you Image Analyst

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2018-2-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/383405-extract-rescale-slope-value#answer_305807

在 MATLAB Online 中打开

for K = 135 : -1 : 1
  petname = sprintf('PET_I1001_PT%03d.dcm', K);
  info(K) = dicominfo(petname);
end
rescale_slopes = [info.RescaleSlope];

If, for some reason you need to find the non-zero entries, then

slice_idx = find(rescale_slopes);

Note: the above code has a limitation that the dicominfo returned by each slice must have exactly the same set of fields. If that assumption is violated then the info(K) assignment will give you an error about assignment between dissimilar structures.

14 个评论
显示 12更早的评论隐藏 12更早的评论

mohd akmal masud 2018-11-1

在 MATLAB Online 中打开

Hi Walter Roberson,

i have this code for extract the location x,y and z. [r,c,slice] = findND(P >= 20000 & P <= 32767). But how to me extract the RescaleSlope for it?

P = zeros(256, 256, 47);
  for K = 20 : 31
     petname = sprintf('PETWB001_PT%03d.dcm', K);
     P(:,:,K) = dicomread(petname);
  end
  numOfPixels = numel(P(:));
  max(P(:));
  sum(P(:));
    mask = (P >= 20000 & P <= 32767);
numPixelsInRange = sum(mask(:));
 sum(P(P>=20000 & P<=32767));
 [r,c,slice] = findND(P >= 20000 & P <= 32767)
 val=P(P(:)>=20000 & P(:)<=32767 );
    if true
      % code
    end

because i use your code below, it extract all the RescaleSlope from slice 1 till 47.

for K = 135 : -1 : 1
  petname = sprintf('PET_I1001_PT%03d.dcm', K);
  thisinfo = dicominfo(petname);
  fn = fieldnames(thisinfo);
  for N = 1 : length(fn)
    thisfield = fn{N};
    info(K).(thisfield) = thisinfo.(thisfield);
  end
  this_slice = dicomread(petname);
  slices{K} = this_slice;
end
rescale_slopes = {info.RescaleSlope};
mask = cellfun(@isempty, rescale_slopes);
rescale_slopes(mask) = {0};       %set missing entries to 0
rescale_slopes = cell2mat(rescale_slopes);
inrange_mask = cellfun(@(M) any(M(:) >= 25000 & M(:) <= 32767), slices);
idx_of_inrange = find(inrange_mask);
rescale_slope_of_inrange = rescale_slopes(inrange_mask);

Walter Roberson 2018-11-2

在 MATLAB Online 中打开

If you have read in all of the values, then

rescale_slopes(20:31)

If for some reason you only want to read in a subset of the slices, then

slices_to_read = 20:31;
num_slices = length(slices_to_read);
for slice_idx = num_slices: -1 : 1
 slice_number = slices_to_read(slice_idx);
  petname = sprintf('PET_I1001_PT%03d.dcm', slice_number);
  info(slice_idx) = dicominfo(petname);
  if ~isfield(info(slice_idx), 'RescaleSlope') || isempty(info(slice_idx).RescaleSlope)
   info(slice_idx).RescaleSlope = 1;
  end
end

The looping from the last backwards towards the first is done for efficiency: it forces the last info() entry to be assigned to first, so afterwards it is not necessary to expand the info() struct the way it would be if you had not pre-allocated the info struct and you were looping forwards.

mohd akmal masud 2019-3-24

在 MATLAB Online 中打开

Sorry all, another question i have but i wrote at this space. please help me

Dear all,

this is my code to view CT image by slice

P = zeros(256, 256, 72);
for K = 1 : 72
    petname = sprintf('I4%03d.dcm', K);
    P(:,:,K) = dicomread(petname);
end
imshow3D(P)

then, this is my code for view SPECT image by slice,

Noted: all my 42 slice SPECT image stored in one file.

[spect map]=dicomread('128x128');
info = dicominfo('128x128');
gp=info.SliceThickness;
spect=(squeeze(spect));%smooth3
aa=size(spect);aa=aa(3);
imshow3D(spect);

Anybody can help me to fuse both SPECT and CT images for all slice?

请先登录，再进行评论。

Answer 2

mohd akmal masud 2018-11-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/383405-extract-rescale-slope-value#answer_345774

在 MATLAB Online 中打开

Hi Walter

I Have This Code

P = zeros(256, 256, 47);
  for K = 20 : 31
     petname = sprintf('PETWB001_PT%03d.dcm', K);
     P(:,:,K) = dicomread(petname);
  end
  numOfPixels = numel(P(:));
  max(P(:));
  sum(P(:));
    mask = (P >= 20000 & P <= 32767);
numPixelsInRange = sum(mask(:));
 sum(P(P>=20000 & P<=32767));
 [r,c,slice] = findND(P >= 20000 & P <= 32767);
 val=P(P(:)>=20000 & P(:)<=32767 )

how to multiply the val with RescaleSlope itself?

i have try the code you given like below. its work, but separately. So have to multiply manually. So how to me make automatically?

for K = 47 : -1 : 1
  petname = sprintf('PETWB001_PT%03d.dcm', K);
  thisinfo = dicominfo(petname);
  fn = fieldnames(thisinfo);
  for N = 1 : length(fn)
    thisfield = fn{N};
    info(K).(thisfield) = thisinfo.(thisfield);
  end
  this_slice = dicomread(petname);
  slices{K} = this_slice;
   if ~isfield(info(K), 'RescaleSlope') || isempty(info(K).RescaleSlope)
     info(K).RescaleSlope = 1;
    end
end
rescale_slopes = {info.RescaleSlope};
mask = cellfun(@isempty, rescale_slopes);
rescale_slopes(mask) = {0};       %set missing entries to 0
rescale_slopes = cell2mat(rescale_slopes);
inrange_mask = cellfun(@(M) any(M(:) >= 25000 & M(:) <= 32767), slices);
idx_of_inrange = find(inrange_mask);
rescale_slope_of_inrange = rescale_slopes(20:31)

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2019-3-8

在 MATLAB Online 中打开

projectdir = 'I:\12. FIZIK NUKLEAR\Jadual Bulanan Kerja PSF 2018\75ct';
numslices = 75;  %assume this
rawslices = cell(numslices, 1);
slices = cell(numslices, 1);
for K = numslices : -1 : 1
  petname = fullfile(projectdir, sprintf('ct%d.dcm', K));
  thisinfo = dicominfo(petname);
  fn = fieldnames(thisinfo);
  for N = 1 : length(fn)
    thisfield = fn{N};
    info(K).(thisfield) = thisinfo.(thisfield);
  end
  this_slice = dicomread(petname);
  if ~isfield(info(K), 'RescaleSlope') || isempty(info(K).RescaleSlope)
     info(K).RescaleSlope = 1;
  end
  if ~isfield(info(K), 'RescaleIntercept') || isempty(info(K).RescaleIntercept)
     info(K).RescaleIntercept = 0;
  end
  rawslices{K} = this_slice;
  slices{K} = double(this_slice) * info(K).RescaleSlope + info(K).RescaleIntercept;
end

Now rawslices is a cell array containing the data as read in from the file (probably int16 datatype), and slices is a cell array containing the data rescaled, as double datatype.

You do not need to do any of that post-processing to find empty values because it is already done inside the reading loop: the isempty() causes empty values to be replaced with known values.

Walter Roberson 2019-3-9

What extension does it have?

This is a Question about DICOM, dealing with the DICOM-specific matter of RescaleSlope and RescaleIntercept. If you are not using DICOM images then your question should have been posted separately and should have included detail of what you were trying to achieve.

mohd akmal masud 2019-3-24

在 MATLAB Online 中打开

Sorry all, another question i have but i wrote at this space. please help me

Dear all,

this is my code to view CT image by slice

P = zeros(256, 256, 72);
for K = 1 : 72
    petname = sprintf('I4%03d.dcm', K);
    P(:,:,K) = dicomread(petname);
end
imshow3D(P)

then, this is my code for view SPECT image by slice,

Noted: all my 42 slice SPECT image stored in one file.

[spect map]=dicomread('128x128');
info = dicominfo('128x128');
gp=info.SliceThickness;
spect=(squeeze(spect));%smooth3
aa=size(spect);aa=aa(3);
imshow3D(spect);