Brightness issue with SART DBT reconstruction

[roshtha]

Expected Behavior

The reconstructed images should have comparable brighness level.

Actual Behavior

I run d19_DBTexample.m for reconstructing images from DBT breast phantom raw data. As I scroll through the reconstructed slices, the images get darker. Also the reconstructed images look blurry and lack details when compared with image outputs from actual DBT system.

Code to reproduce the problem (If applicable)

The geometry definition follows.

% Distances
geo.DSD = 650;  % Distance Source Detector(mm)
geo.DSO = 590;  % Distance Source Origin(mm)

% Detector parameters
geo.nDetector= [3584;2816];% number of pixels(px)
geo.dDetector= [0.083; 0.083]; % size of each pixel(mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector(mm)
% Image parameters
geo.nVoxel= [171;1127;803]; % number of voxels(vx)
geo.dVoxel = [0.25; 0.25; 0.25]; % total size of the image(mm)
geo.sVoxel = geo.nVoxel .* geo.dVoxel; % size of each voxel
% Offsets
Airgap = 15;    %DBT airgap (mm) 
geo.offOrigin = [ 0;0;0]; % Offset of image from origin (mm)
geo.offDetector=[0; 0]; 

geo.accuracy=0.5;
geo.COR=0;

nprojs = 16; % Number of projections
tubeangle = 15; % Angle range
angles=deg2rad(linspace(-tubeangle/2,tubeangle/2,nprojs)); % Angles in rad
geo.rotDetector=[0;0;0];

geo.mode='cone'; % Or 'parallel'. Geometry type.

%% Adapt CT geo to DBT
geo=staticDetectorGeo(geo,angles);

I set

geo.offOrigin = [0;0;0];

as I get a full black image wtih

geo.offOrigin =[((geo.sVoxel(1)/2)-(geo.DSD-geo.DSO)+Airgap);0;geo.sVoxel(3)/2];

Code for loading data. I refer this link https://github.com/CERN/TIGRE/issues/99 and modify input matrix data as

data = load('data1.mat');
img=data.img;
projections = -log(img./single(2^14-1)); %your maximum white level.
projections = rot90(permute(projections,[1 2 3]),1);

Calls SART algorithm as

recSART = SART(projections,geo,angles,2,'OrderStrategy','ordered');

The code for writing images to DICOM

a=permute(recSART,[3 2 1]);
a = imrotate(a,-90);
a = flipud(a);
filestring = '.\SART\';
for k=1:geo.nVoxel(1) 
    a(:,:,k) = a(:,:,k) + abs(min(a(:)));
    a(:,:,k) = (2^14-1).*mat2gray(a(:,:,k));
    a(:,:,k) = uint16(a(:,:,k));
    infoDicom.InstanceNumber = k;
    infoDicom.SOPClassUID = '1.2.840.10008.5.1.4.1.1.2';
    dicomwrite(uint16( a(:,:,k)),[filestring,int2str(k),'.dcm'],infoDicom(:,1), 'CreateMode', 'Copy');
end

Raw data matrix shared at https://drive.google.com/drive/folders/10er-5LklkIFZc3BR0Z_AKhpyFiaivTIC?usp=sharing

I also attach reconstructed image outputs from code and machine. code_output.zip machine_output.zip

Could anyone please provide inputs to understand the problem?

Specifications

• MATLAB version: Matlab 2015b
• OS: Windows 10, 64 bit
• CUDA version: 10.2

[AnderBiguri]

I will have a look at this. Could you also please share your particular geometry? DBT is a bit different to standard CT geometrically, and I am not 100% sure that the function int he demos is correct

[roshtha]

Hi, Thanks for the reply. I am sending geometry and specifications to tigre.toolbox@gmail.com. Please have a look.

[AnderBiguri]

Hi @roshtha ,

If you open staticDetectorGeo() there is an alternative piece of code commented. I think that is the correct way of doing static detector but a user told me that the other one was better. Can you change the code to the commented one and test it again, to see if that is the issue?

[roshtha]

Hi,

Thanks for the reply. I changed the code in staticDetectorGeo(). But problem still persists. From around slice 50, the image intensity diminishes and finally ends up in a dark image.

Plese note that I set geo.offOrigin = [0;0;0]; So value for Airgap is not considered in reconstructing images. Does that cause a problem?

I also attach d19_DBTexample.m modified for geometry parameters. d19_DBTexample.txt

Thanks.

[AnderBiguri]

@roshtha I see. Any reason why you have not added the airgap? I seems like a crucial parameter, its in the demo (coded by someone who work on DBT and compared TIGRE to some other software for validation). Incorrect geometry is going to cause problems in reconstruction, always.

[roshtha]

This was the TIGRE statement for computing geo.offOrigin and the only statement that uses Airgap.
geo.offOrigin =[((geo.sVoxel(1)/2)- (geo.DSD-geo.DSO)+Airgap);0;geo.sVoxel(3)/2]; % Offset of image from origin (mm) But with this statment, I could not get image outputs. (full dark images).

I replaced that with geo.offOrigin = [ 0;0;0];. So Airgap is not considered in reconstruction.

And yes, I see a different author for DBT reconstruction code.

Thanks.

[AnderBiguri]

@roshtha I still don't fully understand how your system works then. Is your center of rotation in the middle of the breast image?

[rodrigovimieiro]

@roshtha as DBT is half-cone beam, you must set the detector offset, otherwise, you will be using the CT geometry. The Origin offset is also useful to move the origin to the pivot, which is above the breast in your geometry.

Please, set these offset parameters and let us see how it is the geometry and why you are getting these black images.

You can use plotgeometry(geo,-pi/6); to see the geometry.

[roshtha]

Thanks @rodrigovimieiro for the reply. I will set detector offset and origin offest parameters and inform you of the progress.

[roshtha]

Hello,

I modified geo.offOrigin as geo.offOrigin = [((geo.sVoxel(1)/2)-(geo.DSD-geo.DSO)+Airgap); 0; 0];

and set geo.offDetector by refering comments from https://github.com/CERN/TIGRE/issues/99 geo.offDetector=[R*sin(angles); ones(1, size(angles, 2))*22];

I could see microcalcifications in reconstructed images. But the images have varying intensities among slices. Does this variation needs to be solved in image level? or is it a problem in setting geometry parameter?

The intensity variation in two continuous slices can be seen from figure below

A slice in which microcalcification is visible:

The geometry definition follows

geo.DSD = 650;  % Distance Source Detector (mm)
geo.DSO = 579;  % Distance Source Origin (mm)
geo.nDetector= [3584;2816];  % number of pixels (px)
geo.dDetector= [0.085; 0.085];  % size of each pixel (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)
geo.nVoxel= [171;1127;803];   % number of voxels(vx)
geo.dVoxel = [0.25; 0.25; 0.25]; % total size of the image(mm)
geo.sVoxel = geo.nVoxel .* geo.dVoxel;  % size of each voxel
Airgap = 17;  % DBT airgap (mm) 

geo.offOrigin = [((geo.sVoxel(1)/2)-(geo.DSD-geo.DSO)+Airgap);0;0];
geo.offDetector=[0; geo.sDetector(2)/2];  

geo.accuracy=0.5;
geo.COR=0;
nprojs = 16;     % Number of projections
tubeangle = 15;  % Angle range
angles=deg2rad(linspace(-tubeangle/2,tubeangle/2,nprojs)); % Angles in rad
geo.rotDetector=[0 0 0];

geo.mode='cone';    
geo=staticDetectorGeo(geo,angles);

The method plotgeometry(geo,0); shows figure A call to plotgeometry(geo,0) after staticDetectorGeo shows error.

Error using * Inner matrix dimensions must agree.

Error in plotgeometry>plotCircle3D (line 89) points=repmat(center',1,size(theta,2))+radius(v(:,1)cos(theta)+v(:,2)*sin(theta));

Error in plotgeometry (line 41)

So the above figure does not consider changes in parameter valuesgeo.offDetector and geo.rotDetector.

I refer this link https://github.com/CERN/TIGRE/issues/99 and modify input matrix data as

data = load('data1.mat');
img=data.img;
projections = -log(img./single(2^14-1)); %your maximum white level.
projections = rot90(projections);

I also use modifications from https://github.com/CERN/TIGRE/issues/99 for function staticDetectorGeo.

R=(geo.DSD-geo.DSO); %Radius of rotation of detector
geo.DSD = geo.DSD-(R - R*cos(angles));
geo.offDetector=[R*sin(angles); ones(1, size(angles, 2))*22];
geo.rotDetector=[zeros(1,size(angles,2));  zeros(1,size(angles,2));  -angles];

Raw data matrix shared at https://drive.google.com/drive/folders/10er-5LklkIFZc3BR0Z_AKhpyFiaivTIC?usp=sharing

I also attach reconstructed image outputs from machine.

machine_output.zip

Thanks.

[AnderBiguri]

Hum, to be honest, I can't answer why your intensity varies. Its not unheard of when the geometry is very far from the standard circular CT (like in DBT), but also its not unheard of when you have wrong geometry parameters. However, the fact that you can see small details on the image suggests that the geometry is somehow correct.

When you compare TIGRE result with the machine output, is it very different?

[rodrigovimieiro]

Thats good news. The MC seems in focus in that slice you showed us.

This intensity variating might be happening, if you are under Matlab, because when you display an image with imshow(slice, []) (I am supposing you are using it), the window level is adjusted automatically, based on the min and max value. From one slice to the other, these values might be changing in the image background, then it will cause this brightness problem.

You can set fixed values to the imshow, e.g. [100 300] (check your slice range), and display all the slices within this range.

Also, you can use ImageJ, and set the window level with a specific selected region.

Note, I am not saying that's exactly the problem. Your image looks fine for me in terms of things in focus.

Moreover, that are a lot of post processing in the images from the equipment.

[AnderBiguri]

@roshtha I am closing this issue as it seems to be quiet now. Do let me know if you have any more issues!