source: issm/branches/trunk-jpl-damage/src/m/model/plot/plot_overlay.m@ 11503

function plot_overlay(md,data,options,plotlines,plotcols,i)
%PLOT_OVERLAY - superimpose radar image to a given field
%
%   Usage:
%      plot_overlay(md,options,plotlines,plotcols,i)
%
%   See also: PLOTMODEL

%process mesh and data
[x y z elements is2d isplanet]=processmesh(md,[],options);
if strcmpi(data,'none'),
        radaronly=1;
        data=NaN*ones(md.mesh.numberofvertices,1);
        datatype=1;
else
        radaronly=0;
        [data datatype]=processdata(md,data,options);
end

%check is2d
if ~is2d, 
        error('buildoverlay error message: overlay not supported for 3d meshes, project on a layer');
end
if datatype==3,
        error('buildoverlay error message: overlay not supported for quiver plots');
end


%radar power
if ~any(isnan(md.radaroverlay.x)) & ~any(isnan(md.radaroverlay.y)) & ~any(isnan(md.radaroverlay.pwr)),
        disp('plot_overlay info: the radar image held by the model is being used');
        xlim=[min(md.radaroverlay.x) max(md.radaroverlay.x)];
        ylim=[min(md.radaroverlay.y) max(md.radaroverlay.y)];
else
        disp('Extracting radar image...');
        %Get xlim and ylim (used to extract radar image)
        xlim=getfieldvalue(options,'xlim',[min(x) max(x)])/getfieldvalue(options,'unit',1);
        ylim=getfieldvalue(options,'ylim',[min(y) max(y)])/getfieldvalue(options,'unit',1);
        options=addfielddefault(options,'xlim',xlim);
        options=addfielddefault(options,'ylim',ylim);
        md=radarpower(md,options);
end

%InterpFromMeshToGrid
xmin=min(md.radaroverlay.x);
ymax=max(md.radaroverlay.y);
xspacing=(max(md.radaroverlay.x)-min(md.radaroverlay.x))/(length(md.radaroverlay.x));
yspacing=(max(md.radaroverlay.y)-min(md.radaroverlay.y))/(length(md.radaroverlay.y));
nlines=length(md.radaroverlay.y);
ncols =length(md.radaroverlay.x);
disp('Interpolating data on grid...');
[x_m y_m data_grid]=InterpFromMeshToGrid(elements,x/getfieldvalue(options,'unit',1),y/getfieldvalue(options,'unit',1),...
        data,xmin,ymax,xspacing,yspacing,nlines,ncols,NaN);

%Process data_grid
pos=find(isinf(data_grid));
if ~isempty(pos),
        disp('Warning: removing Infs from vector (probably log(0)?)');
        data_grid(pos)=NaN;
end
if exist(options,'caxis'),
        caxis_opt=getfieldvalue(options,'caxis');
        data_grid(find(data_grid<caxis_opt(1)))=caxis_opt(1);
        data_grid(find(data_grid>caxis_opt(2)))=caxis_opt(2);
        data_min=caxis_opt(1);
        data_max=caxis_opt(2);
else
        data_min=min(data_grid(:));
        data_max=max(data_grid(:));
end
data_nan=find(isnan(data_grid));

%Generate HSV image
contrast=getfieldvalue(options,'contrast',1);  
transparency=getfieldvalue(options,'alpha',1);
data_grid(data_nan)=data_min; 

colorm=getfieldvalue(options,'colormap','Rignot');
if strcmpi(colorm,'Rignot'),
        %hue (H)
        h_data=(data_grid-data_min)/(data_max-data_min+eps);
        if radaronly,
                h_data(:)=0;
        end
        %saturation (S)
        s_data=max(min((0.1+h_data).^(1/transparency),1),0);
elseif strcmpi(colorm,'Seroussi'),
        %hue (H)
        h_data=1-(data_grid-data_min)/(data_max-data_min+eps)*0.7;
        %h_data=(data_grid-data_min)/(data_max-data_min)*2/3;
        if radaronly,
                h_data(:)=0;
        end
        %saturation (S)
        s_data=max(min((0.1+h_data).^(1/transparency),1),0);
elseif strcmpi(colorm,'redblue')
        data_mean=data_min+(data_max-data_min)/2;
        %hue (H)
        %h_data=0.7*ones(size(data_grid));
        %h_data(find(data_grid>data_mean))=1;
        h_data=1*ones(size(data_grid));
        h_data(find(data_grid>data_mean))=0.7;
        %saturation (S)
        s_data=max(min(abs(data_grid-data_mean)/(data_max-data_mean) ,1),0);
else
        error('colormap not supported yet. (''Rignot'' and ''redblue'' are the only cupported colormaps)');
end

%Saturation is 0 in NaNs
s_data(data_nan)=0;
%intensity (V)
radar=(md.radaroverlay.pwr).^(contrast);
v_data=radar/max(radar(:)); %use radar power as intensity

%Transform HSV to RGB
image_hsv=zeros(size(data_grid,1),size(data_grid,2),3);
image_hsv(:,:,1)=h_data;
image_hsv(:,:,2)=s_data;
image_hsv(:,:,3)=v_data;
image_rgb=hsv2rgb(image_hsv);

%Select plot area 
subplot(plotlines,plotcols,i);

%Plot: 
imagesc(md.radaroverlay.x*getfieldvalue(options,'unit',1),md.radaroverlay.y*getfieldvalue(options,'unit',1),image_rgb);set(gca,'YDir','normal');

%last step: mesh overlay?
if exist(options,'edgecolor'),
        hold on
        A=elements(:,1); B=elements(:,2); C=elements(:,3); 
        patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',zeros(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor'));
end

%Apply options, without colorbar and without grid
options=changefieldvalue(options,'colormap',colorm);           %We used an HSV colorbar
if ~isnan(data_min),
        options=changefieldvalue(options,'caxis',[data_min data_max]); %force caxis so that the colorbar is ready
end
options=addfielddefault(options,'axis','equal off');           %default axis
applyoptions(md,data,options);