Find image intensity and perform fitting - fminsearch

Hello,

I am looking to:

1. Extract intensity along the length of the medium ( xray scan) ( Not sure how)

2. Plot the intensity

3. And use fminsearch to fit the intenisty values in a wave propagation model to determine the constants in the model ( Not sure how)

providing the tiff image link,. since I am unable to upload here: https://ibb.co/cczb54D

My attempt is provided below. Kindly requesting help.

% Objective: 
 %{
   1. Import a CT scan image( CT scan of an Cu wire)
   2. Extract the intensity from the CT scan image(CT scan of an Cu wire) :
   which will be the observed data
   3. Define the exponential decay model describing the Wave propagation through a medium
   4. Perform fminsearch to determine the constants ( imaginary constants)
   in the decay model.
%}

clc;    % Clear the command window.
close all;  % Close all figures 
format long g;
format compact;


%% Importing the CT scan/test image

% fullname = get_full_filename('20220422_5anglestep000.tif')
test_image = imread('test_image.tif');
imshow(test_image)

% Getting the dimensions of the image.
[rows, columns] = size(test_image);


%%  Display histogram

imhist(test_image);
grid on;
title('Histogram of original test image')

%% Plotting intensities along Z axis
z1 = [0 size(test_image,2)];
y1 = [size(test_image,1)/2 size(test_image,1)/2];
c = improfile(test_image,z,y1);



%% Extracting the intensities from the medium

% Image represents the wave propgation through a medium in terms of
% intensity vs thickness(displacement) of the medium


% Looking to extract the intensity as a function of displacment from image

% z: displacement ( Z axis locations)
% wave_obs = intensity (Y axis)




 
%% Defining the equation wave = exp(i(1-del)k*z * exp(-beta)*z

%{
i -  signifying imaginary part
z -  array representing the length of the medium
k -  refractive index
del - constant ( to be determined)
beta -constant ( to be determined)
%}

k=1.5;

%% fminsearch
wave_pred =  fminsearch(wave_prediction);

predicted_wave = wave_pred(z,k);






function wave_prediction = wave_pred(z,k)
    del =5;
    beta =7;
    k=1.5;
    z= 0:0.2:20;
    wave1 = exp(-1i*(1-del)*k*z);
    wave2 = exp(-beta)*z;
    wave_prediction = @(del,beta)wave1.*wave2;
end

Answer 1

[Image Analyst]

You don't need improfile() if you want the intensity along a certain row, so instead of

c = improfile(test_image,z,y1);

do

middleRow = round(size(test_image, 1) / 2) % Find middle row.
c = test_image(middleRow, :) % Get all columns in the middle row.
% Now plot the intensity.
plot(c, 'b-', 'LineWidth', 2);
grid on;
xlabel('X Column');
ylabel('Gray Level');

[Anand Ra]

@Image Analyst Thank you!

The image I have is a CT scan. The middle row is dark, which I think indicates attenuation.Therefore I would expect a decay in the intensity.

Btw,

1. what does the x column signify?
2. Do I have to switch the image ( Swap X and Y) to expect a decay?

[Image Analyst]

@Anand Rathnam not sure I'm following you. x is the horizontal direction and corresponds to columns, or the second dimension when referencing arrays. y is the vertical direction and corresponds to rows, or the first dimension when referencing arrays. So it's array(row, column) or array(y, x) NOT array(x,y). Now some functions like improfile() and plot() use x,y instead of row, column. But other functions, like bwboundaries() and find() return coordinates in (row, column) which is (y, x) instead of as (x,y). Bottom line is you have to be careful when using coordinates to get them right. A very common beginner mistake is to try to get a value like

value = array(x, y)  % WRONG!

instead of

value = array(y, x)  % RIGHT

Your "z" corresponds to x in your improfile call. And your y is y, or row. Your y (row) value is halfway down the image. You z is 0 for the first point, which is not allowed. It has to be 1 since column 1 is the first column. The last z, z(2) is the far right edge of the image. So your call to improfile goes across all columns of the image in the middle row of the image.

[Anand Ra]

@Image Analyst I think that explanation is definitely helpful. I guess I got mixed up.

Xray beam passing through the specimen decays(intensity) due to the change in refractive index ( from air to specimen).

To determine the intensity decay, I think histogram would make sense.

Therefore, I tried plotting the histogram. And I see the expected decay in the intensity.

%% Display histogram

figure(2)

imhist(test_image);

grid on;

title('Histogram of original test image')

Therefore, I tried obtaining the values of those historgam intensities( sum of the intensities)

[ counts , binLocations ] = imhist( test_image )

figure(8)

plot(binLocations, intensity)

My understanding is that these above calculated intensities are along the thickness of the specimen

Is it possible to plot: intensity vs distance( thickness of the medium)

1. intensity is obtained with above command( Hope I am right)
2. How do I calculate the distance( thickness of the medium) - from the number of columns / binlocations?

Ideally looking for the two vectors from the image, so I can perform the data fitting

vector 1 : intensity values ( decaying)

vector 2: Length/Thickness of the medium.

[Image Analyst]

You need to put a calibration standard in there, like an aluminum step wedge, so that you can calibrate gray level vs. thickness.

[Anand Ra]

@Image Analyst Thank you!