Matlab - Raman Mapping

guys, please can you help me in processing matlab code with my data (txt mapping spectra+ filter image for recognition+ that i upload in

https://www.uschovna.cz/zasilka/SYTLSNMDRNVX5WAJ-6BX/

https://www.uschovna.cz/zasilka/SYEGP5642Z8YHCFF-R5X/ (in this link i added also that matlab scripts that i wrote also in this message)

and you can easily free download them on that storage web and use it for main script in "folder".

My mapping code with pca elements and clustering is almost ready, he is fully functional and process, but i need to change some parameters but i dont know how to change it to make it working right to see all clusters for all the particles i need (assume there are 3 or 4 clusters for spectra of Polyethylene (PE), Polyethylene terephtalate (PET), Polytetrafluoroethylene (PTFE) and maybe Polypropylene (PP) ) and proper heat map for all the particles in the image that could be signalising as a spectra look like microplastic. For mapping code i also add 3 additional functions, they are fully functional and need for operating my main mapping code. please could help me urgently, need to answer it asap

My main Raman mapping code:

%% Main mapping script
clear all; close all;
folder               = ''; % write your own file path of my txt spectra
map_files            = dir([folder,filesep,'22.9.23_Rajec_Mapping_100pr_10x_10s_1x*.txt']); %you can leave it this way without changing
main_peak_position   = 3000;
limit_x_axis         = 0;
limit_x_axis_values  = [150 3300];
for i = 1:length(map_files)
    filename = map_files(i).name;
    [X_position, Y_position] = Extract_XY_positions(filename);
    map_files(i).X_position = X_position;
    map_files(i).Y_position = Y_position;
    A = readmatrix([folder,filesep,filename]);
    if limit_x_axis == 1
        index_lower = find(A(:,1) < limit_x_axis_values(2), 1);
        index_upper = find(A(:,1) < limit_x_axis_values(2), 1, 'last');
        map_files(i).Raman_shift     = A(index_lower:index_upper,1);
        map_files(i).Raman_intensity = A(index_lower:index_upper,2);
    else
        map_files(i).Raman_shift     = A(:,1);
        map_files(i).Raman_intensity = A(:,2);
    end
end
%% zjisteni rozmeru mapy
x_positions = sort(unique(cat(1, map_files(:).X_position)));
y_positions = sort(unique(cat(1, map_files(:).Y_position)));
x_positions = [x_positions (1:length(x_positions))'];
y_positions = [y_positions (1:length(y_positions))'];
main_peak_matrix       = zeros(length(x_positions), length(y_positions));
peak_ratio_matrix      = zeros(length(x_positions), length(y_positions));
naive_main_peak_matrix = zeros(length(x_positions), length(y_positions));
%% check of individual spectrum
show_figures = 1;
plot(A(:,2))
for i = 7:length(map_files)
    if show_figures == 1
        check_figure('original','Original Spectrum')
        plot(map_files(i).Raman_shift, map_files(i).Raman_intensity)
        title(['X = ', num2str(map_files(i).X_position),' \mum; Y = ', num2str(map_files(i).Y_position), '\mum']);
    end
    y_filter = mslowess(map_files(i).Raman_shift, map_files(i).Raman_intensity);
    vlnocet_spodni_background = 967;
    vlnocet_horni_background  = 1900;
    x = map_files(i).Raman_shift;
    peak_loc_index = find(x < main_peak_position, 1);
    background_mask = x < vlnocet_spodni_background | x > vlnocet_horni_background;
    pp = polyfit(x(background_mask), map_files(i).Raman_intensity(background_mask), 1);
end
%% plot all spectra 
num_spectra     = length(map_files);
all_Raman_spectra = zeros(length(map_files(1).Raman_intensity), num_spectra);
for i = 1:num_spectra
    all_Raman_spectra(:,i) = map_files(i).Raman_intensity;
end
check_figure('allSpectra','All Spectra')
plot(map_files(1).Raman_shift, all_Raman_spectra)
%% vyhodnoceni mapy
check_figure('native_Raman')
main_peak_matrix = zeros(length(x_positions), length(y_positions));
AAA = (1:length(map_files))';
for i = 1:length(map_files)
    map_index = [ find(map_files(i).X_position == x_positions), ...
        find(map_files(i).Y_position == y_positions) ];
    AAA(i,2) = map_index(1);
    AAA(i,3) = map_index(2);
    main_peak_matrix(map_index(1), map_index(2)) = map_files(i).Raman_intensity(200);
end
show_Raman_image(main_peak_matrix)
colorbar
%% obraz
B = main_peak_matrix;
check_figure('obraz')
imagesc(B)
%% spektra
check_figure('spektra')
for i = 1:num_spectra
    check_figure('spektra')
    plot(map_files(1).Raman_shift, map_files(i).Raman_intensity)
    set(gca,'YLim',[0 25000])
    map_index = [ find(map_files(i).X_position == x_positions), ...
        find(map_files(i).Y_position == y_positions) ];
    check_figure('obraz')
    show_Raman_image(main_peak_matrix)
    hold on
    axis on
    plot(map_index(1), map_index(2), 'r+', 'MarkerSize', 30, 'LineWidth', 2);
    axis off
    hold off
end
%% mean spectrum + peaks
check_figure('summary','All Spectra')
mean_spectrum = mean(all_Raman_spectra,2);
plot(map_files(1).Raman_shift, mean_spectrum)
[pks,locs] = findpeaks(mean_spectrum,'MinPeakProminence',5);
plot(mean_spectrum)
hold on
plot(locs,pks,'r*')
hold off
length(locs)
%% PCA testing
all_Raman_spectra_peaks = transpose(all_Raman_spectra(locs,:));
[coeff,score,latent,tsquared,explained,mu] = pca(all_Raman_spectra_peaks);
plot(coeff(:,1),coeff(:,2),'h')
XXX = [coeff(:,1),coeff(:,2)];
[idx,C] = kmeans(XXX,4);
hold on
plot(C(:,1),C(:,2),'kx','MarkerSize',15,'LineWidth',3)
hold off
%% GMM
plot(score(:,1),score(:,2),'*')
XXX = [score(:,1),score(:,2)];
gm = fitgmdist(XXX,3);
idx = cluster(gm,XXX);
gscatter(XXX(:,1),XXX(:,2),idx);
%% cluster map
check_figure('cluster_Raman')
cluster_map = zeros(length(x_positions), length(y_positions));
for i = 1:length(map_files)
    map_index = [ find(map_files(i).X_position == x_positions), ...
        find(map_files(i).Y_position == y_positions) ];
    cluster_map(map_index(1), map_index(2)) = idx(i);
end
show_Raman_image(cluster_map)
colorbar
%% cluster mean spectrum
cluster_type = 4;
check_figure('cluster_mean_spectrum',['Mean spectrum in cluster ',num2str(cluster_type)]);
selected_spectra_mean = mean(all_Raman_spectra(:,idx == cluster_type),2);
plot(map_files(1).Raman_shift, selected_spectra_mean)

additional essential scripts- functions: show_Raman_image, Extract_XY_positions, check_figure

Extract_XY_positions

function [X_position, Y_position] = Extract_XY_positions(filename)
% return values of position in microns
underscore_positions = strfind(filename,'_');
pos_X               = strfind(filename,'Xµm_');
pos_Y               = strfind(filename,'Yµm_');
diff_vec_x = underscore_positions - pos_X;
diff_vec_x(diff_vec_x<0) = NaN;
[~, min_pos_X] = min(diff_vec_x);
X_position = str2double( filename(underscore_positions(min_pos_X)+1 : underscore_positions(min_pos_X+1)-1) );
diff_vec_y = underscore_positions - pos_Y;
diff_vec_y(diff_vec_y<0) = NaN;
[~, min_pos_Y] = min(diff_vec_y);
Y_position = str2double( filename(underscore_positions(min_pos_Y)+1 : underscore_positions(min_pos_Y+1)-1) );
end

check_figure

function check_figure(figure_name,varargin)
% figure open a new figure, with known description.
% it will always open the figure (the already open one or they will produce a new one)
% check_figure(figure_name) is only to open a figure with this handle
% check_figure(figure_name,figure_description) opens the figure and make a
% string description saved in variable figure_description
command_to_do = sprintf('exist(''%s'',''var'')',figure_name);
if evalin('caller',command_to_do) == 0
    eval(['figure_name2',' = figure;']);
else
    figure_name2 = evalin('caller',figure_name);
end
switch nargin
    case 1
        if isvalid(figure_name2) == 0
            figure_name2 = figure;
        else
            figure(figure_name2);
        end
    case 2
        if isvalid(figure_name2) == 0
            figure_name2 = figure;
            figure_name2.Name = varargin{1};
        else
            figure(figure_name2);
            figure_name2.Name = varargin{1};
        end
end
assignin('caller',figure_name,figure_name2)

show_Raman_image

    function show_Raman_image(image_to_show)
        imagesc(flipud(rot90(image_to_show)))
        colormap('jet')
        v = image_to_show(:);                  % take our pixels
        lims = prctile(v, [1 99]);             % or [min(v), max(v)]
        caxis(lims)
    end

P.S sorry , dont exactly remember what "products" are using, as i remember image processing toolbox, chemometric toolbox and smth, if you run the code and you dont have one, i think its no problem to install from matlab for 5 minutes

Also my cod were functional also before about 3-4 years ago, so elder versions of matlab will easily process scipts

[Mathieu NOE]

hello @Nazir

we have been in contact recently

see function ram1 which is your main code with some comments

all the best

[sudobash]

Hi @Nazir, can you summarize what exactly the issue is that you are facing? Thanks!