Non-rectangular thermal image crop & matching

How to create rectangular IR image and contour with a false color map from the measured non-rectangular image and data from an Excel sheet?
I want to transform the results of the IR image and temperature data from the Excel sheet into a rectangular format. The interrogation region from the measuring devices is non-rectangular. Procedures are attached in the PDF file. I need help matching rectangularized IR images having the original color map with Excel data.

 Accepted Answer

There are two CSV files. One is readable. The other is a duplicate with broken formatting. It should have simply been deleted. There's no point trying to fix it, because the underlying data is exactly identical.
Colormap estimation from a crappy JPG. Data estimation from pseudocolor images. Perspective correction. We've done this before.
% read the file
inpict = imread('thermal.jpg');
% extract a noisy estimate of the colormap
CT0 = imcrop(inpict,[568.51 49.51 30.98 379.98]);
CT0 = reshape(mean(im2double(CT0),2),[],3);
CT0 = flipud(CT0);
% the colorscale limits specified in the image
Trange = [22.82 45.42];
% estimate T from a JPG
datafromjpg = rgb2ind(inpict,CT0,'nodither');
datafromjpg = rescale(datafromjpg,Trange(1),Trange(2));
% these are the coordinates of the box corners
% you can get these using getpts() or impixelinfo() or datatips
boxm = [82 19; 510 25; 488 464; 81 426]; % [x y]
% assert that this is where they're supposed to be
% any coordinates that define a rectangle
boxf = [1 1; 512 1; 512 512; 1 512]; % [x y]
% transform/crop the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d([512 512 3]);
datafromjpg = imwarp(datafromjpg,TF,'fillvalues',255,'outputview',outview);
subplot(1,2,1)
imagesc(datafromjpg)
colormap(parula(256))
cb = colorbar('southoutside');
cb.Ticks = clim();
subplot(1,2,2)
histogram(datafromjpg)
The CSV file appears to be the sensor data from the same camera, but without the annotations and destructive JPG compression. It has the same geometry, so the same transformations apply.
% read the file
datafromcsv = readmatrix('good-vip10t,single=IMG_20231219_214630_TEMP.csv');
% these are the coordinates of the box corners
% you can get these using getpts() or impixelinfo() or datatips
boxm = [82 19; 510 25; 488 464; 81 426]; % [x y]
% assert that this is where they're supposed to be
% any coordinates that define a rectangle
boxf = [1 1; 512 1; 512 512; 1 512]; % [x y]
% transform/crop the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d([512 512 3]);
datafromcsv = imwarp(datafromcsv,TF,'fillvalues',255,'outputview',outview);
figure
subplot(1,2,1)
imagesc(datafromcsv)
colormap(parula(256))
cb = colorbar('southoutside');
cb.Ticks = clim();
subplot(1,2,2)
histogram(datafromcsv)
As to whether those two images are indeed taken at the same time, or whether there's some other explanation for the largest chunk of the error, nobody will ever know.

6 Comments

Hi, DGM.
Thank you very much for your kind comments and code. It’s very helpful. Your comments helped me find the default settings in the thermal camera.
I’d like to cross-check the temperature results obtained by your method with the Excel data from the IR camera. Attached are visible and thermal images as well as trial code. Could you give me some comments and suggestions?
You're picking up some fill values in the imwarp() transformation that are skewing the apparent scale of the data. I just chose to inpaint those. It's just a few values on the image boundary, so they could also just be cropped off.
% read the file
inpict = imread('check-gap_vip10t,vip10t-alignkey=IMG_20231219_215650_IR.jpg');
% extract a noisy estimate of the colormap
CT0 = imcrop(inpict,[568.51 49.51 30.98 379.98]);
CT0 = reshape(mean(im2double(CT0),2),[],3);
CT0 = flipud(CT0); % CT0 = fliplr(CT0);
% the colorscale limits specified in the image
Trange = [22.82 45.42];
% estimate T from a JPG
datafromjpg = rgb2ind(inpict,CT0,'nodither');
datafromjpg = rescale(datafromjpg,Trange(1),Trange(2));
% these are the coordinates of the box corners
% you can get these using getpts() or impixelinfo() or datatips
boxm = [78 1; 77 480; 492 480; 516 2];
% assert that this is where they're supposed to be any coordinates that define a rectangle.
boxf = [1 1; 640 1; 640 480; 1 480]; % [x y]
% transform/crop the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d([480 640 3]);
datafromjpg = imwarp(datafromjpg,TF,'fillvalues',NaN,'outputview',outview);
% some pixels are missing (filled with 'fillvalues')
% we could either inpaint those or crop them off
datafromjpg = regionfill(datafromjpg,isnan(datafromjpg));
subplot(1,2,1)
imagesc(datafromjpg)
colormap(parula(256))
cb = colorbar('eastoutside'); % cb = colorbar('southoutside');
cb.Ticks = clim();
subplot(1,2,2)
histogram(datafromjpg)
%% ---------------------------------------
% Read the table from a CSV file
tbl = readtable("check-gap_IMG_20231219_215650_TEMP_01title.csv");
txt = tbl{:,:};
txt = erase(txt,'"');
txt = cellfun(@str2num,txt,'UniformOutput',false);
datafromcsv = cell2mat(txt);
% transform/crop the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d([480 640 3]);
datafromcsv = imwarp(datafromcsv,TF,'fillvalues',NaN,'outputview',outview);
datafromcsv = regionfill(datafromcsv,isnan(datafromcsv)); % same inpainting
subplot(1,2,1)
imagesc(datafromcsv)
colormap(parula(256))
cb = colorbar('eastoutside');
cb.Ticks = clim();
subplot(1,2,2)
histogram(datafromcsv)
If you want suggestions, don't create a new variable name on each line. Keeping a bunch of similarly-named unused intermediate results around only wastes memory, makes the code difficult to read, and will eventually cause a mistake, since you have to keep track of six copies of the data instead of two. I don't remember, but I think there was one of the wrong copies being used in a call to imwarp(), but I don't remember which one.
Hi, DGM.
Thank you so much for your kind suggestions and code. The code has been modified to compare temperature profiles along a central vertical line. Comparing the temperatures (Figure 7 in the code) shows the difference between the Excel data from the IR camera and the data obtained with your method.
Could you please give me some suggestions and correction tips for improvement in the attached code?
You're simply going to have a hard time getting a good comparison between the two images. The JPG copy has been subject to multiple layers of degradation. It's been quantized in the process of colormapping, then again in JPG encoding, 75% of the overall chroma information (hue and saturation) has been discarded. Then we read it back in and we rely on a small portion of the image (the colorbar) to estimate the original colormap. From that estimate, we again estimate the discrete graylevels of the original image. Due to the spatial downsampling in chroma, pixels in the original image which were once the same graylevel may no longer be represented by the same color, and neither may map closest to the corresponding pixels in the colorbar. The artifacts created near the boundary of the colorbar make the map estimation error tend to be high near the extrema, which in turn tends to skew the whole mapping. I'm kind of surprised that the profiles look as close as they do.
I've heard that some cameras use a nonlinear mapping, but I don't suspect that here. Your CSV appears to be directly read values, though you might double-check that. It wouldn't make sense that the pseudocolor image would be nonlinear, otherwise it would not be practical to read with only two ticklabels.
So is it plausible that the error can be entirely explained by the degradation in the JPG copy? I'm not sure. I really don't know anything about the camera's behavior. If the images are captured sequentially, and the system isn't at steady state, that might explain part of it.
We can try to see if there's a roughly linear relationship between the two pixel populations.
% assume datafromcsv is the reference
% start with a downsampled copy to save time
ks = 0.25;
A = imresize(datafromcsv,ks);
B = imresize(datafromjpg,ks);
[Au,~,idx] = unique(A);
Bu = zeros(size(Au));
for k = 1:numel(Au)
idx = A == Au(k);
thesepix = B(idx);
Bu(k) = median(thesepix);
end
plot(Au,Bu)
grid on
axis equal
xlabel('Values from CSV')
ylabel('Values from JPG')
So it's sort of linear. Is that truly nonlinear, or is that just noise? I don't have a good answer for that.
For the profile line, you might consider sampling a narrow window around that line and then averaging across its width. That might help suppress some of the noise. Alternatively, you might consider just filtering the entire images with some sort of edge preserving filter prior to sampling the profile.
% Temperature profile along the central vertical axis from the bottom.
hw = 10; % window half-width
yaxis_x0 = 320 + (-hw:hw);
Tcsvvertical_y = mean(datafromcsv(:,yaxis_x0),2);
Tjpgvertical_y = mean(datafromjpg(:,yaxis_x0),2);
y240pxl = 480-y;
plot(y240pxl,Tcsvvertical_y,'ro-',y240pxl,Tjpgvertical_y,'bo-')
title('Verticle temperature profile [deg.C]');
xlabel('Distance verticle [pxl]');
ylabel('Temp. vertical')
legend('T_{csv}vertical(y)','T_{jpg}vertical(y)')
... though I don't think it will subdue the larger portion of the difference.
I think the more important question is whether it's necessary to use what we know is a crude estimate (the JPG data) if we have the CSV data. I'm normally inclined to treat pseudocolor images like any graph. It's a visualization of data, but it shouldn't be taken as a reliable store of data.
Thank you so much for the detailed and touching explanation.
Could you please tell me about your favorite reference book on thermal imaging using MATLAB?
I actually don't have any references for thermal imaging. @Image Analyst might know of some good sources, but I'm not the resident professional image processing guy here. Everything I say is based on a casual history of image manipulation in MATLAB and a background in general engineering principles.
That said, you might want to see how much technical documentation your camera manufacturer has available to owners. You appear to have a relatively high-end camera, so with that price might come some better documentation than your typical consumer-market offering.

Sign in to comment.

More Answers (1)

Hello bk park
From what I gather, you are trying to read data from a CSV file and then converting the non-rectangular interrogation region of the temperature measured to a rectangular shape.
This can be achieved by reading the data and first preprocessing it into a usable form. Further, you’d need to isolate the non-rectangular interrogation region using some threshold. Then, the non-rectangular region can be changed to a rectangular shape using padding with 0’s. Another method can be data interpolation. The following code snippet shows how to extract usable data, then the interrogation region and then a rectangular region using both padding and interpolation.
% Read the table from a CSV file
tbl = readtable("good-vip10t,single=IMG_20231219_214630_TEMP_01title.csv");
% Initialize an empty array to store processed temperature data
b = [];
% Loop through each row of the 'Tempeature_IR' column in the table
for i = 1:size(tbl.Tempeature_IR,1)
% Extract the temperature string for the current row, removing quotes
a1 = tbl.Tempeature_IR{i};
a1 = a1(a1~='"');
% Split the string into individual temperature values
a1 = strsplit(a1);
% Convert the temperature strings to numeric values
a1 = str2double(a1);
% Append the numeric temperature values to the array 'b'
b = [b; a1];
end
% Remove the last column from 'b'
b = b(:,1:end-1);
% EXTRACTING INTERROGATION REGION USING THRESHOLDING
% Filter columns of 'b' to keep only those where the first row value is >=
% 40. You can change the threshold.
colsToKeep = b(1, :) >= 40;
result1 = b(:, colsToKeep);
% Initialize the size of the 'result1' matrix
[rows, cols] = size(result1);
% Iterate through each row of 'result1'
for i = 1:rows
% Find the last index in the row where the value is >= 38
idx = find(result1(i, :) >= 38, 1, 'last');
% If such an index is found, set all subsequent elements in the row to 0
if ~isempty(idx)
result1(i, idx+1:end) = 0;
end
end
% PADDING THE OUTSIDE REGION WITH 0's
% Create a mask for values in 'result1' that are >= 38
mask = result1>=38;
% Initialize 'maskedresult' as a copy of 'result1'
maskedresult = result1;
% Set all values in 'maskedresult' that meet the mask condition to 0
maskedresult(mask)=0;
% INTERPOLATION
% Create a meshgrid for the coordinates of 'maskedresult'
[x, y] = meshgrid(1:size(maskedresult,2), 1:size(maskedresult,1));
% Identify indices of zero and non-zero values in 'maskedresult'
zeroInd = find(maskedresult == 0);
nonZeroInd = find(maskedresult ~= 0);
% Extract coordinates and values for interpolation
xq = x(zeroInd); % x-coordinates of zeros
yq = y(zeroInd); % y-coordinates of zeros
xv = x(nonZeroInd); % x-coordinates of non-zeros
yv = y(nonZeroInd); % y-coordinates of non-zeros
v = maskedresult(nonZeroInd); % non-zero values
% Use 'scatteredInterpolant' for interpolation with linear interpolation
% and nearest extrapolation for out-of-bound points
F = scatteredInterpolant(xv, yv, v, 'linear', 'nearest');
% Interpolate values for positions that were zero in 'maskedresult'
interpolatedValues = F(xq, yq);
% Replace zero values in 'maskedresult' with the interpolated values
maskedresult(zeroInd) = interpolatedValues;
% Display the matrix after interpolation
disp('Matrix after interpolation:');
disp(maskedresult);
The above operation give the following heatmaps:
Picture 1: Padded with 0’s
Picture 2: Interpolated.
For further understanding, refer to the links to the MATLAB documentation given below:
  1. scatteredInterpolantfunction- https://www.mathworks.com/help/releases/R2022b/matlab/ref/scatteredinterpolant.html
I hope you find the above explanation and suggestions useful!

1 Comment

Hi, G. Pant.
Thank you very much for your kindness and code. It’s very helpful.

Sign in to comment.

Products

Release

R2022b

Asked:

on 15 Feb 2024

Commented:

DGM
on 11 Mar 2024

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!