Applying Gaussian 3D filter to 32 bit float 3D Micro-CT image

Hi All,

I have a Micro-CT data file to process it through Matlab as shown in below. My array is a 3D array and it is 32 bit float. but gaussian filters which are introduced in mathworks are not supported to float data type

fid = fopen('E:\RE\re_C_7mm_6um_1d42s_01\re_C_7mm_6um_1d42s.vol');
voxels = fread(fid, '*float32');  %read as 32 bit float. I'm assuming they are 32 bit float(CT data)
fclose(fid);
voxels = reshape(voxels,[1920,1920,1536]);
d=imbinarize(voxels,95); % Thresholding the volume
d1=imgaussfilt3(d)  % This  lines gives me an error that it supports only for uint,int, single or double.check the error

What I should do to filter my 3D array using the gausfilt3 command? Appreciate your valuable advices. I want to remove the P(see the attached image)with or without gausfilt3.Appreciate your valuable adivces!

Error is shown below. A is the file for the filtering. If you check the gaussfilt3 command in mathworks A is defined. it is not my own name.

Error using imgaussfilt3
Expected A to be one of these types:

uint8, uint16, uint32, int8, int16, int32, single, double

Instead its type was logical.

Error in imgaussfilt3>parseInputs (line 285)
validateattributes(A, supportedClasses, supportedImageAttributes, mfilename, 'A');

Error in imgaussfilt3 (line 116)
[A, options] = parseInputs(args{:});

Error in volreading (line 55)
d1=imgaussfilt3(d)
 

Answer 1

[Image Analyst]

Try imfilter() with your own custom kernel that is a Gaussian shape.

[Hege]

Thank you so much Image Analyst. I tried. But there is no good results. code is successfully processed for the [3 3 3] and [9 9 9] kernel sizes as well. But I am not sure the way what I am doing to remove the P(in the Image) is correct. As you know, in the Micro-CT process, fibre sample is rotated 360 degrees.because of that P is created(probably it is air mixed with noise). Actually, I was struggling to remove P by applying different techniques such as medfilt(3), bwareaopen etc. Nothing worked.Finally, imfilter() option just tried and did not give good results. I really appreciate if you could please advice me!!

My expectation is to get a clear 3D image after removing the noise and then skeletonize the fibres(skeletonization to get the medial axises of fibres).

Really appreciate your valuable advices!

[Image Analyst]

Not sure what you consider noise. Have you tried bwareaopen() to remove small particles? Do you just want to smooth the fiber surface by blurring the binary image and re-thresholding?

[Hege]

Thank you so much for your kind attention Image Analyst. I really appreciate it.

I tried bwareaopen().But it didnt give good results for the 3D Volume(3D array) which means it was difficult to identify whether atleast there was a change. Actually, I want to remove everything, except fibres.That is my expectation. I tried lots of operations.medfilt3(), bwareaopen(),imfilter() etc. Actually, below command worked for the thresholding. It removes the voxels values below 95. I got this '95 value' after doing some experiments.

d=imbinarize(voxels,95); % Thresholding the volume

As a result of this, voxels values higher than 95 is remained and there are many things that I want(fibres) and that I dont want(rest of the things).

I am just thinking one thing now, Can I threshold a volume using a upper value and a lower value using the 'imbinarize' command? As an example, I want the voxels between 95-105.So I can check what is happening atleast.

Your valuable suggestions are highly appreciated!

[Image Analyst]

You can threshold on intensity of course. But as of now, there is no bwareafilt() equivalent for 3-D to threshold based on blob size. But you can get the volumes with regionprops3() and find volumes in a certain range and get rid of them that way.

[Hege]

Thank you very much Image Analyst!

I tried regionprops3() and I received a huge set of data(Please see the attached image. I attached a small portion of the table as the table size is 863618x3 table. Reason is fibers are broken.That is why I received this amount of voxels.Actually, I want to remove the cylindrical shape(P) and the sample holder. So, I am thinking about these voxels numbers to identify and how to remove it. But I noted that only one voxel value has 311793705(first value in the image table).Others are comapratively small. If you check the image you will notice it.Thinking how I should use the table to remove the unwanted things using these values.

Actually, doubt in my head is, I generated the table using below commands.

d=imbinarize(Micro-CT file,114.6);

e=regionprops3(d);

but in the table there voxel values less than 114.6 as well.something 1,2 etc.

In the mean time, I found an experimental solution but worrying how I should apply it to the 3D Image. If I can binarize the voxels between 114.6-220 clearly, I can remain the fibers and others are gone.Something like a badpass filter.What I am trying to say is I want the voxels between 114.6 and 220. As a result of this, voxels out of the margin should be zero. Basically, thinking about how to apply a lower threshold value and upper threshold value to my 3D array.

Any of your suggestions are highly appreciated!.Thank you in Advance.

[Image Analyst]

Is your threshold on blob volume? Or intensity value?

[Hege]

Thank you so much for your attention Image Analyst. Actually, I found those values using the myVGL software which is supported to the Micro-CT system. So, it is based on the intensity value.

[Hege]

Hi Image Analyst, Appreciate if you could please give me an advice.

If we want to convert a 2D image to double data type from different data type(ex=uint8) we can use the 'im2double' .but what we can use for the 3D image if we want to convert the 3D image to double type from 32 bit float?

Currrently, I am developing a script to with upper and lower threshold values. Appreciate your suggestions. But I have a problem with covnverting the 3D image in to double data type. I really appreciate your valuable suggestions and advices regarding what I am developing now.Thank you in Advance!

d=voxels;
%d1=im2double(d);------------------------I have a doubt in here for the 3D array??
for i=1:1920;
    for j=1:1920;
        for k=1:1536;
            if (i&&j&&k)> 95
              d(i,j,k)=1;
             else  d(i,j,k)=0; 
              if (i&&j&&k)<220
                   d(i,j,k)=1;
                else  d(i,j,k)=0;
              end
             end
       end
    end
    end

[Image Analyst]

(i&&j&&k) will just be true (1) or false (0), and in this case since i, j, and k all more than 1, it will always be true, so not sure what you're doing there. To convert voxels to double, you can do

d = double(voxels);

[Hege]

Thank you Image Analyst.

Honestly, I tried the d=double(voxels) command first. may be then something wrong in my script.

Actually, what I am trying to do is ,voxels values more than 95 should be 1 and if not it should be zero. Further, voxels values less than 220 should be 1 and if not it should be zero. That is what I am trying to do at the moment. Appreciate if you can show if there is any mistake of way I am trying to do thresholding.Thank you

[Hege]

Thank you for your support Image Analyst. I thresholeded the 3D volume and heading for the skeletonization. But I got the values from Micro-CT experimental works. Thank you in Advance