Map data from a conical surface in an image to a rectangle in x-theta coordinates.

Hi all,

I have 3 images taken from different azumithal angles: top (0 deg), side (90 deg), and bottom (180 deg) (see attached). I'm writing a script to average the images over time, crop the cones out of each image, "unwrap them" to get them in a rectangular shape, and then stitch them together to get a compositie image on all sides. So far, I'm able to:

1. Import the series of images, average them over time to get a mean image from the top, bottom, and side cameras.
2. Crop the cone from the image
3. Use edge detection to get x-y coordinates from the cropped cone

What I need help with most is the "unwrapping" part to covert the image's cartesian coordinates in x-y to x-theta.

Below my code so far.

Thanks!

% Oil Flow Analysis Code

clc; clear; close all;

oil = 'D:\PSWT Data_June 2023\Oil Flow\';

addpath(oil);

%% Add file paths

top = [oil 'TopCam']; % top images

bottom = [oil 'BotCam']; % bottom images

side = [oil 'SouthCam']; % side on images

addpath(top,bottom,side);

%% Set default parameters

fps = 1000; % fps

Pref = 14.7; % reference wind-off pressure (psi)

runn = 'Run6499';

%% Load image sets

% Load images

top_i = imageDatastore([top '\' runn '\*.tif']);

bot_i = imageDatastore([bottom '\' runn '\*.tif']);

side_i = imageDatastore([side '\' runn '\*.tif']);

% Load color lists

top_color = [top '\' runn '\' runn '_color.lst'];

bot_color = [bottom '\' runn '\' runn '_color.lst'];

side_color = [side '\' runn '\' runn '_color.lst'];

% Get number of frames for top, bottom, and side images

n_top = numel(top_i.Files);

n_bot = numel(bot_i.Files);

n_side = numel(side_i.Files);

%%% Import image matrices %%%

% Load top image matrices

for i = 1:n_top

[im_top{i}] = imread(string(top_i.Files(i)));

end

% Load bottom image matrices

for i = 1:n_bot

im_bot{i} = imread(string(bot_i.Files(i)));

end

% Load side image matrices

for i = 1:n_side

im_side{i} = imread(string(side_i.Files(i)));

end

%% Get mean image and covert to 16-bit

% Average image matrices over last 100 frames when oil flow is steady

top_m = mean(cat(3,im_top{end-100:end}),3);

bot_m = mean(cat(3,im_bot{end-100:end}),3);

side_m = mean(cat(3,im_side{end-100:end}),3);

% Covert to 16-bit

top_m = uint16(top_m);

bot_m = uint16(bot_m);

sid_m = uint16(side_m);

%% Crop cone

num_edges = 5; % num of edges of polygon

% Acquire edge locations

figure, imshow(top_m,[0 14000])

edge_locs = ginput(num_edges); % Select edges for mask

xi = edge_locs(:,1); % x-coordinate locations

yi = edge_locs(:,2); % y-coordinate locations

BW_top = roipoly(top_m,xi,yi);

top_crop = uint16(BW_top);

%% Edge detection

% Get coordinates for edges

endPointsImage = bwmorph(top_crop, 'endpoints');

[rows, columns] = find(endPointsImage);

% Check to see if coordinates make sense

plot(columns,rows,'r','LineWidth',1.5);

Answer 1

[Alan]

Hi,

According to my understanding you want to use the 3 projections (top, left, and side) to construct a composite image. You could check out the warp() function to apply geometric transformations to an image and stich them later on: https://in.mathworks.com/help/images/ref/images.geotrans.warper.images.geotrans.warper.warp.html

To convert x-y coordinates of the image to polar coordinates, you could use the cart2pol() function: https://in.mathworks.com/help/matlab/ref/cart2pol.html