brdf_bate1.m

BRDF model,BRDF model,BRDF model,BRDF model,BRDF model,BRDF model

% 求某一波段的二向反射率
% file_time 时间
% Band      波段号
% geom      角度信息 [geom_szen;geom_vzen;geom_relaz]
% --------------------------------------------------------
% =====          带有路径的 求BRDF     ===================
%---------------------------------------------------------
function Ref = Brdf_bate1(path,file_time,Band,geom)
% Ref = BrdfBand(2008257,6,geom);
[data_iso,data_vol,data_geo,rows,columns] = LayerData(path,file_time,Band);
%title_str = Titlestr(file_time,Band,geom)
[k1,k2] = brdf_forward(geom);
Ref = data_iso*0.001 + data_vol*0.001*k1 + data_geo*0.001*k2;
[m,n]=size(Ref);
for i=1:m
    for j=1:n
       if Ref(i,j)<0|Ref(i,j)>1;
            %Ref(i,j)=(Ref(i-1,j-1)+ Ref(i-1,j)+Ref(i-1,j+1)+ Ref(i,j-1)+ Ref(i,j)+ Ref(i,j+1)+ Ref(i+1,j-1)+ Ref(i+1,j)+ Ref(i+1,j+1))/9;
     Ref(i,j)=0;
        end
   end
end



function [data_iso,data_vol,data_geo,rows,columns] = LayerData(path,file_time,Band)
%a = Band;
path1 = path;
path2 = 'MCD43A1_';
path3 = int2str(file_time);
path4 = '.hdf';
hdf_path = [path1,path2,path3,path4];
%hdf_path = 'D:\MODIS\MCD43A1\MCD43A1.A2008153.hdf';
fileinfo = hdfinfo(hdf_path,'eos');
 rows = fileinfo.Grid.Rows;
 columns = fileinfo.Grid.Columns;
if Band>0 & Band<=11
    c=[1:3:28;2:3:29;3:3:30];
    
    Layer_Name1 = fileinfo.Grid.DataFields(c(1,Band),1).Name;
    Layer_Name2 = fileinfo.Grid.DataFields(c(2,Band),1).Name;
    Layer_Name3 = fileinfo.Grid.DataFields(c(3,Band),1).Name;
    
%data_iso_int16 = hdfread(hdf_path ,Layer_Name1, 'Index', {[1  1],[1  1],[570  481]});
%data_vol_int16 = hdfread(hdf_path ,Layer_Name2, 'Index', {[1  1],[1  1],[570  481]});
%data_geo_int16 = hdfread(hdf_path ,Layer_Name3, 'Index', {[1  1],[1  1],[570  481]});
data_iso_int16 = hdfread(hdf_path ,Layer_Name1, 'Index', {[1  1],[1  1],[rows,columns]});
data_vol_int16 = hdfread(hdf_path ,Layer_Name2, 'Index', {[1  1],[1  1],[rows,columns]});
data_geo_int16 = hdfread(hdf_path ,Layer_Name3, 'Index', {[1  1],[1  1],[rows,columns]});
  data_iso = double(data_iso_int16);          
  data_vol = double(data_vol_int16);           
  data_geo = double(data_geo_int16); 

disp( Layer_Name1);
disp( Layer_Name2);
disp( Layer_Name3);
%readed = ' 数据已读取！';
%dis_v = [Band, readed];
%disp(dis_v);
else disp('Eorr!');
end
function [rosskernel,likernel] = brdf_forward(geom)

 %  this function runs the RossThickLiSparseReciprocal model in the 
 %  forward mode and returns the calculated reflectance.
 %
 %  Parameters:
 %  geom_t geom                   view geometry
 %  param_t params                BRDF kernel weights parameters 
 %

%  float phi, cosphi, sinphi, tv, costv, sintv, ti, costi, sinti;
%  float cosphaang, phaang, sinphaang, rosskernel, tantv, tanti;
%  float likernel, refl;


M_PI = 3.1415926535;
M_PI_2 =  M_PI/2;
M_PI_4 = M_PI/4;
BR = 1.0;                
HB = 2.0;
geom_szen = geom(1,1);
geom_vzen = geom(2,1);
geom_relaz = geom(3,1);
  ti = geom_szen;
  tv = geom_vzen;
  phi = geom_relaz;

  ti = D2R(ti); 
  tv = D2R(tv); 
  phi = D2R(phi); 
        
  cosphi = cos(phi);
  sinphi = sin(phi);
  costv = cos(tv);
  costi = cos(ti);
  sintv = sin(tv);
  sinti = sin(ti);
        
  
[cosphaang,phaang,sinphaang] = GetPhaang(costv,costi,sintv,sinti,cosphi);     
  rosskernel = ((M_PI_2 - phaang) * cosphaang + sinphaang)./(costi + costv)-M_PI_4;
        
  tantv = sintv./costv;
  tanti = sinti./costi;
        
 likernel = LiKernel(HB,BR,tantv,tanti,sinphi,cosphi);
 
function  [cosres,res,sinres] = GetPhaang(cos1,cos2,sin1,sin2,cos3)

  %This func calculates the phase angle 

  cosres = cos1*cos2 + sin1*sin2*cos3;
  res = acos( max(-1., min(1.,cosres)) );
  sinres = sin(res);
function y=D2R(p)
     M_PI = 3.1415926535;
   y = M_PI*p/180;
function [res] = GetDistance(tan1,tan2,cos3)

% This func calculates the D distance */

  float temp;

  temp = tan1*tan1+tan2*tan2-2.*tan1*tan2*cos3;
  res = sqrt(max(0.,temp));
function [overlap,temp] = GetOverlap(hbratio,distance,cos1,cos2,tan1,tan2,sin3)

  %This func calculates the Overlap distance */
  M_PI = 3.1415926535;
  temp = 1./cos1 + 1./cos2;
   cost = hbratio*sqrt(distance*distance+tan1*tan1*tan2*tan2*sin3*sin3)./temp;
   cost = max(-1., min(1.,cost));
   tvar = acos(cost);
   sint = sin(tvar);
   overlap = 1./M_PI *(tvar-sint*cost)*(temp);
   overlap = max(0.,overlap);

function  [sinp,cosp,tanp] = GetpAngles(brratio,tan1)
%This func calculates the 'prime' angles

  float angp;

    tanp = brratio*tan1;
  if(tanp < 0);
      tanp = 0.;
  end
  angp = atan(tanp);
  sinp = sin(angp);
  cosp = cos(angp);
  
function result = LiKernel(hbratio,brratio,tantv,tanti,sinphi,cosphi)


  %This func calculates the LiSparseModis kernel */

 % float sintvp, costvp, tantvp, sintip, costip, tantip;
 % float phaangp, cosphaangp, sinphaangp, distancep, overlap, temp;

 [sintvp,costvp,tantvp] = GetpAngles(brratio,tantv);
 [sintip,costip,tantip] = GetpAngles(brratio,tanti);
 [cosphaangp,phaangp,sinphaangp] = GetPhaang(costvp,costip,sintvp,sintip,cosphi);
  distancep = GetDistance(tantvp,tantip,cosphi);
 [overlap,temp] = GetOverlap(hbratio,distancep,costvp,costip,tantvp,tantip,sinphi);

  result = overlap - temp + 0.5 * (1.+cosphaangp)./costvp./costip;
function  title_str = Titlestr(file_time,Band,geom)
time_hdf = int2str(file_time);
yy =  time_hdf(1:4);
mmdd = time_hdf(5:7);
data_hdf = str2num(mmdd);
mm = datestr(data_hdf,5);
dd = datestr(data_hdf,7);
time_str = [yy '年' mm '月' dd '日'];
geom_szen = int2str(geom(1,1));

geom_vzen = int2str(geom(2,1));
geom_relaz =int2str(geom(3,1));
geom_str = [geom_szen '° ' geom_vzen '° ' geom_relaz '° '];
band_str = int2str(Band);
title_str = [time_str ' ' geom_str ' ' 'BAND : ' band_str];