MATLAB Answers

1

How to remove artifact (“shadow” causing blurry areas) from 3D CT-data

Asked by Maria Huguett on 4 Apr 2019
Latest activity Commented on by Maria Huguett on 10 Apr 2019
Dear all,
I’m dealing with the segmentation of CT data, in order to get the porosity of my sediment sample. I already used the Otsu segmentation method, to get the threshold value and on this way determine how many voxels correspond to air and how many voxels are sediment grain. For sampling the sediment, I used paraffin, to get the sediment grains “stick” together, getting a compact and stable sample. After sampling, the sediment was scanned with x-ray. My sediment has also 3 phases: air (pores), grains and paraffin. However I already applied Otsu, considering paraffin as air. I would like to provide more information about my data:
My data is a 3D matrix of size 1981x1981x2009. The voxels of this matrix are gray intensity levels, which range from 0 till 65535 (16 bit). The whole matrix is not only sediment. The matrix is like a big 3D box and at the center of this box you can find the sediment. That means that sediment is surrounded of external air (black). I'm only interested in the sediment itself for porosity analysis. I attached an image (2D) of the sediment. It is a slice of the sediment in the Z-axis, at layer 990. In MATLAB I'm able to see slices of my sediment at different layers in each axis. When I get this slices (like the one attached), I can see some "shadows" in my sediment slice (or regions in my sediment, at which the sediment appears a little bit blurred, for example at the top). I guess, this could be artifacts, but I'm not sure. With Otsu as segmentation method, I’m getting big porosity values (which doesn’t make sense for my task). I guess this “shadow” is affecting Otsu, such that Otsu is not working well by choosing the threshold (This could be possible?)
Before looking for another segmentation method, I would like to remove this “shadow”/artifact of my whole 3D data and after that try the segmentation with Otsu again. Maybe after removing this artifact I get accurate Porosity results. Do you know a way to make this work for a 1981x1981x2009 big matrix? (and not only for one image, like the one attached here).
Or maybe I need a segmentation method, which could be more accurate, because considering 3 phases (as my sediment has grain, pores and paraffin)? If this is the case, you know a method to segment 3 phases and how to implement it in matlab? I would be very grateful for your help and advices, thanks a lot for your time!
slice_z_axis_layer_990.jpg

  0 Comments

Sign in to comment.

1 Answer

Answer by Sean de Wolski
on 9 Apr 2019

I'd use a multipass set of global and adaptive thresholds. You can adjust the sensitivities and thresholds as necessary.
I = imread('snip.png');
I = I(:,:,1); % gray
snip.png
Remove really small values
I(I<25) = 0;
Adaptive threshold
T = adaptthresh(I);
BW = imbinarize(I,T);
Now copy remove these and identify everything else.
I2 = I;
I2(BW) = 0;
BW2 = imbinarize(I2, "global");
Create label image.
L = uint8(BW2);
L(BW) = 2;
View it.
imshow(L, parula(3))
labelled.png

  1 Comment

Dear Sean,
thanks for your reply. I was trying your demo and it works. Now, I want to apply that for my data, so I added the first 4 lines. In the line: I(I<4500) = 0; At that line, I set the (random) gray value 4500 as threshold. When I set a bigger value, I get more pores (and less grains) and setting a small gray value there are less pores (and more grains). Is there a way that matlab calculates this grey value "automatically"? Because choosing it by myself affects the number of pores I will get drawn in the image.
fileID = fopen('sediment_sample.raw','r'); % my data
B = fread(fileID, [283*283*287 1], 'uint16'); % Read 283*283*287 matrix elements in a column vector
C = reshape(B, 283, 283, 287); % Reshape B to a 3D Matrix
I=C(:,:,143); % Get a slice of the sediment at Layer 143
I(I<4500) = 0;
% Adaptive threshold
T = adaptthresh(I);
BW = imbinarize(I,T);
I2 = I;
I2(BW) = 0;
BW2 = imbinarize(I2,'global');
% Create label image.
L = uint16(BW2);% my data are 16 bit (uint16)
L(BW) = 2;
%View it.
imshow(L, parula(3))
this is my output. Here I get two colors: purple for pores and yellow for grains. But how can I get the color corresponding for paraffin (like the blue one at your example)? Thanks a lot :)
sediment_segmentation.jpg

Sign in to comment.