OpenCVで平均濃度の求め方

質問

OpenCVを使って画像の平均濃度を求めたいのですが、どうやるのでしょうか。
0～256レベルで、例えば0～85の間での平均濃度の求め方です。
濃度ヒストグラムと同じなのでしょうか？

回答No.6

machongola · ベストアンサー

こんばんは。
　調査を再開した所、以下の様なものを見つけました。

　http://www.mis.med.akita-u.ac.jp/~kata/image/equalize.html
　http://www.icnet.ne.jp/~nsystem/simd_tobira/hist_stretch.html

　暗い濃度に偏ったヒストグラムを補正する手法です。
　此れを利用して補正をしてから4値化を行った所、くっきりとした変換が出来ました。
　画像は一番上のURLに在ったものを御拝借しています。参考になれば。

//ヒストグラム配列最大個数
const int HISTOGRAM_MAX = 256;

//データ構造
struct DATA
{
//ヒストグラム配列
int histogram[HISTOGRAM_MAX];

//最小濃度
int min;

//最大濃度
int max;
};

//データの初期化
static void InitData(DATA* p)
{
const DATA zero = {{0}, 255, 0};
*p = zero;
}

//first - (last - 1)までの累積
static int Ruiseki(const int histogram[], int first, int last)
{
int sum = 0;
for(int i = first; i < last; ++i)
sum += histogram[i];
return sum;
}

//first - (last - 1)までのヒストグラムで最も大きい濃度位置を返す
static int Greater(const int histogram[], int first, int last)
{
int max = 0;
int idx = 0;
for(int i = first; i < last; ++i)
if(max < histogram[i])
{
max = histogram[i];
idx = i;
}
return idx;
}

//first - 255までの濃度位置に対する累積と平均累計を比較して閾値を決める
static int DecideThreashold(const int histogram[], int first, int heikinRuiseki)
{
for(int i = first; i < HISTOGRAM_MAX; ++i)
{
//first - (first + i)の累積を取る
const int ruiseki = ::Ruiseki(histogram, first, i);

//平均累積と一致
if(ruiseki == heikinRuiseki)
{
return i;
}
//平均累積を超えた
else if(ruiseki > heikinRuiseki)
{
//濃度位置が0超過の時
if(i > 0)
{
//濃度位置-1に対する累積と現在の濃度位置に対する累積のうち、平均累積から誤差の低い方を取る
if(::abs(heikinRuiseki - ruiseki) > ::abs(heikinRuiseki - ::Ruiseki(histogram, first, i - 1)))
return i - 1;
}
return i;
}
}

//255超過　first - 255までのヒストグラムで最も大きい濃度位置を返す
return ::Greater(histogram, first, HISTOGRAM_MAX);
}

//ヒストグラムの形成
static void CalcHistogram(DATA* p, const IplImage* ipl)
{
for(int y = 0; y < ipl->height; ++y)
{
for(int x = 0; x < ipl->width; ++x)
{
const int pos = (x * ipl->nChannels) + (ipl->widthStep * y);
const int value = (unsigned char)ipl->imageData[pos];

//濃度に対応したヒストグラムをカウントする
p->histogram[value]++;

//最小濃度であれば記憶
if(p->min > value)
p->min = value;

//最大濃度であれば記憶
if(p->max < value)
p->max = value;
}
}
}

//ヒストグラムの補正
static void Equalizing(const DATA* p, IplImage* ipl)
{
for(int y = 0; y < ipl->height; ++y)
{
for(int x = 0; x < ipl->width; ++x)
{
const int pos = (x * ipl->nChannels) + (ipl->widthStep * y);
const int value = (unsigned char)ipl->imageData[pos];

//補正前濃度に対応した累積
const int hvalue= ::Ruiseki(p->histogram, 0, value);

//最大濃度に対応した累積
const int hmax = ::Ruiseki(p->histogram, 0, p->max);

//最小濃度に対応した累積
const int hmin = ::Ruiseki(p->histogram, 0, p->min);

//計算しておく
const int range = hmax - hmin;

//補正結果
const int result = ((hvalue - hmin) * 255) / (range == 0 ? 1 : range);

//255を超過させないで代入
ipl->imageData[pos] = result > 255 ? 255 : result;
}
}
}

//4値化
static void Yonchika(const DATA* p, IplImage* ipl)
{
//総累積を3で割る
const int heikinRuiseki = ::Ruiseki(p->histogram, 0, HISTOGRAM_MAX) / 3;

//t1, t2, t3の閾値を得る
const int t1 = ::DecideThreashold(p->histogram, 0, heikinRuiseki);
const int t2 = ::DecideThreashold(p->histogram, t1, heikinRuiseki);
const int t3 = ::DecideThreashold(p->histogram, t2, heikinRuiseki);

for(int y = 0; y < ipl->height; ++y)
{
for(int x = 0; x < ipl->width; ++x)
{
const int pos = (x * ipl->nChannels) + (ipl->widthStep * y);
const int value = (unsigned char)ipl->imageData[pos];

//閾値と比較して4値化する
if(value >= 0 && value < t1)
ipl->imageData[pos] = 85*0;
else if(value >= t1 && value < t2)
ipl->imageData[pos] = 85*1;
else if(value >= t2 && value < t3)
ipl->imageData[pos] = 85*2;
else if(value >= t3 && value < HISTOGRAM_MAX)
ipl->imageData[pos] = 85*3;
else
;//エラー
}
}
}

//お試し
int main()
{
DATA data;
//データの初期化
::InitData(&data);
//IPLをグレースケールで読み込む
IplImage* ipl = ::cvLoadImage("test.bmp", 0);
//IPLからヒストグラムを形成する
::CalcHistogram(&data, ipl);
//IPLを補正する
::Equalizing(&data, ipl);

//データの初期化
::InitData(&data);
//補正後のIPLからヒストグラムを形成する
::CalcHistogram(&data, ipl);
//4値化する
::Yonchika(&data, ipl);
//セーブする
::cvSaveImage("yonchika.bmp", ipl);
//IPLの解放
::cvReleaseImage(&ipl);
return 0;
}

お礼

ありがとうございます。
確かにこれだとしっかり変換できます。
こんな方法があるとは考えつきませんでした。
とても参考になりました！

回答No.5

machongola

こんばんは。補足頂きました。

　constに関して値を変更しない事を明示的にするつもりで行っているだけですので、在っても無くても構わないです。
　当方はC++をメインで使用しているので、値を変更しない場合にconstを付けるクセがあります。

　unsingedの方ですが、IplImage構造体のimageDataメンバがchar*型なので、参照すると符号付のchar型です。(0-255)の間で処理を進める為に、符号を外しました。

回答No.4

machongola

こんにちは。補足頂きました。

　そう言う事でしたか。
　(1)(2)に付いてですが、

　一応、以下も参考にはしていたのですが、
　http://mikilab.doshisha.ac.jp/dia/research/person/shuto/research/0605/two-level.html

　例えば、

(1)濃度ヒストグラム[0-85]迄を加算して行く。

　は

　histogram[0] = 5
　histogram[1] = 11
　・
　・
　・
　histogram[255] = 24

　とあった場合、

　count = 5 + 11 + 省略 + 24

　と言う事で、

(2)（濃度ヒストグラム[0-85]　×　濃度[0-85]）を加算して行く。

　は

　histogram[0] = 5
　histogram[1] = 11
　・
　・
　・
　histogram[255] = 24

　とあった場合、

　histogram[255]で言えば、濃度255が 24個ある訳ですから、255 x 24 と言う事です。

　color = (0 x 5) + (1 x 11) + (省略) + (255 x 24)

　そして最後に

　avg = color / count

　で平均を算出したつもりでした。
　若しかしたら、此れでは無理なのかもしれません・・・。

補足

大変遅れて申し訳ありません。
そういった平均なのですね。
学校で平均濃度をならったはずなのですが、ノートがないのでこういう計算方法であってたかな？と思った所存です。

ところでプログラムのところどころに
const int
const unsigned char Y
(unsigned char)img->imageData[pos];
などの修飾子が見受けられますが、このプログラムではどういった意味をもっているのでしょうか。
constは指定した変数が定数であることを指定し、unsignedは負の数がないことだったと思いますが、修飾子をつける意味はあるのでしょうか？
なくても大丈夫なのでしょうか？

回答No.3

machongola

こんばんは。補足頂きました。

・return (count > 0) ? (color / count) : 0;　について

　if(count > 0)
　　return color / count;
　else
　　return 0;

　と同義です。こうしている理由は、countの数字が「0」であった時、

　color / 0

　と言う事に成ります。
　「0」で割り算をすると、プログラムが停止してしまうので、其れを回避する為に、この様にしています。
　OpenCVの4値化とは関係の無い話ですが、「0」で割り算をしても停止しない様にする方法もある様です。

・4値化について

　当方の書き方が悪かった部分があるのですが、Yの出力は [0-255] 迄で、0 を含めて 256段階と言うのが正確な解釈です。

　そのソースで言えば、

　t_3 <= Y <= 256

　ですので Y=256 の時も t_3 の閾値を使用する事になります。
　と言う事は、Y=256（真っ白）である筈が、t_3 の数字に置き換わる為、真っ白に成らない（明るいグレー）事が多いと言う事です。
　しかし、逆に Y=256 の時に、別の数字を入れてしまうと4値化の枠を超えてしまいます（4色よりも多い色彩を持つ事になる）。
　如何しても色彩の断層を一律にするのであるならば、

　if(0<=Y && Y<=t_1){
　　img->imageData[pos] = 85 * 0;
　}else if(t_1<=Y && Y<=t_2){
　　img->imageData[pos] = 85 * 1;
　}else if(t_2<=Y && Y<=t_3){
　　img->imageData[pos] = 85 * 2;
　}else if(t_3<=Y && Y<=256){
　　img->imageData[pos] = 85 * 3;
　}

　とすれば出来ます（この手法で問題なければ）。

・余談ですが、以前立てられていた質問の

　http://oshiete1.goo.ne.jp/qa4885645.html

　(1)～(3)は、こう言う事でしょうか。

//4値化
static void Yonchika(const IplImage* img, const int histogram[], int t1, int t2, int t3)
{
for(int y = 0; y < img->height; ++y)
{
for(int x = 0; x < img->width; ++x)
{
const int pos = (x * img->nChannels) + (img->widthStep * y);
const unsigned char Y = (unsigned char)img->imageData[pos];

const int a1 = ::Avg(histogram, 0, t1);
const int a2 = ::Avg(histogram, t1, t2);
const int a3 = ::Avg(histogram, t2, t3);
const int a4 = ::Avg(histogram, t3, 256);

if(Y >= 0 && Y < t1)
img->imageData[pos] = a1;
else if(Y >= t1 && Y < t2)
img->imageData[pos] = a2;
else if(Y >= t2 && Y < t3)
img->imageData[pos] = a3;
else if(Y >= t3 && Y < 256)
img->imageData[pos] = a4;

t1 = (a1 + a2) / 2;
t2 = (a2 + a3) / 2;
t3 = (a3 + a4) / 2;
}
}
}

//テスト
void test()
{
//グレースケール変換して読み込み
int histogram[256] = {0};
IplImage* img = ::cvLoadImage("test.bmp", 0);

//ヒストグラムの作成
::CreateHistogram(img, histogram);

//[0-84][85-119][120-239][240-255]でスタート
::Yonchika(img, histogram, 85, 120, 240);

//ファイルへセーブとIPL解放
::cvSaveImage("yonchika.bmp", img);
::cvReleaseImage(&img);
}

補足

ありがとうございます。
以前質問させていただいた時に教えてもらったプログラムでも4値化はできるのですが、１つ問題がありました。
それは、ヒストグラムに偏りがある(例えば白い部分がなくほぼ黒に近い画像)と4値化されないということです。
なので、今回は画像の平均濃度を求めて、その閾値により確実に4値化させることがしたかったのです。

ところで平均濃度なのですが
(1)濃度ヒストグラム[0-85]迄を加算して行く。
(2)（濃度ヒストグラム[0-85]　×　濃度[0-85]）を加算して行く。
後は(2)を(1)で割ればよい。
ということなのですが、よろしければこの原理について教えていただいてもよろしいでしょうか。

回答No.2

machongola

こんばんは。補足頂きました。

・const unsigned char Y = (unsigned char)img->imageData[pos];
・++histogram[Y];
　イメージをグレースケールに変換して読み込んでいるので、イメージデータの画素は全て[0-255]に変換されています。
　この数字が色彩の暗い・明るいに該当する為、例えば　Y=23　の時

　++histogram[23]　であれば

　濃度23の発見個数をカウントアップすると言う事です。此れを繰り返していくとhistogram[0-255]の中にグラフが出切る筈です。

・histgram[]とcountの型が一致しないのが原因かなと思ったりもしたのですが、どう思われますか？
　はい。unsigned char histogram[256]では完全なサイズ不足でした。
　此れでは各濃度の発見個数が255個迄に成ってしまうので、int histogram[256]に修正しました。

　ディスプレイして確認して見ましたが、当方の平均濃度算出に対する解釈が短絡的過ぎるのかもしれません。
　此れで不可能であるならば、ギブアップです。

　一応画像の説明をすると、
　波形は、山の絵（色彩あり）のヒストグラムです。
　白い線は　s = (t1 + t2) / 2　で算出し続けて見つけた閾値です。
　山の絵（モノクロ）は、見つけた閾値を境にして2値化したものです。

//ヒストグラム[first-last]迄の平均濃度を求める
static int Avg(const int histogram[], int first, int last)
{
int count = 0;
int color = 0;
for(int i = first; i < last; ++i)
{
count += histogram[i];
color += i * histogram[i];
}
return (count > 0) ? (color / count) : 0;
}

//[0-s]と[s-256]領域を計算して閾値を得る
static int FindThreshold(const int histogram[], int sCur)
{
//最有力候補の閾値(一番動きが少なかった物）
int sPri = 0;

//histogram[0-256]迄
for(int i = 0; i < 256; ++i)
{
//[0-s]の平均濃度
const int t1 = ::Avg(histogram, 0, sCur);

//[s-256]の平均濃度
const int t2 = ::Avg(histogram, sCur, 256);

//新たな閾値を作成する
const int sNew = (t1 + t2) / 2;

//新たな閾値と一歩前の閾値を引いて差を出す
const int d1 = ::abs(sNew - sCur);

//新たな閾値と最有力候補の閾値を引いて差を出す
const int d2 = ::abs(sNew - sPri);

//新たな閾値の動きが全く無いので此処で確定
if(d1 == 0)
return sNew;

//新たな閾値の方が最有力候補の閾値よりも動きが小さいので、新たな閾値を最有力候補にする
if(d1 < d2)
sPri = sNew;

//t1, t2を求める為に代入する
sCur = sNew;
}
return sPri;
}

//histogram[0-256]の作成
static void CreateHistogram(const IplImage* img, int histogram[])
{
for(int y = 0; y < img->height; ++y)
{
for(int x = 0; x < img->width; ++x)
{
//x, yに対応した画素の位置計算
const int pos = (x * img->nChannels) + (img->widthStep * y);

//画素を取る
const unsigned char Y = (unsigned char)img->imageData[pos];

//ヒストグラムを加算する
++histogram[Y];
}
}
}

//確認
void test()
{
//グレースケール変換で読み込む
int histogram[256] = {0};
IplImage* img = ::cvLoadImage("test.bmp", 0);

//ヒストグラム作成
::CreateHistogram(img, histogram);

//閾値を探させる
const int threshold = ::FindThreshold(histogram, 85);

//閾値を元に2値化をして見る
::cvThreshold(img, img, threshold, 255, CV_THRESH_BINARY);

//ファイルへ書き込む
::cvSaveImage("nichika.bmp", img);

//IPLの解放
::cvReleaseImage(&img);
}

補足

ありがとうございます。
こちらの方でもできました。

質問ですが
return (count > 0) ? (color / count) : 0;
これはどういう意味なのでしょうか。

あと、自分なりにアレンジして、4値化できるプログラムを作ったのですが、4値化ですと閾値は３つ必要ですよね。
新閾値がt1,t2,t3と出てきたとして、[0～t1][t1～t2][t2～t3][t3～256]の領域ができます。
img->imageData[pos]を参照してその値がどの領域の中にあるのかを調べて、それに適当な値をimageDataに戻すとき、その戻す値をどうすればいいでしょうか。

for(int y = 0; y < img->height; ++y){
for(int x = 0; x < img->width; ++x){
const int pos = (y * img->widthStep) + (x * img->nChannels);
const unsigned char Y = (unsigned char)img->imageData[pos];
if(0<=Y && Y<=t_1){

img->imageData[pos] = 0;
}else if(t_1<=Y && Y<=t_2){

img->imageData[pos] = t_1;
}else if(t_2<=Y && Y<=t_3){

img->imageData[pos] = t_2;
}else if(t_3<=Y && Y<=256){

img->imageData[pos] = t_3;
}
}
}
これが4値値に振り分けるときの部分ですが、この場合imageData[pos]に戻す値をimageData[pos]の値が[0～t1]なら0,[t1～t2]ならt1,[t2～t3]ならt2,[t3～256]ならt3としてますが、じゃあ256はどうするの？ということです。
わかりにくくてすみませんが、いかがでしょうか。
よろしくお願いします。

回答No.1

machongola

こんにちは。

(1)濃度ヒストグラム[0-85]迄を加算して行く。
(2)（濃度ヒストグラム[0-85]　×　濃度[0-85]）を加算して行く。

　後は(2)を(1)で割ればよいのでは。

//histogram[first-last]迄の平均濃度を求める
static int Avg(const unsigned char histogram[], int first, int last)
{
//histogram[first-last]の合計数
int count = 0;

//histogram[first-last]の合計濃度
unsigned long color = 0;

//first-last迄回転する
for(int i = first; i <= last; ++i)
{
//数を加算する
count += histogram[i];

//濃度を加算する
color += i * histogram[i];
}

//平均を返す
return (count > 0) ? (color / count) : 0;
}

//確認
void test()
{
//ビットマップをグレースケール変換でロード
IplImage* img = ::cvLoadImage("test.bmp", 0);

//濃度ヒストグラム
unsigned char histogram[256] = {0};

//濃度ヒストグラムの作成
for(int y = 0; y < img->height; ++y)
{
for(int x = 0; x < img->width; ++x)
{
const int pos = (x * img->nChannels) + (img->widthStep * y);
const unsigned char Y = (unsigned char)img->imageData[pos];
++histogram[Y];
}
}

//濃度ヒストグラム[0-85]迄の平均濃度を計算する
const int result = ::Avg(histogram, 0, 85);

//IPLの解放
::cvReleaseImage(&img);
}

　では駄目でしょうか。

補足

ありがとうございます。
やってみたのですが、うまくいきませんでした。

個人的に
const unsigned char Y = (unsigned char)img->imageData[pos];
++histogram[Y];
この辺がよくわからなかったのと、histgram[]とcountの型が一致しないのが原因かなと思ったりもしたのですが、どう思われますか？

あと、質問の補足を致しますと、256レベルで閾値sを決めますと、[0～s]と[s～256]の２つの領域ができます。
それぞれの領域で平均濃度を求め、その値をt1,t2とします。
s=(t1+t2)/2を計算し、新たな閾値sとして、sの値がほとんど変わらなくなるまで繰り返してその閾値を求めたいのですが平均濃度で足踏み状態というわけなのですが、いかがでしょうか。
よろしくお願いします。