K1号館のエレベータ内映像における有彩色の分析

画像内の彩度が0.1以上のピクセルを双六角錐色モデル上にプロットした結果を示す。 また、色相のヒストグラム(0〜359[degree]を72分割)を算出して、最大頻度の領域を示した。

図1(a):エレベータA(青) 図1(b):双六角錐色モデル上での色配置

図2(a):エレベータB(赤) 図2(b):双六角錐色モデル上での色配置

図3(a):エレベータC(オレンジ) 図3(b):双六角錐色モデル上での色配置

図4(a):エレベータD(緑) 図4(b):双六角錐色モデル上での色配置

表1:エレベータ内映像において最大頻度を有する色相領域

対象名配列添え字範囲[degree]
エレベータA(青) 47 235〜240
エレベータB(赤) 67 335〜340
エレベータC(オレンジ) 4 20〜25
エレベータD(緑) 37 185〜190

ColorModel2.hpp
#pragma once
	
#include <cmath>
	
struct rgb {
	unsigned char r,g,b;
	
	rgb() {}
	rgb(unsigned char r, unsigned char g, unsigned char b) : r(r), g(g), b(b) {}
};
	
double max2(double a, double b)
{
	if (a >= b) return a;
	
	return b;
}
	
double max3(double a, double b, double c)
{
	return max2(max2(a,b),c);
}
	
double min2(double a, double b)
{
	if (a <= b) return a;
	
	return b;
}
	
double min3(double a, double b, double c)
{
	return min2(min2(a,b),c);
}
	
inline unsigned char clamp(int x){
	if (x > 255) return (unsigned char)255;
	else if (x < 0) return (unsigned char)0;
	
	return (unsigned char)x;
}
	
double calcX(double h, double I, double S)
{
	double _pi = atan(1.0)*4;
	double _2pi = _pi*2;
	
	double dh = h;
	if (dh < 0) {
		dh = h + _2pi;
	}
	else if (dh > _2pi) {
		dh = h - _2pi;
	}
	
	double M2;
	if (I <= 0.5) {
		M2 = I*(1+S);
	}
	else {
		M2 = I+S-I*S; 
	}
	
	double M1 = 2*I-M2;
	
	double X;
	if (dh < _pi/3) {
		X = (M1 + (M2 - M1)*dh)/(_pi/3);
	}
	else if (dh < _pi) {
		X = M2;
	}
	else if (dh < _pi*4/3) {
		X = (M1 + (M2 - M1)*(_pi*4/3 - dh))/(_pi/3);
	}
	else {
		X = M1;
	}
	
	return X;
}
	
// 双六角錐(十二面体)モデルによるHSI
struct  hsi {
	double H, S, I;
	
	hsi(double H, double S, double I) : H(H), S(S), I(I) {}
	
	hsi(const rgb &c) {
		double B = c.b/255.0;
		double G = c.g/255.0;
		double R = c.r/255.0;
	
		double maxI = max3(R, G, B);
		double minI = min3(R, G, B);
		I = (maxI + minI)/2;
	
		if (I <= 0.5) {
			S = (maxI - minI)*(maxI + minI);
		}
		else {
			S = (maxI - minI)*(2 - (maxI + minI));
		}
	
		if (maxI != minI) {
			double rr = (maxI - R)/(maxI - minI);
			double gg = (maxI - G)/(maxI - minI);
			double bb = (maxI - B)/(maxI - minI);
	
			if (R == maxI) {
				H = atan(1.0) * 4 / 3 * (bb - gg);
			}
			else if (G == maxI) {
				H =  atan(1.0) * 4 / 3 * (2 + rr - bb);
			}
			else {
				H =  atan(1.0) * 4 / 3 * (4 + gg - rr);
			}
	
			if (H < 0) {
				H = H + atan(1.0)*8;
			}
		}
		else {
			H = 0;
		}
	}
	
	rgb to_rgb() {
		if (S == 0) {
			unsigned char _i = clamp((int)(I*255));
			return rgb(_i, _i, _i);
		}
	
		double _pi = atan(1.0)*4;
	
		int _r = (int)(calcX(H+_pi*2/3, I, S) * 255);
		int _g = (int)(calcX(H, I, S) * 255);
		int _b = (int)(calcX(H-_pi*2/3, I, S) * 255);
	
		return rgb(clamp(_r), clamp(_g), clamp(_b));
	}
};

main.cpp
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

#include <iostream>

#include "ColorModel2.hpp"

#ifdef _DEBUG
//Debugモードの場合
#pragma comment(lib,"opencv_core231d.lib")
#pragma comment(lib,"opencv_imgproc231d.lib")
#pragma comment(lib,"opencv_highgui231d.lib")
#else
//Releaseモードの場合
#pragma comment(lib,"opencv_core231.lib")
#pragma comment(lib,"opencv_imgproc231.lib")
#pragma comment(lib,"opencv_highgui231.lib")
#endif

using namespace std;

/*
char *TARGET_FILENAME = "a_blue.jpg";
char *SAVE_FILENAME = "a_blue_hsi.jpg";
char *TITLE = "HSIの図(a_blue.jpg)";
*/

/*
char *TARGET_FILENAME = "b_red.jpg";
char *SAVE_FILENAME = "b_red_hsi.jpg";
char *TITLE = "HSIの図(b_red.jpg)";
*/

/*
char *TARGET_FILENAME = "c_orange.jpg";
char *SAVE_FILENAME = "c_orange_hsi.jpg";
char *TITLE = "HSIの図(c_orange.jpg)";
*/

char *TARGET_FILENAME = "d_green.jpg";
char *SAVE_FILENAME = "d_green_hsi.jpg";
char *TITLE = "HSIの図(d_green.jpg)";

int main()
{
	cv::Mat disp0 = cv::Mat(401, 401, CV_8UC3);
	disp0 = cv::Scalar(255, 255, 255);

	int width = disp0.size().width;
	int height = disp0.size().height;

	int cx = 200, cy = 200;
	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {

			double r = sqrt(1.0*(x-cx)*(x-cx) + 1.0*(y-cy)*(y-cy));
			double _s = r/200.0;

			if (_s <= 1.0) {

				double _h = atan2(1.0*(cy - y), 1.0*(x - cx));
				
				if (_h < 0) _h += (atan(1.0)*8.0);

				hsi hsi_c0(_h, _s, 0.5);
				rgb rgb_c0 = hsi_c0.to_rgb();

				disp0.at<cv::Vec3b>(y,x)[2] = rgb_c0.r;
				disp0.at<cv::Vec3b>(y,x)[1] = rgb_c0.g;
				disp0.at<cv::Vec3b>(y,x)[0] = rgb_c0.b;

			}
		}
	}

	cv::Mat image = cv::imread(TARGET_FILENAME);
	int width2 = image.size().width;
	int height2 = image.size().height;

	int bin = 72;
	int div = 360/72;
	int *hist = new int[bin];
	for (int i = 0; i < bin; i++) {
		hist[i] = 0;
	}

	for (int y = 0; y < height2; y++) {
        for (int x = 0; x < width2; x++) {
			unsigned char _r = image.at<cv::Vec3b>(y,x)[2];
			unsigned char _g = image.at<cv::Vec3b>(y,x)[1];
			unsigned char _b = image.at<cv::Vec3b>(y,x)[0];

			
			rgb rgb_c0(_r, _g, _b);
			hsi hsi_c0(rgb_c0);

			if (hsi_c0.S >= 0.1) {
				int px = (int)(hsi_c0.S * 200 * cos(hsi_c0.H) + 200);
				int py = (int)(200 - hsi_c0.S * 200 * sin(hsi_c0.H));

				if (px >= 0 && px < 401 && py >= 0 && py < 401) {
					disp0.at<cv::Vec3b>(py,px)[2] = 0;
					disp0.at<cv::Vec3b>(py,px)[1] = 0;
					disp0.at<cv::Vec3b>(py,px)[0] = 0;
				}

				int H_angle = (int)(hsi_c0.H * 45.0 / atan(1.0));
				if (H_angle < 0) {
					H_angle += 360;
				}
				hist[H_angle / div]++;
			}
		}
	}

	int max_bin = 0;
	int max_f = hist[0];
	for (int i = 0; i < bin; i++) {
		cout << "[" << i << "] " << hist[i] << endl;
		if (max_f < hist[i]) {
			max_f = hist[i];
			max_bin = i;
		}
	}

	cout << "max_bin = " << max_bin << endl;

	cv::Point ptCenter(cx, cy);
	double angle1 = max_bin * div * atan(1.0) / 45.0;
	double angle2 = (max_bin+1) * div * atan(1.0) / 45.0;
	cv::Point pt1((int)(200.0*cos(angle1)+cx), (int)(cy-200.0*sin(angle1)));
	cv::Point pt2((int)(200.0*cos(angle2)+cx), (int)(cy-200.0*sin(angle2)));
	cv::line(disp0, ptCenter, pt1, cv::Scalar(0, 0, 0));
	cv::line(disp0, ptCenter, pt2, cv::Scalar(0, 0, 0));

	delete [] hist;

	cv::imwrite(SAVE_FILENAME, disp0);
	cv::imshow(TITLE, disp0);
	cv::imshow(TARGET_FILENAME, image);
	cv::waitKey(0);
	return 0;
}