Extracting the indices in a 2D map associted with lines drawn over the map

Hello,

I have a 2D map of some geophysical variable. Visually, to mark the regions of interest and draw cross-sections, I have drawn several lines that indicate the various regions of the map. Now, I would like to extract the indices associated with each of these lines. I would then like to average the geophysical variable once I know the indices (which is simple). I just don't know how to extract the indices given these lines.

Any guidance will be appreciated!

Answer 1

[Image Analyst]

Use improfile() to draw your lines. Two of the outputs will be the vector of x and vector of y values for the line you drew.

[Sai Prasanth]

Sounds great, let me try it and update!

[Sai Prasanth]

Hello, I must admit that I am actually struggling to get this to work - the improfile doesn't work interactively when I say:

I = imread('wmax.png');

[cx,cy,c] = improfile(I)

The "Select Line Segments Interactively" section in Matlab's improfile is too vague and they don't have any example to illustrate this. Can you please elaborate further?

[Image Analyst]

@Sai Prasanth, try this -- it works fine:

grayImage = imread('moon.tif');
imshow(grayImage);
g = gcf;
g.WindowState = 'maximized'; % Maximize window.
message = sprintf('Left click once.  Then move to other end and right click');
title(message, 'FontSize', 20);
uiwait(helpdlg(message));
[cx,cy,c] = improfile()
hold on;
plot(cx, cy, 'r-', 'LineWidth', 2); % Draw where they clicked.

[Sai Prasanth]

Thank you!! This is terrific, only trouble is that I am having a hard time interpretting the cx and cy. That is - say, I have a 2D geophysical variable that I am plotting and I need the indices across one of these lines. Now, I plot the variable, save it as a png. I use the above code to read it as a png and I get a bunch of cx and cy. Now, these cx and cy and the dimensions of the grey image that the above code generates is nothing like the original 2D array. The whole point of drawing these lines was to be able to pick the array indices so that I could access the values at those specific locations to draw a histogram/average etc. How do I get across this issue?

[Image Analyst]

I don't understand what you're saying. Maybe if you say this before you draw:

axis('on', 'image');

you'll see that the cx and cy DO have coordinate values that correspond to the actual pixel location in the image. If you have some other calibration system, then you'll have to convert the pixel locations from a digital image into something else, like latitude or longitude or kilometers or whatever. To convert to kilometers - that's trivial. To get latitude and longitude, there might be something in the Mapping Toolbox that might help (I don't have it).

So getting back to your question: "The whole point of drawing these lines was to be able to pick the array indices so that I could access the values at those specific locations to draw a histogram/average etc." I believe the cx and cy are the indices. Of course they're fractional because the profile is along a straight line. If you wanted integer coordinates to line up with rows and columns, then you can simply cast to integer by rounding:

cx = round(cx); % Round to nearest integer.
cy = round(cy); % Round to nearest integer.

Now the line will be jaggy. If you want to get the intensity along the jaggy line, simply run along the line

for k = 1 : length(cx)
    intensity(k) = grayImage(cy(k), cx(k)); % Make sure y comes first, NOT x!
end

[Sai Prasanth]

Thanks for this - but I am actually looking for an index to index conversion in a sense. To give you a more complete picture, if you see the plot above, you'll see that I am plotting this 2-D variable (say, Ice-water) as a function of lat-lon. The variable, Ice-Water, as well as the lat-lon are all of the dimension 699 x 699. Since it seems like improfile only operates on images, I am converting this 2D plot into a .png image. This is where the problem starts. Now, when I use your code, the array 'grayImage' has dimensions of 875x1167x3. So the cx and cy are indeed indices but they correspond to the indices of grayImage. What I need are the indices that correspond to the original Ice-Water variable of the dimension 699x699. Please see the attachments before and after the image conversion and using your above code. Without this index conversion, I am completely stuck.

[Image Analyst]

Well you have the spatial calibration right there. So simply compute how many lat or lon degrees there are per pixel and multiply by it to convert to degrees: degrees = pixels * (degrees / pixel). You see, pixels cancel out and you're left with degrees!

[rows, columns, numberOfColorChannels] = size(grayImage) % Units of pixels
% Get the vertical spatial calibration
yl = ylim()
realWorldDistance = abs(yl(2) - yl(1)) % Units of degrees
degreesPerPixelVertically = realWorldDistance / rows
% Get the horizontal spatial calibration
xl = xlim()
realWorldDistance = abs(xl(2) - xl(1)) % Units of degrees
degreesPerPixelHorizontally = realWorldDistance / columns
% Now convert cx and cy from pixels to real world units;
cx = cx * degreesPerPixelHorizontally; % Convert to degrees
cy = cy * degreesPerPixelVertically; % Convert to degrees

[Sai Prasanth]

Thank you! :) Here, the x and y limits were just the image indices (since it's only an image, it doesn't have the lat-lon as indices), so the xlim and ylim are just index units (no units) and the degreesPerPixelHorizontally will just be 1. I can certainly load the limits from the original array though.

Edit: It is not as trivial and the index limits are messed up when there is a conversion of a plot of 2D array into an image. There are some other fringe indices outside the axis, for example.

[Sai Prasanth]

I guess, my problem would instantly be solved if there was a way to extract these indices from a line drawn on a plot - without the need for converting it to an image. Since improfile only works on images, the indices of an image are totally different.

Please see in the below code that creates a line on a simple pcolor plot - I would just like to know the indices that it passes through (without converting to an image).

f3 = figure('DefaultAxesFontSize',18, 'DefaultAxesTitleFontWeight','bold')
 set(gcf,'color','w');
 set(0,'DefaultAxesFontName','Times New Roman','DefaultTextFontName','Times New Roman','DefaultAxesFontWeight', 'bold'); 
 pcolor(wrfvars_time1.XLONG,wrfvars_time1.XLAT,log10(Int_icepath))
 xlabel('Longitude') 
 ylabel('Latitude') 
 title('log_{10}(Vertically Integrated Ice-Water) in g/m^2') 
 shading flat
 caxis([-2 5])
 xlim([-36 -25])
 ylim([6 16.15])
 colormap(map)
 colorbar; 
   
x1 = -34.49; 
x2 = -31.58; 
y1 = 11.01;
y2 = 11.67;  

plot([x1 x2], [y1 y2],'LineWidth',3,'Color','red');

[Sai Prasanth]

@Image Analyst : Hi, can you please help me out with this one? I'd like to extract the indices from a line without going through an image.