关于遥感影像BIL、BIP、BSQ你知道多少？给一个二进制文件你会读取嘛~

由于微信公众号改变了推送规则，为了每次新的推送可以在第一时间出现在您的订阅列表中，记得将本公众号设为星标或置顶喔~

  介绍了遥感影像通用数据格式BIP、BIL、BSQ，并给出MATLAB、C++从原始二进制数据读取影像的代码~

🌿前言

  遥感影像的通用存储格式有BIP、BIL、BSQ三种，通常在编程处理图像时都会使用GDAL、Arcpy、Rasterio等库来进行读取，如果因某些原因使用不了这些库，又该如何解决呢？

  本文先简单介绍了BIP、BIL、BSQ三种通用格式的区别，随后分析了如何正确读取二进制存储遥感数据，最后给出了MATLAB和C++读取参考代码。

  文末有获取项目所有文件的方式~

🍀BIP&BIL&BSQ

BIP

  BIP(Band Interleaved by Pixel)，可以理解为按照像素存储，先依次存储每个波段的第一个像素，再存储每个波段的第二个像素…

BIL

  BIL(Band Interleaved by Line)，可以理解为按照行进行存储，先依次存储每个波段的第一行，再存储每个波段的第二行…

BSQ

  BSQ(Band Sequential)，即按照顺序存储，依次存第一个波段所有像素，再存第二个波段所有像素…

一图理解

  下图应该能够很好的反映出存储的方式(按照对应的数字顺一遍就很好理解)

了解更多可以参考此链接：https://blog.csdn.net/m0_46387817/article/details/108467947

🌸数据读取思路

重要两点

  有关如何正确地读取影像数据，笔者认为有两点是至关重要的：

1. 弄清所读楚遥感影像的存储是BIP还是BIL还是BSQ
2. 知道所读影像的数据类型：float、double

  如果数据格式不正确，就会导致数据读取的有问题，例如一个存储类型是8位的数据按照16位的数据类型读取，就会直接导致读取到的数组大小发生改变，这个错误还容易发现。

  如果数据格式存储的8位无符号整型，恰好又按照8位的整型去读取，读出来的数组大小一致，但是值已经是错误的值了。

读取步骤

1. 将数据从二进制文件中按照存储类型读到一维数组

  需要根据影像的元数据文件或者一些属性信息，了解影像单个像元存储的数据类型，再利用一些读取二进制文件的函数例如MATLAB中的fopen和fread函数、C++中的ifstream函数将数据正确地读到一维数组。

2. 根据BIP、BIL、BSQ对一维数组进行解析

  根据索引关系将一维数组解析为对应的波段数组。

索引关系

  关于如何建立一维数组和波段数组（三维）的关系，即建立一个$ind与ri、ci、ki$的函数：
$$ind=\text{i}\times ri+\text{j}\times ci+\text{k}\times bi$$
  式中,$ind$表示一维数组中的索引，$ri,ci,bi$分别代表第$ri$行、第$ci$列、第$bi$个波段,$i,j,k$是待求量。

  笔者也是思考了很久，简单总结出如下两种思考方法(其实两种很类似,也挺绕的😴)：

方法一

  即固定任意两个变量的值为0，变化第三个变量的值，观察其规律，解求得到其对应的待求量。

方法二

  先假设图像为一维(即单波段、一行值)，得到一个规律；再逐渐增加维度，观察变换确立关系。

索引关系

🍗MATLAB代码

  笔者根据自己理解，编写了parserData的MATLAB函数。并利用ENVI软件制作了BIL、BIP、BSQ格式的RGB图片文件，使用MATLAB读取数据文件经过解析后简单绘制可视化。

读取效果

BIP读取

BIL读取

BSQ读取

parserData.m

function result = parseData(lineData, row, column, bandcount, type)
    % 解析一维数组到BIP、BIL、BSQ格式
    totalsize = length(lineData);
    result = zeros([row, column, bandcount]);

    % 解析BIP格式
    if upper(type) == "BIP"
        for i = 1:bandcount
            tmp = lineData(i:bandcount:totalsize);
            result(:, :, i) = reshape(tmp, column, row)';
        end
    % 解析BIL格式
    elseif upper(type) == "BIL"
        tmp = reshape(lineData, bandcount*column, row);
        for i = 1:bandcount
            result(:, :, i) = tmp(column*(i - 1)+1:column*i, :)';
        end
    % 解析BSQ格式
    elseif upper(type) == "BSQ"
        tmp = reshape(lineData, column, row*bandcount)';
        for i = 1:bandcount
            result(:, :, i) = tmp(row*(i - 1)+1:row*i, :);
        end
    end
end

demo.m

clear, clc, close all;
type = "BIP";
filepath = strcat("../", type, "/tree", type);
filehead = strcat(filepath, ".hdr");

% 1.读取数据到一维数组
fileID = fopen(filepath, 'r');
if (fileID == -1)
    disp(strcat(filepath, "---OpenFailed"));
end
fileData = fread(fileID, "uint8");

types = ["BIP", "BIL", "BSQ"];

% 2.解析数据，并绘图
for i = 1:length(types)
    figure
    data = parseData(fileData, 465, 300, 3, types(i));
    imshow(uint8(data));
    resultName = strcat("../", type, "/result", types(i), ".png");
    title(strcat(type, "->", types(i)), 'FontSize', 14, 'FontName', 'Times New Roman')
    print(gcf, '-dpng','-r300', resultName);
end

🍟C++代码

  C++主要是利用ifstream函数读取二进制文件，存储到vector数组，编写一个模板类Parser以便解析一维数组。

  由于笔者能力有限，并未做C++的图行可视化，通过编写test函数生成简单一维数组调用类进行解析以测试类的可靠性，通过输出读取图像第一个像元值来与ENVI软件、取色值软件所探测得到的RGB做对比以测试读取的正确性。

Parser.hpp

#pragma once
#include <vector>
using namespace std;

template<typename T>
class Parser {
public:
    vector<T> lineData;
public:
    Parser(const vector<T> lineData);
    void BIP(int r, int c, int b,
        vector<vector<vector<T>>>& result);
    void BIL(int r, int c, int b,
        vector<vector<vector<T>>>& result);
    void BSQ(int r, int c, int b,
        vector<vector<vector<T>>>& result);
};

template<typename T>
Parser<T>::Parser(const vector<T> lineData) {
    this->lineData = lineData;
}

template<typename T>
void Parser<T>::BIP(int r, int c, int b,
    vector<vector<vector<T>>>& result) {
    result.clear();
    for (int bi = 0; bi < b; bi++) {
        vector<vector<T>> bandData(r, vector<T>(c));
        for (int ri = 0; ri < r; ri++) {
            for (int ci = 0; ci < c; ci++) {
                int ind = ci * b + ri * b * c + bi;
                bandData[ri][ci] = lineData[ind];
            }
        }
        result.push_back(bandData);
    }
}
template<typename T>
void Parser<T>::BIL(int r, int c, int b,
    vector<vector<vector<T>>>& result) {
    result.clear();
    for (int bi = 0; bi < b; bi++) {
        vector<vector<T>> bandData(r, vector<T>(c));
        for (int ri = 0; ri < r; ri++) {
            for (int ci = 0; ci < c; ci++) {
                int ind = ci + ri * b * c + bi * c;
                bandData[ri][ci] = lineData[ind];
            }
        }
        result.push_back(bandData);
    }
}
template<typename T>
void Parser<T>::BSQ(int r, int c, int b,
    vector<vector<vector<T>>>& result) {
    result.clear();
    for (int bi = 0; bi < b; bi++) {
        vector<vector<T>> bandData(r, vector<T>(c));
        for (int ri = 0; ri < r; ri++) {
            for (int ci = 0; ci < c; ci++) {
                int ind = bi * r * c + ri * c + ci;
                bandData[ri][ci] = lineData[ind];
            }
        }
        result.push_back(bandData);
    }
}

Demo.cpp

#include <stdio.h>
#include <fstream>
#include <vector>
#include "Parser.hpp"

using namespace std;

template <typename T>
void print(vector < vector<vector<T>>>& result) {
    for (auto& band : result) {
        for (auto& row : band) {
            for (auto& pixel : row) {
                printf("%d ", pixel);
            }
            printf("\n");
        }
        printf("\n");
    }
}

void test() {
    vector<int> lineData(24);
    for (int i = 0; i < 24; i++) {
        lineData[i] = i + 1;
    }
    int row = 2, column = 3, bandcount = 2;
    Parser<int> parser(lineData);
    vector<vector<vector<int>>> result;
    parser.BIP(row, column, bandcount, result);
    //parser.BIL(row, column, bandcount, result);
    //parser.BSQ(row, column, bandcount, result);
    print(result);
}

int main() {
    test(); //简单测试Parser类的函数
    string filepath = "./../../BIP/treeBIP"; //二进制文件路径
    int row = 465, column = 300, bandcount = 3; //自己定义的行、列、波段数

    int size = row * column * bandcount;
    typedef uint8_t datatype;//影像存储数据类型

    //1.读取二进制文件到一维数组
    ifstream inFile(filepath, ios::in | ios::binary);
    if (!inFile) {
        printf("Error: Cannot open file\n");
        return 0;
    }
    vector<datatype> lineData(size);
    inFile.read(reinterpret_cast<char*>(&lineData[0]), size * sizeof(datatype));
    inFile.close();

    //2.解析一维数组
    vector<vector<vector<datatype>>> result;
    Parser<datatype> parser(lineData);
    parser.BIP(row, column, bandcount, result);

    printf("%d %d %d",result[0][0][0],result[1][0][0],result[2][0][0]);
    return 0;
}

🍕项目文件获取

有关本文所涉及的案例代码，均可通过关注本公众号，后台私信回复“241107”获取~

🌹结语

1. 笔者为测绘遥感方向的学习者，愿意结交志同道合的伙伴，欢迎关注公众号与笔者一起交流学习~

2. 关于影像的大小(size)和存储格式等信息，均可从影像数据的元数据文件获取~

3. 在MATLAB代码中，读取后的数组索引应该为[row][column][band]，在C++代码中，索引为[band][row][column]，这是由于MATLAB和C++数组取值便捷性而设计的~

• 路虽远，行则将至；事虽难，做则必成。希望认真学习的您能够有所收获~

• 本公众号的原创成果，在未经允许的情况下，请勿用于任何商业用途！

• 如果本文有幸可以帮到您，欢迎您的点赞👍、收藏⭐与关注❤；您的点赞👍、收藏⭐与关注❤是我创作的最大动力~