g_viz_field.m 3.83 KB
Newer Older
Pascal Bourgault's avatar
Pascal Bourgault committed
1
2
3
4
5
6
7
8
9
10
11
function [h_figure,h_pcolor,h_datetext] = g_viz_field(imgFname,rectFile,varargin)
%G_VIZ_FIELD Generates a map with a georectified image
%       G_VIZ_FIELD Generates a map with the georectified image in imgFname
%                   Converts in grayscale if needed.
%       Inputs:
%           imgFname,   filename of the image to georectify
%           rectFile,   .mat file created by g_rect
%       Parameters (name, value)
%           showTime,   Default 0, if 1, displays the timestamp of the
%                       image in the figure's title
%           showLand,   Default 0, if 'f' or 'h', displays the land contour
12
%                       with m_gshhs_f or m_gshhs_h
Pascal Bourgault's avatar
Pascal Bourgault committed
13
14
15
16
17
18
19
20
%           landcolor,  Default [241 235 144]/255, Color of the land on the
%                       map
%
%       Outputs:
%           h_figure, h_pcolor and h_datetext are the handles to the
%               corresponding objects in the figure, for use in g_viz_anim
%

21
22
23
24
25
26
% Set some plotting parameters
ms = 10;  % Marker Size
fs = 10;  % Font size
lw = 2;  % Line width

% Option
Pascal Bourgault's avatar
Pascal Bourgault committed
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
show_time = 0;
show_land = 0;
land_color = [241 235 144]/255;
if length(varargin) > 1
    for i=1:2:length(varargin)
        switch lower(varargin{i})
            case 'showtime'
                show_time = varargin{i+1};
            case 'showland'
                show_land = varargin{i+1};
            case 'landcolor'
                land_color = varargin{i+1};
        end
    end
end

43
% Load the georectification file
Pascal Bourgault's avatar
Pascal Bourgault committed
44
load(rectFile);
Daniel Bourgault's avatar
Daniel Bourgault committed
45

46
47
% Determine the region to plot, delimited by the GCP and the camera
% position. Add a factor (fac) all around
Daniel Bourgault's avatar
Daniel Bourgault committed
48
fac = 0.1;
49
50
51
52
53
54
55
56
57
lon_min = min([lon_gcp,LON0]);
lon_max = max([lon_gcp,LON0]);
lat_min = min([lat_gcp,LAT0]);
lat_max = max([lat_gcp,LAT0]);

lon_min = lon_min - fac*abs(lon_max - lon_min);
lon_max = lon_max + fac*abs(lon_max - lon_min);
lat_min = lat_min - fac*abs(lat_max - lat_min);
lat_max = lat_max + fac*abs(lat_max - lat_min);
Daniel Bourgault's avatar
Daniel Bourgault committed
58
59
60

rgb0 = double(imread(imgFname))/255;

Pascal Bourgault's avatar
Pascal Bourgault committed
61
pp = size(rgb0,3);
Daniel Bourgault's avatar
Daniel Bourgault committed
62
if pp == 3
Pascal Bourgault's avatar
Pascal Bourgault committed
63
  int = g_rgb2gray(rgb0); 
Daniel Bourgault's avatar
Daniel Bourgault committed
64
65
66
67
68
69
else
  int = rgb0;     
end
clear rgb0;
int = int';

70
71
% Normalize the image if wanted
%int = int - nanmean(nanmean(int));
Daniel Bourgault's avatar
Daniel Bourgault committed
72

73
74
%figure('Renderer', 'painters', 'Position', [100 100 900 700])
figure
Daniel Bourgault's avatar
Daniel Bourgault committed
75
76
77
m_proj('mercator','longitudes',[lon_min lon_max],'latitudes',[lat_min lat_max]);
hold on;

Pascal Bourgault's avatar
Pascal Bourgault committed
78
79
cmap = contrast(int,256);
colormap(cmap);
80
h = m_pcolor(LON,LAT,int);
Daniel Bourgault's avatar
Daniel Bourgault committed
81
82
shading('interp');

83
% Uncomment one of these lines if you want the coastline to be plotted.
Pascal Bourgault's avatar
Pascal Bourgault committed
84
if show_land
85
86
87
    if strcmpi(show_land,'fjord')
        m_usercoast('fjord_coastlines.mat', 'patch', land_color)
    elseif strcmpi(show_land, 'f')
Pascal Bourgault's avatar
Pascal Bourgault committed
88
89
90
91
92
        m_gshhs_f('patch',land_color)
    else
        m_gshhs_h('patch',land_color)
    end
end
Daniel Bourgault's avatar
Daniel Bourgault committed
93

Pascal Bourgault's avatar
Pascal Bourgault committed
94
95
96
97
98
99
100
101
102
103
104
if show_time
    info = imfinfo(imgFname);
    if isField(info,'DateTime')
        date = info.DateTime;
    elseif isField(info,'Comment')
        date = info.Comment;
    else
        date = '';
    end
    ht = title(date);
end
Daniel Bourgault's avatar
Daniel Bourgault committed
105
106
107
108
m_plot(LON0,LAT0,'kx','markersize',ms,'linewidth',lw);  % Camera location

%% Plot GCPs and ICPs.
for n=1:length(i_gcp)
109
110
111
112
113
    
  % Plot the original GCPs that may have elevations
  m_plot(lon_gcp0(n),lat_gcp0(n),'ko','markersize',ms,'linewidth',lw);
  
  % Plot the projected GCPs that are actually used for the minimization
Daniel Bourgault's avatar
Daniel Bourgault committed
114
  m_plot(lon_gcp(n),lat_gcp(n),'bo','markersize',ms,'linewidth',lw);
115
116
  
  % Plot the Image Control Points once georectified
Daniel Bourgault's avatar
Daniel Bourgault committed
117
118
119
120
121
122
123
124
  m_plot(LON(i_gcp(n),j_gcp(n)),LAT(i_gcp(n),j_gcp(n)),'rx','MarkerSize',ms,'linewidth',lw);
end

%title([datestr(mtime,31),' UTC']);
%title([time,' UTC']);

clear LON LAT X Y

125
m_grid('box','fancy','fontsize',fs);
126
127
%m_gshhs('f','patch', 'gray');
%print('g_rect_drone_1', '-dpng')
Pascal Bourgault's avatar
Pascal Bourgault committed
128
129
130
131
132
133
134
if nargout > 1
    h_figure = gcf;
    h_pcolor = h;
    if nargout == 3
        h_datetext = ht;
    end
end