乘积变换实现uint16图像

function test()
 msidata = read_ENVIimage('qb_boulder_msi.img');%获取多光谱影像数据，测试数据为4个波段，bsq的存储方式，分辨率为1024*1024
 pandata=read_ENVIimage('qb_boulder_pan.img');%获取全色波段的图像矩阵，大概为4096*4096，相对于分辨率提高了4倍
 
 figure
 imshow(msidata(:,:,1:3)/500);%uint16转换到0-1
 figure
 imshow(pandata/700);%uint16转换到0-1
 Out=MultiplyFusion(msidata,pandata);
 figure
 imshow(Out(:,:,1:3));
end

%%
function [Out]=MultiplyFusion(msi,pan)
 [row1,col1,bands1]=size(msi);
 [row2,col2]=size(pan);
 rowscale=row2/row1;
 colscale=col2/col1;
 Out=zeros(row2,col2,bands1);
 msi=double(msi);
 pan=double(pan);
 %对多波段影像进行重采样
%  meanpan=MyMean(pan);
 ReampleMSI=zeros(row2,col2,bands1);
 for i=1:bands1
     ReampleMSI(:,:,i)=ReSample(msi(:,:,i),rowscale,colscale);
     Out(:,:,i)=ReampleMSI(:,:,i)/500.*pan/700;%uint16转换到0-1  %uint16转换到0-1
%      Out(:,:,i)=Normalization(Out(:,:,i));%归一化的思想并不是最好的，后面整体结果偏暗
 end
 
end
 
 
 
%%
function [B]=ReSample(A,row_scale,col_scale)
%实现图像的重采样,传入的参数是一个灰度波段的图像,k为将灰度图像采样的倍数
[row,col] = size(A);
B=zeros(row*row_scale,col*col_scale);  %这里针对实验数据是行和列增加4倍
%下面采样最邻近采样的方法进行采样，这个方法实际上针对于上采样来说就是将一个区域所有的值赋值为一样的
for i=1:row
    for j=1:col
        %上面的循环主要是找到大范围的点，下面的循环是进行具体值得赋值
        for n=0:row_scale-1 %0-3
            for m=0:col_scale-1 %0-3
                B(i*row_scale-n,j*col_scale-m)=A(i,j);
            end
        end
    end
end
%上面采样完毕之后就返回结果
end

%%
function B=Normalization(A)
%主要是针对于单波段影像进行处理
[row,col]=size(A);
minDN=1000;
maxDN=0;
for i=1:row
    for j=1:col
        if A(i,j)>maxDN
            maxDN=A(i,j);
        end
        if A(i,j)<minDN
            minDN=A(i,j);
        end
    end
end
%接下来就是将图像值进行归一化处理(0,1)
B=1/(maxDN-minDN)*(A-minDN);
end

%%
function B=MyMean(A)
%主要是针对于单波段影像进行处理
[row,col]=size(A);
sum=0;
for i=1:row
    for j=1:col
        sum=sum+A(i,j);
    end
end
%接下来就是将图像值进行归一化处理(0,1)
B=sum/row/col;
end

function data=read_ENVIimage(imagefile)
%读取.img格式的二五年间，前提是图像以',img'后缀名
if length(imagefile)>=4  %如果位于同级目录下的长度大于4才进行操作
    %主要的思路就是利用提高的img格式的文件寻找位于同级目录下的hdr文件
    switch strcmp(imagefile(length(imagefile)-3:end),'img')
        case 0
            hdrfilename=strcat(imagefile(1:(length(imagefile)-4)),'.hdr');
    end
else
    hdrfilename=strcat(imagefile,'hdr');
end
%以只读的方式打开头文件
fid = fopen(hdrfilename,'r');
info=fread(fid,'char=>char');
info=info';%转化为行向量
fclose(fid);

%查找列数
a=strfind(info,'samples = ');
b=length('samples = ');
c=strfind(info,'lines');
samples=[];
for i=a+b:c-1
    samples=[samples,info(i)];
end
samples=str2num(samples);
%查找行数
a=strfind(info,'lines   = ');
b=length('lines   = ');
c=strfind(info,'bands');
lines=[];
for i=a+b:c-1
    lines=[lines,info(i)];
end
lines=str2num(lines);
%查找波段数
a=strfind(info,'bands   = ');
b=length('bands   = ');
c=strfind(info,'header offset');
bands=[];
for i=a+b:c-1
    bands=[bands,info(i)];
end
bands=str2num(bands);
%查找数据类型
a=strfind(info,'data type = ');
b=length('data type = ');
c=strfind(info,'interleave');
datatype=[];
for i=a+b:c-1
    datatype=[datatype,info(i)];
end
datatype=str2num(datatype);
precision=[];
switch datatype
    case 1
        precision='uint8=>uint8';%头文件中datatype=1对应ENVI中数据类型为Byte，对应MATLAB中数据类型为uint8
    case 2
        precision='int16=>int16';%头文件中datatype=2对应ENVI中数据类型为Integer，对应MATLAB中数据类型为int16
    case 12
        precision='uint16=>uint16';%头文件中datatype=12对应ENVI中数据类型为Unsighed Int，对应MATLAB中数据类型为uint16
    case 3
        precision='int32=>int32';%头文件中datatype=3对应ENVI中数据类型为Long Integer，对应MATLAB中数据类型为int32
    case 13
        precision='uint32=>uint32';%头文件中datatype=13对应ENVI中数据类型为Unsighed Long，对应MATLAB中数据类型为uint32
    case 4
        precision='float32=>float32';%头文件中datatype=4对应ENVI中数据类型为Floating Point，对应MATLAB中数据类型为float32
    case 5
        precision='double=>double';%头文件中datatype=5对应ENVI中数据类型为Double Precision，对应MATLAB中数据类型为double
    otherwise
        error('invalid datatype');%除以上几种常见数据类型之外的数据类型视为无效的数据类型
end
%查找数据格式
a=strfind(info,'interleave = ');
b=length('interleave = ');
c=strfind(info,'byte order');
interleave=[];
for i=a+b:c-1
    interleave=[interleave,info(i)];
end
    interleave=strtrim(interleave);%删除字符串中前后的空格
%读取图像文件
fid = fopen(imagefile, 'r');
data = multibandread(imagefile ,[lines, samples, bands],precision,0,interleave,'ieee-le');
% data=double(data);%将影像uint16格式转化为double格式便于显示
% R=data(:,:,1);
% R=data(:,:,2);
% R=data(:,:,3);
% figure
% imagesc(R);%uint8对应double[0,1] uint16对应于double并不是[0,1]的范围，所以显示出来就是一片白
data=double(data);
% imshow(uint8(data(:,:,1)));
end