128-bit long double floating-point data type

本文深入探讨了IBM AIX操作系统中128位长双精度浮点数的数据类型,该类型提供比默认64位长双精度更高的精度,能处理多达31位有效数字。文章详细讲解了使用128位长双精度浮点数编译程序的方法,其与IEEE754标准的符合性,实现格式,以及数值宏值的特殊考虑。

https://www.ibm.com/support/knowledgecenter/en/ssw_aix_71/com.ibm.aix.genprogc/128bit_long_double_floating-point_datatype.htm

The AIX® operating system supports a 128-bit long double data type that provides greater precision than the default 64-bit long double data type. The 128-bit data type can handle up to 31 significant digits (compared to 17 handled by the 64-bit long double). However, while this data type can store numbers with more precision than the 64-bit data type, it does not store numbers of greater magnitude.

The following special issues apply to the use of the 128-bit long double data type:

  • Compiling programs that use the 128-bit long double data type
  • Compliance with the IEEE 754 standard
  • Implementing the 128-bit long double format
  • Values of numeric macros

Compiling programs that use the 128-bit long double data type

To compile C programs that use the 128-bit long double data type, use the xlc128 command. This command is an alias to the xlc command with support for the 128-bit data type. The xlc command supports only the 64-bit long double data type.

The standard C library, libc.a, provides replacements for libc.a routines which are implicitly sensitive to the size of a long double. Link with the libc.a library when compiling applications that use the 64-bit long double data type. Link applications that use 128-bit long double values with both the libc128.a and libc.a libraries. When linking, be sure to specify the libc128.a library before the libc.a library in the library search order.

Compliance with IEEE 754 standard

The 64-bit implementation of the long double data type is fully compliant with the IEEE 754 standard, but the 128-bit implementation is not. Use the 64-bit implementation in applications that must conform to the IEEE 754 standard.

The 128-bit implementation differs from the IEEE standard for long double in the following ways:

  • Supports only round-to-nearest mode. If the application changes the rounding mode, results are undefined.
  • Does not fully support the IEEE special numbers NaN and INF.
  • Does not support IEEE status flags for overflow, underflow, and other conditions. These flags have no meaning for the 128-bit long double inplementation.
  • The 128-bit long double data type does not support the following math APIs: atanhlcbrtlcopysignlexp2l,expm1lfdimlfmalfmaxlfminlhypotlilogblllrintlllroundllog1pllog2llogbllrintllroundlnanl,nearbyintlnextafterlnexttowardnexttowardfnexttowardlremainderlremquolrintlroundlscalblnl,scalbnltgammal, and truncl.

Implementing the 128- bit long double format

A 128-bit long double number consists of an ordered pair of 64-bit double-precision numbers. The first member of the ordered pair contains the high-order part of the number, and the second member contains the low-order part. The value of the long double quantity is the sum of the two 64-bit numbers.

Each of the two 64-bit numbers is itself a double-precision floating-point number with a sign, exponent, and significand. Typically the low-order member has a magnitude that is less than 0.5 units in the last place of the high part, so the values of the two 64-bit numbers do not overlap and the entire significand of the low-order number adds precision beyond the high-order number.

This representation results in several issues that must be considered in the use of these numbers:

  • The exponent range is the same as that of double precision. Although the precision is greater, the magnitude of representable numbers is the same as 64-bit double precision.
  • As the absolute value of the magnitude decreases (near the denormal range), the additional precision available in the low-order part also decreases. When the value to be represented is in the denormal range, this representation provides no more precision than the 64-bit double-precision data type.
  • The actual number of bits of precision can vary. If the low-order part is much less than 1 ULP of the high-order part, significant bits (either all 0's or all 1's) are implied between the significant of the high-order and low-order numbers. Certain algorithms that rely on having a fixed number of bits in the significand can fail when using 128-bit long double numbers.

Values of numeric macros

Because of the storage method for the long double data type, more than one number can satisfy certain values that are available as macros.The representation of 128-bit long double numbers means that the following macros required by standard C in the values.h file do not have clear meaning:

  • Number of bits in the mantissa (LDBL_MANT_DIG)
  • Epsilon (LBDL_EPSILON)
  • Maximum representable finite value (LDBL_MAX)

Number of bits in the mantissa

The number of bits in the significand is not fixed, but for a correctly formatted number (except in the denormal range) the minimum number available is 106. Therefore, the value of the LDBL_MANT_DIG macro is 106.

Epsilon

The ANSI C standard defines the value of epsilon as the difference between 1.0 and the least representable value greater than 1.0, that is, b**(1-p), where b is the radix (2) and p is the number of base b digits in the number. This definition requires that the number of base b digits is fixed, which is not true for 128-bit long double numbers.

The smallest representable value greater than 1.0 is this number:

0x3FF0000000000000, 0x0000000000000001

The difference between this value and 1.0 is this number:

0x0000000000000001, 0x0000000000000000
0.4940656458412465441765687928682213E-323

Because 128-bit numbers usually provide at least 106 bits of precision, an appropriate minimum value for p is 106. Thus, b**(1-p) and 2**(-105) yield this value:

0x3960000000000000, 0x0000000000000000
0.24651903288156618919116517665087070E-31

Both values satisfy the definition of epsilon according to standard C. The long double subroutines use the second value because it better characterizes the accuracy provided by the 128-bit implementation.

Maximum long double value

The value of the LDBL_MAX macro is the largest 128-bit long double number that can be multiplied by 1.0 and yield the original number. This value is also the largest finite value that can be generated by primitive operations, such as multiplication and division:

0x7FEFFFFFFFFFFFFF, 0x7C9FFFFFFFFFFFFF
0.1797693134862315907729305190789002575e+309
1. Image raw data file format (MRD) Imaging data is stored in a file with an extension .MRD. The file format consists of: 1) a 256 byte header 2) a 256 byte text description 3) the data proper 4) an ASCII string containing the sample information filename 5) an ASCII copy of the parameter file (.PPR file) used 256 byte header Bytes (hex) 'C' data type Usage 0-3 Long (4 byte integer) Dimension 1, number of data samples 4-7 Long Dimension 2, number of views (2D data) 8-B Long Dimension 3, number of secondary views (3D data) C-F Long Dimension 4, number of slices ... Unspecified 12-13 Int (2 byte integer) Data type code ... Unspecified 98-9B Long Dimension 5, number of echoes 9C-9F Long Dimension 6, number of experiments The data type code indicates the format of individual data elements, as follows: 0x00 unsigned char 1 byte 0x01 signed char 1 byte 0x02 short 2 bytes 0x03 int 2 bytes 0x04 long 4 bytes 0x05 float 4 bytes 0x06 double 8 bytes 0x10 bit mask indicating complex data as real, imaginary pairs of base data type All MRS .MRD files are at present signed char, int or floating point, either real or complex. The default is complex floating point, data type 0x15. 256 byte text Bytes (hex) 100-1FF 'C' data type 256 bytes ASCII This is normally unused Data Bytes (hex) 200 --> 'C' data type Usage Text, 0 terminated string Complex or single value Signed char, int or float Usage Data Multi-dimensional data are stored sequentially, the index of dimension 1 (samples) varying most rapidly, followed by the other dimensions in order i.e. samples (1DFT), views (2DFT), secondary views (3DFT), slices, echoes, and experiments. Sample file name This is the fully qualified path name of the sample information file, as a zero terminated ASCII string, currently padded to 120 bytes. Parameters The contents of the .PPR parameter file. This is ASCII text, in lines delimited by <CR> <LF> (0x0d, 0x0a). The format is described below. The current values of the instrument parameters are also appended to the PPR section of the MRD file in the format :_<win.ini parameter identifier> <value> <CR> <LF> :END <CR> <LF> 2. Image files MRS reconstructed image files are in a format very similar to the first part of that used for raw data (.MRD). They consist of: 1) a 256 byte header 2) a 256 byte text description 3) the data proper 256 byte header Bytes (hex) 'C' data type Usage 0-3 Long (4 byte integer) Number of pixels in X direction 4-7 Long Number of pixels in Y direction 8-B Long Number of pixels in Z direction (currently 1) C-F Long Dimension 4 (currently 1) ... Unspecified 12-13 Int (2 byte integer) Data type code ... Unspecified 30-33 Float (4 bytes) Scaling factor from highest pixel to 4095 ... Unspecified 6D Unsigned char (1 byte) Number of bits per pixel (currently 12) ... Unspecified The data type code indicates the format of individual data elements, as follows: 0x00 unsigned char 1 byte 0x01 signed char 1 byte 0x02 short 2 bytes 0x03 int 2 bytes 0x04 long 4 bytes 0x05 float 4 bytes 0x06 double 8 bytes 0x10 bit mask indicating complex data as real, imaginary pairs of base data type All MRS .SUR files are at present integer, indicated by data type 0x03. 256 text description At present this is unused and contains the ASCII text "There is no script !". Data Data follows, at present of integer type. The pixel positions increments in the X dimension until the next increment in the Y dimension.请根据这些信息读取MRD文件,显示出图像
07-04
;------------------------------------------------------------------------------------------------- ; REMOVE THICK CLOUD IN SATELLITE IMAGES ; ; This code can be used for whole TM scene ; Developed and coded by Zhu Xiaolin ; Department of Geography,the Ohio State University ; Email:zhuxiaolin55@gmail.com ; updated:2017-11-21 ; Debug history: ; 2013-10-23 correct:1) the bug of similar pixel searching ; 2) the bug of using wrong values in cloudy image ; 2017-11-21:correct a nonstop loop bug. It happens when processing large cloud patches and ; no similar pixel selected for a cloudy pixel ; 2021-05-09:correct a bug of similar pixel searching for line-shape clouds ; ; Please cite: Zhu, X., Gao, F., Liu, D. and Chen, J. A modified ; neighborhood similar pixel interpolator approach for removing ; thick clouds in Landsat images. IEEE Geoscience and Remote ; Sensing Letters, 2012, 9(3), 521-525 ; NOTE: the efficiency may be low for large image.If so, please use CLOUD_REMOVE_FAST, ; in which some process was simplified, but the accuracy is still satisfactory ;------------------------------------------------------------------------------------------------- ;function of open file Pro GetData,ImgData = ImgData,ns = ns,nl = nl,nb = nb,Data_Type = Data_Type,$ FileName = FileName,Map_info = map_Info, Fid = Fid Filter = ['all file;*.*'] Envi_Open_File,FileName,R_Fid = Fid Envi_File_Query,Fid,ns = ns,nl = nl,nb = nb,Data_Type = Data_Type map_info = envi_get_map_info(fid=Fid) dims = [-1,0,ns - 1 ,0,nl - 1] case Data_Type Of 1:ImgData = BytArr(ns,nl,nb) ; BYTE Byte 2:ImgData = IntArr(ns,nl,nb) ; INT Integer 3:ImgData = LonArr(ns,nl,nb) ; LONG Longword integer 4:ImgData = FltArr(ns,nl,nb) ; FLOAT Floating point 5:ImgData = DblArr(ns,nl,nb) ; DOUBLE Double-precision floating 6:ImgData = COMPLEXARR(ns,nl,nb); complex, single-precision, floating-point 9:ImgData = DCOMPLEXARR(ns,nl,nb);complex, double-precision, floating-point 12:ImgData = UINTARR(ns,nl,nb) ; unsigned integer vector or array 13:ImgData = ULONARR(ns,nl,nb) ; unsigned longword integer vector or array 14:ImgData = LON64ARR(ns,nl,nb) ;a 64-bit integer vector or array 15:ImgData = ULON64ARR(ns,nl,nb) ;an unsigned 64-bit integer vector or array EndCase For i = 0,nb-1 Do Begin Dt = Envi_Get_Data(Fid = Fid,dims = dims,pos=i) ImgData[*,*,i] = Dt[*,*] EndFor End ;------------------------------------------------------------------- ; main body of the program ;------------------------------------------------------------------- pro CLOUD_REMOVE t0=systime(1) ;set parameters ;------------------------------------------------ num_class=4 ;set the estimated number of classes in image min_pixel=20 ;set the sample size of similar pixels extent1=1 ;set the range of cloud neighborhood DN_min=0.0 ;set the range of DN DN_max=1.0 patch_long=500 ;set the block size when process large images, 500~1000 recommend ;------------------------------------------------- FileName1 = Dialog_PickFile(title = 'Open the cloudy image:') envi_open_file,FileName1,r_fid=fid envi_file_query,fid,ns=ns,nl=nl,nb=nb,dims=dims map_info = envi_get_map_info(fid=fid) orig_ns=ns orig_nl=nl n_ns=ceil(float(ns)/patch_long) n_nl=ceil(float(nl)/patch_long) ind_patch=intarr(4,n_ns*n_nl) for i_ns=0,n_ns-1,1 do begin for i_nl=0,n_nl-1,1 do begin ind_patch[0,n_ns*i_nl+i_ns]=i_ns*patch_long ind_patch[1,n_ns*i_nl+i_ns]=min([ns-1,(i_ns+1)*patch_long-1]) ind_patch[2,n_ns*i_nl+i_ns]=i_nl*patch_long ind_patch[3,n_ns*i_nl+i_ns]=min([nl-1,(i_nl+1)*patch_long-1]) endfor endfor temp_file=FileName1 tempoutname=temp_file+'_temp_cloud' pos=indgen(nb) for isub=0,n_ns*n_nl-1,1 do begin dims=[-1,ind_patch[0,isub],ind_patch[1,isub],ind_patch[2,isub],ind_patch[3,isub]] envi_doit, 'resize_doit', fid=fid, pos=pos, dims=dims, interp=0, rfact=[1,1], $ out_name=tempoutname+strtrim(isub+1,1), r_fid=r_fid1 envi_file_mng, id=r_fid1, /remove endfor envi_file_mng, id=fid, /remove ;----------------------------------------------------------- FileName2 = Dialog_PickFile(title = 'Open the input clear image:') envi_open_file,FileName2,r_fid=fid tempoutname=temp_file+'_temp_clear' pos=indgen(nb) for isub=0,n_ns*n_nl-1,1 do begin dims=[-1,ind_patch[0,isub],ind_patch[1,isub],ind_patch[2,isub],ind_patch[3,isub]] envi_doit, 'resize_doit', fid=fid, pos=pos, dims=dims, interp=0, rfact=[1,1], $ out_name=tempoutname+strtrim(isub+1,1), r_fid=r_fid1 envi_file_mng, id=r_fid1, /remove endfor envi_file_mng, id=fid, /remove ;------------------------------------------------------------ FileName3 = Dialog_PickFile(title = 'Open the cloud mask image:') envi_open_file,FileName3,r_fid=fid tempoutname=temp_file+'_temp_cloud_mask' pos=indgen(1) for isub=0,n_ns*n_nl-1,1 do begin dims=[-1,ind_patch[0,isub],ind_patch[1,isub],ind_patch[2,isub],ind_patch[3,isub]] envi_doit, 'resize_doit', fid=fid, pos=pos, dims=dims, interp=0, rfact=[1,1], $ out_name=tempoutname+strtrim(isub+1,1), r_fid=r_fid1 envi_file_mng, id=r_fid1, /remove endfor envi_file_mng, id=fid, /remove ;------------------------------------------------------------------- tempoutname=temp_file+'_temp_filled' print,'there are total',n_ns*n_nl,'patches' for isub=0,n_ns*n_nl-1,1 do begin ;open each block FileName = temp_file+'_temp_cloud' GetData,ImgData=ImgData,ns = ns1,nl = nl1,nb = nb,Data_Type = Data_Type,FileName = FileName+strtrim(isub+1,1),Fid = Fid1 fine1=ImgData FileName = temp_file+'_temp_clear' GetData,ImgData=ImgData,ns = ns2,nl = nl2,nb = nb,FileName = FileName+strtrim(isub+1,1),Fid = Fid2 fine2=ImgData FileName = temp_file+'_temp_cloud_mask' GetData,ImgData=ImgData,FileName = FileName+strtrim(isub+1,1),Fid = Fid3 cloud=ImgData ;----------------------------------------------------------------------- cloud_posite=where(cloud ne 0,c_c) if (c_c gt 0) then begin num_cloud_path=max(cloud[cloud_posite]) ; for icloud=1,num_cloud_path,1 do begin ; remove each cloud patches cloud_area=where(cloud eq icloud,num_cloud1) if (num_cloud1 gt 0) then begin cloud_pixel_locate=intarr(2,num_cloud1) cloud_pixel_locate[1,*]=fix(cloud_area/ns1) cloud_pixel_locate[0,*]=cloud_area-ns1*cloud_pixel_locate[1,*] left_cloud=min(cloud_pixel_locate[0,*]) ;find out the cloud area right_cloud=max(cloud_pixel_locate[0,*]) up_cloud=min(cloud_pixel_locate[1,*]) down_cloud=max(cloud_pixel_locate[1,*]) cloud_area=0 ;clear it cloud_pixel_locate=0 ;clear it ;neigborhood of cloud area: rectangle which can completely cover cloud left_region=max([0,left_cloud-extent1]) right_region=min([ns1-1,right_cloud+extent1]) up_region=max([0,up_cloud-extent1]) down_region=min([nl1-1,down_cloud+extent1]) a_region=right_region-left_region+1 b_region=down_region-up_region+1 x_center=a_region/2.0 ; the position of cloud center y_center=b_region/2.0 sub_fine1=fine1[left_region:right_region,up_region:down_region,*] sub_fine2=fine2[left_region:right_region,up_region:down_region,*] sub_cloud=cloud[left_region:right_region,up_region:down_region] ;---------------------remove background----------------- fine1_value=sub_cloud fine1_value[*,*]=0 fine2_value=fine1_value for ib=0,nb-1,1 do begin fine1_value=mean(fine1_value+sub_fine1[*,*,ib]) fine2_value=mean(fine2_value+sub_fine2[*,*,ib]) endfor ind_com=where(fine1_value eq 0 or fine2_value eq 0, c_com) if (c_com gt 0) then begin sub_cloud[ind_com]=-1 ;mask out background endif fine1_value=0 fine2_value=0 similar_th_band=fltarr(nb) for iband=0,nb-1,1 do begin similar_th_band[iband]=stddev(sub_fine2[*,*,iband])*2.0/float(num_class) ;compute the threshold of similar pixel endfor ind_cloud=where(sub_cloud eq icloud,num_cloud) ind_clear=where(sub_cloud eq 0,num_clear) xy_clear=fltarr(2,num_clear) ;the location of each clear pixel xy_clear[1,*]=fix(ind_clear/(a_region)) xy_clear[0,*]=ind_clear-(a_region)*xy_clear[1,*] fine2_clear=fltarr(num_clear,nb) fine1_clear=fltarr(num_clear,nb) for ib=0,nb-1,1 do begin fine2_clear[*,ib]=(sub_fine2[*,*,ib])[ind_clear] fine1_clear[*,ib]=(sub_fine1[*,*,ib])[ind_clear] endfor for ic=0.0,num_cloud-1.0,1.0 do begin ri=float(fix(ind_cloud[ic]/(a_region))) ; the location of target pixel ci=float(ind_cloud[ic]-(a_region)*ri) center_dis=((x_center-ci)^2+(y_center-ri)^2)^0.5; the distance between target pixel and center of cloud s_dis=((xy_clear[0,*]-ci)^2+(xy_clear[1,*]-ri)^2)^0.5 sub_on_common=fltarr(num_clear,nb) for iband=0,nb-1,1 do begin sub_on_common[*,iband]= fine2_clear[*,iband]-sub_fine2[ci,ri,iband] endfor rmsei=mean(sub_on_common^2,DIMENSION=2)^0.5+0.0001 order_clear=sort(rmsei,/L64) ;order all common pixels by rmse with the target pixel iclear=0.0 isimilar=0.0 fine1_similar=fltarr(min_pixel,nb) fine2_similar=fltarr(min_pixel,nb) rmse_similar=fltarr(min_pixel) dis_similar=fltarr(min_pixel) similar_rmse12=fltarr(min_pixel) maxi_num_seach=min([num_clear,20.0*min_pixel]) ;limit the searching pixels to increase the efficiency while (isimilar le min_pixel-1 and iclear le maxi_num_seach-1.0 ) do begin indicate_similar=0 for iband=0,nb-1,1 do begin if (abs(fine2_clear[order_clear[iclear],iband]-sub_fine2[ci,ri,iband])le similar_th_band[iband]) then begin indicate_similar=indicate_similar+1 endif endfor if (indicate_similar eq nb) then begin for ib=0, nb-1,1 do begin fine1_similar[isimilar,ib]=(fine1_clear[*,ib])[order_clear[iclear]] fine2_similar[isimilar,ib]=(fine2_clear[*,ib])[order_clear[iclear]] endfor rmse_similar[isimilar]=((total((fine2_clear[order_clear[iclear],*]-sub_fine2[ci,ri,*])^2))/float(nb))^0.5 dis_similar[isimilar]=s_dis[order_clear[iclear]] similar_rmse12[isimilar]=((total((fine2_clear[order_clear[iclear],*]-fine1_clear[order_clear[iclear],*])^2))/float(nb))^0.5 isimilar=isimilar+1 endif iclear=iclear+1.0 endwhile ; if similar pixel larger than 0 if (isimilar gt 0 ) then begin ind_null=where(dis_similar gt 0) rmse_similar1=(rmse_similar[ind_null]-min(rmse_similar[ind_null]))/(max(rmse_similar[ind_null])-min(rmse_similar[ind_null])+0.000001)+1.0 dis_similar1=(dis_similar[ind_null]-min(dis_similar[ind_null]))/(max(dis_similar[ind_null])-min(dis_similar[ind_null])+0.000001)+1.0 C_D=(rmse_similar1)*dis_similar1+0.0000001 weight=(1.0/C_D)/total(1.0/C_D) ;compute the time weight W_T1=center_dis/(center_dis+mean(dis_similar[ind_null])) W_T2=mean(dis_similar[ind_null])/(center_dis+mean(dis_similar[ind_null])) for iband=0,nb-1,1 do begin predict_1=total((fine1_similar[*,iband])[ind_null]*weight) predict_2=sub_fine2[ci,ri,iband]+total(((fine1_similar[*,iband])[ind_null]-(fine2_similar[*,iband])[ind_null])*weight) if (predict_2 gt DN_min and predict_2 lt DN_max) then begin sub_fine1[ci,ri,iband]=W_T1*predict_1+W_T2*predict_2 endif else begin sub_fine1[ci,ri,iband]=predict_1 endelse endfor endif else begin ; if no similar pixel, use all similar pixels for iband=0,nb-1,1 do begin diff=mean(fine1_clear[*,iband]-fine2_clear[*,iband]) sub_fine1[ci,ri,iband]=sub_fine2[ci,ri,iband]+diff endfor endelse endfor fine1[left_region:right_region,up_region:down_region,*]=sub_fine1 endif endfor endif print,'finished ',isub+1,' patch' ; Envi_Write_Envi_File,fine1,Out_Name = tempoutname+strtrim(isub+1,1) envi_file_mng, id=Fid1, /remove, /delete envi_file_mng, id=Fid2, /remove, /delete envi_file_mng, id=Fid3, /remove, /delete endfor ;-------------------------------------------------------------------------------------- ;mosiac all the filled patch mfid=intarr(n_ns*n_nl) mdims=intarr(5,n_ns*n_nl) mpos=intarr(nb,n_ns*n_nl) pos=indgen(nb) x0=intarr(n_ns*n_nl) y0=intarr(n_ns*n_nl) for isub=0,n_ns*n_nl-1,1 do begin envi_open_file, tempoutname+strtrim(isub+1,1), r_fid= sub_fid if (sub_fid eq -1) then begin envi_batch_exit return endif envi_file_query, sub_fid, ns=sub_ns, nl=sub_nl mfid[isub] = sub_fid mpos[*,isub] = indgen(nb) mdims[*,isub] = [-1,0, sub_ns-1,0, sub_nl-1] x0[isub]=ind_patch[0,isub] y0[isub]=ind_patch[2,isub] endfor xsize = orig_ns ysize = orig_nl pixel_size = [1.,1.] use_see_through = replicate(1L,n_ns*n_nl) see_through_val = replicate(0L,n_ns*n_nl) ; out_name=Dialog_PickFile(Title = 'Enter the filename of the restored image') out_name=FileName1+'cloud_remove' envi_doit, 'mosaic_doit', fid=mfid, pos=mpos, $ dims=mdims, out_name=out_name, xsize=xsize, $ ysize=ysize, x0=x0, y0=y0, georef=0,MAP_INFO=map_info, $ out_dt=Data_Type, pixel_size=pixel_size, $ background=0, see_through_val=see_through_val, $ use_see_through=use_see_through for i=0,n_ns*n_nl-1,1 do begin envi_file_mng, id=mfid[i], /remove, /delete endfor print, 'time used', floor((systime(1)-t0)/3600), 'hour',floor(((systime(1)-t0) mod 3600)/60),'m',(systime(1)-t0) mod 60,'s' end这是完整代码,代码运行的环境是ENVI5.6和IDL 8.8.0 (Win32 x86_64 m64).,根据这些信息解决上面的问题
10-08
% 文件路径设置 hdr_file = 'C:\Users\11218\Desktop\aagg.hdr'; spe_file = 'C:\Users\11218\Desktop\aagg.spe'; % 指定波段索引 red_band = 117; % 红色波段 green_band = 71; % 绿色波段 blue_band = 26; % 蓝色波段 try % ============================================= % 步骤1: 手动解析HDR文件 % ============================================= fid = fopen(hdr_file, 'r'); hdr_text = textscan(fid, '%s', 'Delimiter', '\n'); fclose(fid); hdr_text = hdr_text{1}; % 初始化元数据结构体 hdr_info = struct(); % 解析关键参数 for i = 1:length(hdr_text) line = strtrim(hdr_text{i}); % 跳过空行和注释 if isempty(line) || line(1) == '#' || line(1) == ';' continue; end % 分割键值对 [key, value] = strtok(line, '='); if isempty(value) [key, value] = strtok(line); else value = strtrim(value(2:end)); % 去除等号 end key = lower(strtrim(key)); % 处理特殊格式的值 if contains(value, '{') && contains(value, '}') value = value(2:end-1); % 去除花括号 end % 存储元数据 switch key case 'samples' hdr_info.samples = str2double(value); case 'lines' hdr_info.lines = str2double(value); case 'bands' hdr_info.bands = str2double(value); case 'data type' hdr_info.data_type = str2double(value); % 转换为MATLAB数据类型 switch hdr_info.data_type case 1 % 8-bit byte matlab_type = 'uint8'; case 2 % 16-bit signed integer matlab_type = 'int16'; case 3 % 32-bit signed long integer matlab_type = 'int32'; case 4 % 32-bit floating point matlab_type = 'single'; case 5 % 64-bit double matlab_type = 'double'; case 12 % 16-bit unsigned integer matlab_type = 'uint16'; case 13 % 32-bit unsigned long integer matlab_type = 'uint32'; otherwise error('未知的数据类型: %d', hdr_info.data_type); end hdr_info.matlab_type = matlab_type; case 'interleave' hdr_info.interleave = lower(value); case 'byte order' hdr_info.byte_order = str2double(value); end end % ============================================= % 步骤2: 读取SPE文件数据 % ============================================= % 设置字节顺序 if isfield(hdr_info, 'byte_order') if hdr_info.byte_order == 0 byte_order = 'ieee-le'; % 小端序 else byte_order = 'ieee-be'; % 大端序 end else byte_order = 'n'; % 使用系统默认字节顺序 end % 打开并读取二进制数据 fid = fopen(spe_file, 'r', byte_order); total_elements = hdr_info.samples * hdr_info.lines * hdr_info.bands; data = fread(fid, total_elements, ['*' hdr_info.matlab_type]); fclose(fid); % ============================================= % 步骤3: 重塑数据为三维数组 % ============================================= data = reshape(data, hdr_info.samples, hdr_info.lines, hdr_info.bands); % 根据interleave类型调整维度顺序 switch hdr_info.interleave case 'bsq' % Band Sequential (波段顺序) % 无需调整,MATLAB默认是[行,列,波段] data = permute(data, [2, 1, 3]); case 'bil' % Band Interleaved by Line (按行交叉) data = reshape(data, hdr_info.samples, hdr_info.bands, hdr_info.lines); data = permute(data, [3, 1, 2]); case 'bip' % Band Interleaved by Pixel (按像素交叉) data = reshape(data, hdr_info.bands, hdr_info.samples, hdr_info.lines); data = permute(data, [3, 2, 1]); otherwise error('未知的interleave类型: %s', hdr_info.interleave); end % ============================================= % 步骤4: 提取RGB波段并合成图像 % ============================================= % 提取波段数据 red = data(:, :, red_band); green = data(:, :, green_band); blue = data(:, :, blue_band); % 创建RGB合成图像 rgb = cat(3, red, green, blue); % 转换为double类型以便处理 rgb_double = double(rgb); % ============================================= % 步骤5: 对比度增强 % ============================================= % 线性拉伸 (2%-98%) p_low = 2; p_high = 98; % 计算百分数 low_vals = prctile(rgb_double(:), p_low); high_vals = prctile(rgb_double(:), p_high); % 应用线性拉伸 rgb_stretched = (rgb_double - low_vals) / (high_vals - low_vals); rgb_stretched(rgb_stretched < 0) = 0; rgb_stretched(rgb_stretched > 1) = 1; % ============================================= % 步骤6: 显示并保存图像 % ============================================= figure; imshow(rgb_stretched); title(sprintf('RGB Composite (R:%d, G:%d, B:%d)', red_band, green_band, blue_band)); catch ME fprintf('处理过程中出错: %s\n', ME.message); % 显示详细错误信息 fprintf('错误发生在: %s (行号 %d)\n', ME.stack(1).name, ME.stack(1).line); rethrow(ME); end 删去对错误的相关处理,简化代码
最新发布
11-26
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值