extract rescale slope value

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???

[Image Analyst]

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

[mohd akmal masud]

Thank you Image Analyst

Answer 1

[Walter Roberson]

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.

[mohd akmal masud]

Thank You Walter Roberson. I will try

[mohd akmal masud]

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

rescale_slopes = [info.RescaleSlope];

I have try that above, but the error appear like below

   Subscripted assignment between dissimilar structures
  Error in imagescale(line 4)
  info(K)=dicominfo(petname);

[Walter Roberson]

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
end

rescale_slopes = {info.RescaleSlope};
mask = cellfun(@isempty, rescale_slopes);
rescale_slopes(mask) = {0};       %set missing entries to 0
rescale_slopes = cell2mat(rescale_slopes);
slice_idx = find(rescale_slopes);   %find non-zero entries, which would exclude the missing entries and any that happened to be 0 anyhow

[mohd akmal masud]

wowww amazing!!! its work. thank you walter roberson.

[mohd akmal masud]

if i just want info rescale slope that image have value pixel intensity from K >= 25000 & K <= 32767, is it true i write the command below:??

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 
end 
rescale_slopes = {info.RescaleSlope } ;
rescale_slopes(mask) = {K >= 25000 & K <= 32767}

[Walter Roberson]

To clarify: you want to know the rescale slope for images that have at least one pixel in the range 25000 to 32767 ?

[mohd akmal masud]

yes correct

[Walter Roberson]

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);

[mohd akmal masud]

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);

[mohd akmal masud]

But slice 1 till 19 not have the pixel value from 20000 to 32767..

[Walter Roberson]

After the line

    slices{K} = this_slice;

please add

    if ~isfield(info(K), 'RescaleSlope') || isempty(info(K).RescaleSlope)
     info(K).RescaleSlope = 1;
    end

[mohd akmal masud]

if i just want rescale slope from slice 20-31, is it correct just change below

for K = 135 : 20 : 31

[Walter Roberson]

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]

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]

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)

[mohd akmal masud]

Hi Walter,

I i have 75 ct slice(the hame file start with ct10, ct 20, ct30...........ct75), how could i write the code for looping?

75 ct slice inside this folder for example ( I:\12. FIZIK NUKLEAR\Jadual Bulanan Kerja PSF 2018\75ct)

[Walter Roberson]

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.

[mohd akmal masud]

But my image does not have extension name (.dcm)

[Walter Roberson]

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]

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?