Main Content

大きなイメージの統計値の計算

ライブ スクリプトを開く

この例では、blockproc を使用して大きなイメージの統計値を計算し、その情報を使用してイメージをブロック単位で処理する精度を高める方法を示します。関数 blockproc は、イメージ ブロックに演算を適用したり、結果を集計したり、結果を新しいイメージとして返す際に最適です。多くのイメージ処理アルゴリズムは、イメージに関する "グローバルな" 情報を必要としますが、1 つのイメージ ブロックのみを検討する場合はこのような情報を利用できません。こうした制約は、大きすぎてメモリに完全には読み込めないイメージを処理する場合に問題となる可能性があります。

この例では、複数のスペクトルのカラー合成イメージの強調の例と類似する作業を行いますが、blockproc を使用して大きなイメージ用に調整されています。この例では、Erdas LAN ファイル rio.lan の可視帯域を強調します。これらのタイプのブロック処理手法は、通常、大きなイメージで効果的ですが、この例では小さなイメージを使用して概念を説明します。

手順 1: トゥルーカラー合成の作成

blockproc を使用して、 rio.lan からデータ (Landsat Thematic Mapper 撮影イメージを Erdas LAN ファイル形式で含んでいるファイル) を読み取ります。blockproc には TIFF ファイルと JPEG2000 ファイルのみを読み取る組み込みサポートがあります。他のタイプのファイルを読み取るには、イメージ アダプター クラスを記述して特定のファイル形式の I/O をサポートしなければなりません。この例では、あらかじめ組み込まれているイメージ アダプター クラス、LanAdapter を使用します。この関数は LAN ファイルの読み取りをサポートします。イメージ アダプター クラスの記述方法の詳細については、サポートされない形式のイメージ ファイルに対するブロック処理の実行を参照してください。

Erdas LAN 形式には帯域 3、2、および 1 にそれぞれ可視赤色、緑色、および青色スペクトルが含まれています。blockproc を使用して可視帯域を RGB イメージに抽出します。

rio.lan に関連付けられた LanAdapter オブジェクトを作成します。

input_adapter = LanAdapter("rio.lan");

可視の R 帯、G 帯、B 帯を選択します。

input_adapter.SelectedBands = [3 2 1];

ブロック データを変更せずに返すシンプルなブロック関数を作成します。

identityFcn = @(block_struct) block_struct.data;

初期のトゥルーカラー イメージを作成します。

truecolor = blockproc(input_adapter,[100 100],identityFcn);

強調前の結果を表示します。

imshow(truecolor)
title("Truecolor Composite (Not Enhanced)")

結果のトゥルーカラー イメージは、複数のスペクトルのカラー合成イメージの強調の例にある paris.lan のイメージに似ています。RGB イメージはぼんやりしており、コントラストがほとんどありません。

手順 2: イメージの強調 - 最初の試行

最初に、blockproc を使用してデータをダイナミック レンジ全体に引き伸ばします。この最初の試みでは、データの各ブロックに stretchlim と imadjust を個別に呼び出す新しい関数ハンドルのみを定義します。

adjustFcn = @(block_struct) imadjust(block_struct.data, ...
    stretchlim(block_struct.data));
truecolor_enhanced = blockproc(input_adapter,[100 100],adjustFcn);

imshow(truecolor_enhanced)
title("Truecolor Composite with Blockwise Contrast Stretch")

結果が間違っていることがすぐにわかります。問題は、関数 stretchlim が入力イメージ上のヒストグラムを計算し、この情報を使用してストレッチ範囲を計算することです。各ブロックは隣接するブロックから分離して調整されるため、ブロックはそれぞれのローカル ヒストグラムから異なる範囲が計算されます。

手順 3: ヒストグラムのアキュムレータ クラスの確認

イメージのダイナミック レンジ全般にわたるデータの分布を確認するには、3 つの可視帯域のヒストグラムをそれぞれ計算します。

十分に大きなイメージを処理する場合は、imhist を呼び出すだけではイメージ ヒストグラムを作成できません。徐々にヒストグラムを作成する 1 つの方法は、イメージを移動するときに各ブロックのヒストグラムを合計するクラスを blockproc で使用することです。

HistogramAccumulator クラスを調べます。

type HistogramAccumulator
% HistogramAccumulator Accumulate incremental histogram.
%   HistogramAccumulator is a class that incrementally builds up a
%   histogram for an image.  This class is appropriate for use with 8-bit
%   or 16-bit integer images and is for educational purposes ONLY.

%   Copyright 2009 The MathWorks, Inc.

classdef HistogramAccumulator < handle
   
    properties
        Histogram
        Range
    end
    
    methods
        
        function obj = HistogramAccumulator()
            obj.Range = [];
            obj.Histogram = [];
        end
        
        function addToHistogram(obj,new_data)
            if isempty(obj.Histogram)
                obj.Range = double(0:intmax(class(new_data)));
                obj.Histogram = hist(double(new_data(:)),obj.Range);
            else
                new_hist = hist(double(new_data(:)),obj.Range);
                obj.Histogram = obj.Histogram + new_hist;
            end
        end
    end
end

このクラスは関数 hist を呼び出す単純なラッパーで、データをヒストグラムに徐々に追加できるようにします。このラッパーは blockproc 固有ではありません。HistogramAccumulator クラスの次の簡単な使用例を確認します。

HistogramAccumulator オブジェクトを作成します。

hist_obj = HistogramAccumulator;

サンプル イメージを 2 等分します。

full_image = imread("liftingbody.png");
top_half = full_image(1:256,:);
bottom_half = full_image(257:end,:);

ヒストグラムをインクリメンタルに計算します。

addToHistogram(hist_obj,top_half);
addToHistogram(hist_obj,bottom_half);
computed_histogram = hist_obj.Histogram;

関数 imhist の結果と比較します。

normal_histogram = imhist(full_image);

結果を確認します。ヒストグラムは数値的にまったく同じです。

figure
subplot(1,2,1)
stem(computed_histogram,"Marker","none")
title("Incrementally Computed Histogram")
subplot(1,2,2)
stem(normal_histogram',"Marker","none")
title("imhist Histogram")

手順 4: blockproc での HistogramAccumulator クラスの使用

HistogramAccumulator クラスを blockproc で使用して、rio.lan のデータの赤色帯域のヒストグラムを作成します。データの各ブロックに addToHistogram メソッドを呼び出す blockproc の関数ハンドルを定義できます。ヒストグラムを表示すると、利用可能なダイナミック レンジの一部にデータが集中しているのがわかります。他の可視帯域にも同様の分布があります。これが、元のトゥルーカラー合成がぼんやりと見える一因です。

HistogramAccumulator オブジェクトを作成します。

hist_obj = HistogramAccumulator;

blockproc 関数ハンドルを設定します。

addToHistFcn = @(block_struct) addToHistogram(hist_obj, block_struct.data);

赤チャネルのヒストグラムを計算します。addToHistFcn 関数ハンドルは出力を何も生成しないことに注意してください。blockproc に渡された関数ハンドルは何も返さないため、blockproc も何も返しません。

input_adapter.SelectedBands = 3;
blockproc(input_adapter,[100 100],addToHistFcn);
red_hist = hist_obj.Histogram;

結果を表示します。

figure
stem(red_hist,"Marker","none")
title("Histogram of Red Band (Band 3)")

手順 5: コントラスト ストレッチを使用したトゥルーカラー合成の強調

イメージに適切なコントラスト ストレッチを行うことができます。通常のインメモリ ワークフローの場合、関数 stretchlim を使用するだけで imadjust の引数を計算できます (複数のスペクトルのカラー合成イメージの強調の例と同様です)。しかし、これまで見てきたように、大きなイメージを処理する場合、stretchlim はフル イメージ ヒストグラムに依存するため簡単には blockproc に合わせて使用できません。

各可視帯域のイメージ ヒストグラムを計算したら、imadjust の正しい引数を手作業で計算します (stretchlim の場合と同様です)。

まず、緑色帯域と青色帯域のヒストグラムを計算します。

hist_obj = HistogramAccumulator;
addToHistFcn = @(block_struct) addToHistogram(hist_obj,block_struct.data);
input_adapter.SelectedBands = 2;
blockproc(input_adapter,[100 100],addToHistFcn);
green_hist = hist_obj.Histogram;

hist_obj = HistogramAccumulator;
addToHistFcn = @(block_struct) addToHistogram(hist_obj,block_struct.data);
input_adapter.SelectedBands = 1;
blockproc(input_adapter,[100 100],addToHistFcn);
blue_hist = hist_obj.Histogram;

各ヒストグラムの CDF を計算します。

computeCDF = @(histogram) cumsum(histogram) / sum(histogram);
findLowerLimit = @(cdf) find(cdf > 0.01, 1, "first");
findUpperLimit = @(cdf) find(cdf >= 0.99, 1, "first");

red_cdf = computeCDF(red_hist);
red_limits(1) = findLowerLimit(red_cdf);
red_limits(2) = findUpperLimit(red_cdf);

green_cdf = computeCDF(green_hist);
green_limits(1) = findLowerLimit(green_cdf);
green_limits(2) = findUpperLimit(green_cdf);

blue_cdf = computeCDF(blue_hist);
blue_limits(1) = findLowerLimit(blue_cdf);
blue_limits(2) = findUpperLimit(blue_cdf);

関数 imadjust の引数を準備します。

rgb_limits = [red_limits' green_limits' blue_limits'];

[0, 1] の範囲にスケールします。

rgb_limits = (rgb_limits - 1) / (255);

新しい adjustFcn を作成してグローバルなストレッチ範囲を適用し、blockproc を使用してトゥルーカラー イメージを調整します。

adjustFcn = @(block_struct) imadjust(block_struct.data,rgb_limits);

完全な RGB データを選択します。

input_adapter.SelectedBands = [3 2 1];
truecolor_enhanced = blockproc(input_adapter,[100 100],adjustFcn);

結果を表示します。結果のイメージが大幅に改善され、データはダイナミック レンジに大きく広がりました。blockproc を使用すると、メモリにイメージ全体が読み込まれることはありません。

imshow(truecolor_enhanced)
title("Truecolor Composite with Corrected Contrast Stretch")

参考

クラス

関数

関連する例

詳細