ちょいめも

物理/Python/Cの雑記帳

pandas、matplotlibでヒストグラム表示

#csv読み込んでヒストグラム表示
#第一引数:csvフルパス

import sys
import os
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd

argvs = sys.argv #引数取得
argc = len(argvs) #引数の個数

if argc != 2:
	print('引数1 : csvファイル フルパス\n')
	input('終了します。何かキーをおしてください。。。\n')
	sys.exit()

dirname, filename = os.path.split(argvs[1])
df = pd.read_csv(argvs[1]) #csv読み込み

#plt.style.use('ggplot') 
#font = {'family' : 'meiryo'}
#matplotlib.rc('font', **font)

#----------グラフ1つ描画----------
df.plot( y=['a', 'e'], bins=10, alpha=0.5, figsize=(5,4), kind='hist') #αは透過 0~1


#----------グラフ複数描画----------
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 10), sharex = False) #sharexはx軸を共有するかどうか

df.plot(y=['a', 'e'], bins=10, alpha=1.0, kind='hist', stacked=True, ax=axes[0,0]) #stacked:積み上げ
axes[0,0].set_title('ae')
axes[0,0].set_xlabel("X")
ave = df['a'].mean()
axes[0,0].text(ave,7,"ave:%.1f" % ave, color='blue', bbox=dict(facecolor='none', edgecolor='blue', pad=5.0)) #平均の位置にテキストを表示したり

df.plot(y=['b', 'c'], bins=10, alpha=1.0, kind='hist', stacked=True, ax=axes[0,1])
axes[0,1].set_title('bc')
axes[0,1].set_xlabel("X")
ave = df['b'].mean()
axes[0,1].text(0.1,0.5,"ave:%.1f" % ave, transform = axes[0,1].transAxes,
		color='blue', bbox=dict(facecolor='none', edgecolor='blue', pad=5.0)) #textの位置をグラフ軸座標で描画。最大1

df.plot(y=['a', 'b', 'c', 'd', 'e'], bins=10, alpha=1.0, kind='hist', stacked=True, ax=axes[0,2])
axes[0,2].set_title('abcde')
axes[0,2].set_xlabel("X")

df.plot(y=['d'], bins=10, alpha=1.0, kind='hist', stacked=True, ax=axes[1,0])
axes[1,0].set_title('d')
axes[1,0].set_xlabel("X")

df.plot(y=['e'], bins=10, alpha=1.0, kind='hist', stacked=True, ax=axes[1,1])
axes[1,1].set_title('e')
axes[1,1].set_xlabel("X")

axes[1,2].axis('off') #空白

fig.tight_layout()  # タイトルとラベルが被るのを解消

plt.savefig(dirname + '\histogram.png') #複数グラフをcsvと同じ場所に保存


matplotlib覚書

よく忘れるので覚書

import matplotlib.pyplot as plt

#----------折れ線グラフを1つ描画----------
x = [1,2,3,4,5]
y1 = [2,3,5,7,8]
y2 = [2,4,7,10,13]
plt.plot(x, y1, label="1")
plt.plot(x, y2, label="2")
plt.legend()  # 凡例をグラフにプロット
plt.title("sample : $ y = x $") #数式を入れる場合はTeX形式
plt.xlabel("X")
plt.ylabel("Y")
plt.xlim(1, 5)  # yを0-5000の範囲に限定
plt.ylim(0, 15)  # yを0-5000の範囲に限定
plt.hlines([4, 10], 2, 4, linestyles="dashed")  # y=4,10 x=2~4に破線を描画

plt.savefig('image01.png') #グラフをエクスポート
plt.show() #ここで、ここまでのグラフの設定、描画を終了

#----------散布図----------

keisankekka = len(y2)
plt.plot(x, y1, ".", label="len:%.1f" % keisankekka) #ラベルを計算結果にしたり
plt.plot(x, y2, "x", label="2")
plt.legend()  # 凡例をグラフにプロット
plt.text(1.5, 10, 'chushaku') #注釈を入れる、グラフのx、y座標を指定する
plt.savefig('image02.png') #グラフをエクスポート

#----------折れ線グラフを複数描画----------
fig, axes = plt.subplots(nrows = 2, ncols = 2, figsize=(5, 5), sharex = False) #sharex:x軸共有しない

#グラフ1
axes[0,0].plot(x, y1, label="0,0")
axes[0,0].legend() # 凡例を表示
axes[0,0].set_title('1')
axes[0,0].set_xlabel('x')
axes[0,0].set_ylabel('y')
axes[0,0].set_xlim(1, 5)
axes[0,0].set_ylim(0, 10)

axes[0,1].plot(x, y2, label="0,1")
axes[0,1].legend() # 凡例を表示
axes[0,1].set_title('2')
axes[0,1].set_xlabel('x')
axes[0,1].set_ylabel('y')
axes[0,1].set_xlim(1, 3)
axes[0,1].set_ylim(0, 15)

axes[1,0].plot(x, y2, label="1,0")
axes[1,0].legend() # 凡例を表示
axes[1,0].set_title('3')
axes[1,0].set_xlabel('x')
axes[1,0].set_ylabel('y')
axes[1,0].set_xlim(1, 3)
axes[1,0].set_ylim(0, 15)

axes[1,1].axis('off') #空白

fig.tight_layout()  # タイトルとラベルが被るのを解消
plt.savefig('image03.png') #グラフをエクスポート

xlsxを開いてcsvに値を保存

import sys
import os
import re
import glob
import csv
import win32com.client

class func():
		
	def __init__(self, load_path=''):
		
		self.dir, self.file = os.path.split(load_path)
		self.load_path = load_path

	def read_xlsx(self):

		xlApp = win32com.client.Dispatch("Excel.Application")
		xlApp.DisplayAlerts = False
		
		wb = xlApp.Workbooks.Open(self.load_path)
		sheet = wb.Worksheets("Sheet1")
		
		a = sheet.Range("A1").Value
		b = sheet.Range("B1").Value

		wb.Close()
		xlApp.Quit()
		wb = None
		xlApp = None

		if not os.path.isfile('summary.csv'):
			f = open('summary.csv', 'w')
			f.close()
			
		f = open('summary.csv', 'a', newline='')
		csvWriter = csv.writer(f)
		csvWriter.writerow([a,b])
		f.close()
		
if __name__ == '__main__':
	
	argvs = sys.argv 
	argc = len(argvs) #引数の個数
	
	if argc != 2:
		sys.exit()
	
	if re.search('.xlsx',argvs[1]):
		try:
			sheet = func(argvs[1])
			sheet.read_xlsx()
		except:
			print('Error : Input file is not supported.')
			sys.exit()
	else:
		try:
			xlsx_list = glob.glob('%s/*.xlsx' % argvs[1])
			
			if len(xlsx_list) != 0:
				for i in xlsx_list:
					sheet = func(i)
					sheet.read_xlsx()
			else:
				print('Error : Input file/folder is not supported.')
				sys.exit()
			
		except:
			print('Error : Input file is not supported.')
			sys.exit()

csvファイルをExcelで開いて保存する

#カンマ区切りが正しくないcsvファイルを一度エクセルで開いて
#保存しなおすことで修正する

import win32com.client

xlApp = win32com.client.Dispatch("Excel.Application")
xlApp.DisplayAlerts = False

#ファイルを開く
wb = xlApp.Workbooks.Open("D:\\python\\test_err.csv")

#名前をつけて保存
wb.SaveAs("D:\\python\\test_mod.csv")

#ファイルを閉じる
wb.Close()

#エクセル終了
xlApp.Quit()

#以下をしないとプロセスにExcelが残る
wb = None
xlApp = None

画像をシフトして出力と、複数画像を平均して出力(python)

画像を縦横任意のpx数シフトして出力

img_shift.py

# -*- coding: utf-8 -*-
#画像をシフトして保存する
#第一引数:画像ファイルパス
#第二引数:縦シフトpx (rows)
#第三引数:横シフトpx (cols)

import cv2
import numpy as np
import sys
import os

argvs = sys.argv #引数取得
argc = len(argvs) #引数の個数

if argc != 4:
	print('img_shiftについて\n')
	print('引数1 : 画像ファイルパス\n引数2 : 縦シフトpx (rows)\n引数3 : 横シフトpx (cols)\n')
	input('終了します。何かキーをおしてください。。。\n')
	sys.exit()

img = cv2.imread(argvs[1]) #画像読み込み
rows, cols, channels = img.shape

M = np.float32([[1,0,argvs[3]],[0,1,argvs[2]]]) #変換
dst = cv2.warpAffine(img,M,(cols,rows)) #cols,rowsは出力画像サイズ

filename, ext = os.path.splitext(os.path.basename(argvs[1]))
cv2.imwrite(filename + '_shift_row' + argvs[2] + 'col' + argvs[3] + ext, dst) #変換画像保存

出力結果
f:id:choimemo:20180311164701j:plain
↓ img_shift.py test1.jpg 10 20
f:id:choimemo:20180311164705j:plain

フォルダ内に含まれるJPG画像を平均して出力

img_ave.py

# -*- coding: utf-8 -*-
#フォルダ内のjpg画像を平均して別名で保存
#第一引数:画像フォルダパス(jpg画像のみ抽出して平均)

import cv2
import numpy as np
import glob
import sys

argvs = sys.argv #引数取得
argc = len(argvs) #引数の個数

if argc != 2:
	print('img_aveについて\n')
	print('引数1 : 画像ファイルが入ったフォルダのパス\n')
	input('終了します。何かキーをおしてください。。。\n')
	sys.exit()

jpg_list = glob.glob(argvs[1] + '/*.jpg') #フォルダ内ファイル取得

if len(jpg_list) != 0:
	img0 = cv2.imread(jpg_list[0]) #画像読み込み サイズチェック
	rows, cols, channels = img0.shape
	img_ave = np.zeros((rows, cols, channels)) #平均値入れる箱を初期化

	for i in jpg_list:
		img_tmp = cv2.imread(i).astype(int)
		if img_tmp.shape != img0.shape:
			print(i + ' : 画像ファイルサイズが異なります。\n全て同じサイズにしてください。\n')
			input('終了します。何かキーをおしてください。。。\n')
			sys.exit()
		else:
			img_ave = img_ave + img_tmp

	img_ave = img_ave / len(jpg_list)
	cv2.imwrite('img_ave.jpg', img_ave) #平均画像保存
													
else:
	print('エラー:JPGファイルのみの対応です。。。')
	input('終了します。何かキーをおしてください。。。\n')
	sys.exit()

出力結果
f:id:choimemo:20180311164701j:plainf:id:choimemo:20180311164753j:plainf:id:choimemo:20180311164808j:plainf:id:choimemo:20180311164810j:plainf:id:choimemo:20180311164812j:plain

f:id:choimemo:20180311164815j:plain

CIELABとかの覚書

#include <stdio.h>
#include <math.h>
#define _USE_MATH_DEFINES
#define DEG (180 / M_PI)
#define RAD (M_PI / 180)

//CIEXYZからCIELabに変換 Ymax = 100で入力
void XYZ2Lab(double X, double Y, double Z, double* L, double* a, double* b){
	
	double fx, fy, fz;
	double X_Xn, Y_Yn, Z_Zn;
	double Xn = 95.046;	//D65白色点
	double Yn = 100.;
	double Zn = 108.906;
	
	X_Xn = X / Xn;
	Y_Yn = Y / Yn;
	Z_Zn = Z / Zn;
	
	if (X_Xn > pow((6. / 29.), 3.)){
		fx = pow(X_Xn, 1. / 3.);
	}else{
		fx = (pow((29. / 3.), 3.) * X_Xn + 16.) / 116.;
	}

	if (Y_Yn > pow((6. / 29.), 3.)){
		fy = pow(Y_Yn, 1. / 3.);
	}else{
		fy = (pow((29. / 3.), 3.) * Y_Yn + 16.) / 116.;
	}

	if (Z_Zn > pow((6. / 29.), 3.)){
		fz = pow(Z_Zn, 1. / 3.);
	}else{
		fz = (pow((29. / 3.), 3.) * Z_Zn + 16.) / 116.;
	}

	*L = 116. * fy - 16.;
	*a = 500. * (fx - fy);
	*b = 200. * (fy - fz);
	
}

//CIEXYZからsRGBに変換 Ymax = 100で入力 8bitRGB
void XYZ2sRGB(double X, double Y, double Z, double* R, double* G, double* B){

	int i, j;
	double R_tmp, G_tmp, B_tmp;
	double M_inv[3][3] = 	{{3.240970, -1.537383, -0.498611},
							{-0.969244, 1.875968, 0.041555},
							{0.055630, -0.203977, 1.056972}};	// XYZ→RGBの変換行列 sRGB(D65)
	
	//ガンマ補正前のRGB
	R_tmp = M_inv[0][0] * X / 100. + M_inv[0][1] * Y / 100. + M_inv[0][2] * Z / 100.;
	G_tmp = M_inv[1][0] * X / 100. + M_inv[1][1] * Y / 100. + M_inv[1][2] * Z / 100.;
	B_tmp = M_inv[2][0] * X / 100. + M_inv[2][1] * Y / 100. + M_inv[2][2] * Z / 100.;
	
	//ガンマ補正
	if(R_tmp <= 0.0031308){
		*R = (12.92 * R_tmp) * 255.; 
	}else{
		*R = (1.055 * pow(R_tmp, (1. / 2.4)) - 0.055) * 255.;
	}

	if(G_tmp <= 0.0031308){
		*G = (12.92 * G_tmp) * 255.;
	}else{
		*G = (1.055 * pow(G_tmp, (1. / 2.4)) - 0.055) * 255.;
	}
	
	if(B_tmp <= 0.0031308){
		*B = (12.92 * B_tmp) * 255.;
	}else{
		*B = (1.055 * pow(B_tmp, (1. / 2.4)) - 0.05) * 255;
	}
	
}

//RGBからHTMLカラーコードに変換する
void RGB2HTMLcolorcode(int R, int G, int B, char* code, int* dec){

	sprintf(code, "0x%02x%02x%02x", R, G, B);
	*dec = (R << 16) + (G << 8) + B;

}

//ΔEabを計算する
double deltaEab(double L1, double a1, double b1, double L2, double a2, double b2){
	
	return sqrt(pow(L2 - L1, 2) + pow(a2 - a1, 2) + pow(b2 - b1, 2));
	
}

//ΔE00を計算する
double deltaE00(double L1, double a1, double b1, double L2, double a2, double b2){
	
	double kL, kC, kH;
	double dLp, Lbar;
	double Cs_1, Cs_2, Cbar;
	double ap_1, ap_2;
	double Cp_1, Cp_2, Cpbar, dCp;
	double hp_1, hp_2;
	double dsmallhp, dHp, Hpbar;
	double T;
	double S_L, S_C, S_H;
	double R_C, dtheta, R_T;
	double term1, term2, term3, term4;
	
	kL = 1.;	//kL
	kC = 1.;	//kC
	kH = 1.;	//kH
	dLp = L2 - L1;		//ΔL'
	Lbar = (L1 + L2) / 2.;		//Lbar
	Cs_1 = sqrt(pow(a1, 2.) + pow(b1, 2.));		//C*_1
	Cs_2 = sqrt(pow(a2, 2.) + pow(b2, 2.));		//C*_2
	Cbar = (Cs_1 + Cs_2) / 2.;		//Cbar
	ap_1 = a1 + (a1 / 2.) * (1. - sqrt(pow(Cbar, 7.) / (pow(Cbar, 7.) + pow(25., 7.))));	//a'_1
	ap_2 = a2 + (a2 / 2.) * (1. - sqrt(pow(Cbar, 7.) / (pow(Cbar, 7.) + pow(25., 7.))));	//a'_2
	Cp_1 = sqrt(pow(ap_1, 2.) + pow(b1, 2.));		//C'_1
	Cp_2 = sqrt(pow(ap_2, 2.) + pow(b2, 2.));			//C'_2
	Cpbar = (Cp_1 + Cp_2) / 2.;		//C'bar
	dCp = Cp_2 - Cp_1;		//ΔC'
	hp_1 = DEG * atan2(b1, ap_1);		//h'_1
	hp_2 = DEG * atan2(b2, ap_2);		//h'_2
	//Δh'
	if (Cp_1 == 0. || Cp_2 == 0.){
		dsmallhp = 0.;
	}else if(fabs(hp_1 - hp_2) <= 180.){
		dsmallhp = hp_2 - hp_1;
	}else if((fabs(hp_1 - hp_2) > 180.) && (hp_2 <= hp_1)){
		dsmallhp = hp_2 - hp_1 + 360.;
	}else{
		dsmallhp = hp_2 - hp_1 - 360.;
	}
	dHp = 2. * sqrt(Cp_1 * Cp_2) * sin(RAD * (dsmallhp/ 2.));		//ΔH'
	//H'bar
	if (Cp_1 == 0. || Cp_2 == 0.){
		Hpbar = hp_1 + hp_2;
	}else if(fabs(hp_1 - hp_2) <= 180.){
		Hpbar = (hp_1 + hp_2) / 2.;
	}else if((fabs(hp_1 - hp_2) > 180.) && (hp_1 + hp_2 < 360)){
		Hpbar = (hp_1 + hp_2 + 360.) / 2.;
	}else{
		Hpbar = (hp_1 + hp_2 - 360.) / 2.;
	}
	T = 1. - 0.17 * cos(RAD * (Hpbar - 30.)) + 0.24 * cos(RAD * (2. * Hpbar)) + 0.32 * cos(RAD * (3. * Hpbar + 6.)) - 0.20 * cos(RAD * (4. * Hpbar - 63.));		//T
	S_L = 1. + (0.015 * pow((Lbar - 50.), 2.)) / (sqrt(20. + pow((Lbar - 50.), 2.)));		//S_L
	S_C = 1. + 0.045 * Cpbar;		//S_C
	S_H = 1. + 0.015 * Cpbar * T;		//S_H
	R_C = 2. * sqrt(pow(Cpbar, 7.) / (pow(Cpbar, 7.) + pow(25., 7.)));		//R_C
	dtheta = 30. * exp(-1. * pow(((Hpbar - 275.) / 25.), 2.));		//Δθ
	R_T = -1. * sin(RAD * (2 * dtheta)) * R_C;		//R_T
	term1 = pow(dLp / (kL * S_L), 2.);		//第1項
	term2 = pow(dCp / (kC * S_C), 2.);		//第2項
	term3 = pow(dHp / (kH * S_H), 2.);		//第3項
	term4 = R_T * (dCp / (kC * S_C)) * (dHp / (kH * S_H));		//第4項
	
	return sqrt(term1 + term2 + term3 + term4);
}

int main(void){

	double _L, _a, _b;
	double R, G, B;
	char code[8];
	int dec;
	
	XYZ2Lab(30., 30., 30., &_L, &_a, &_b);
	printf("%f\n%f\n%f\n", _L, _a, _b);
	printf("\n");
	
	XYZ2sRGB(30., 30., 30., &R, &G, &B);
	printf("%f\n%f\n%f\n", R, G, B);
	printf("\n");
	
	RGB2HTMLcolorcode(15, 15, 15, code, &dec);
	printf("%s\n", code); //文字列16進数で表示
	printf("%d\n", dec);  //10進数で表示
	printf("\n");
	
	//ΔEab計算
	printf("%f\n", deltaEab(83, -50, -4, 55, -15, 7));
	//ΔE00計算
	printf("%f\n", deltaE00(83, -50, -4, 55, -15, 7));
	
	return 0;

}

EXCELでΔE00の計算シートも作成、Labの色も表示するマクロコード追加
DeltaE.xlsm - Google ドライブ

以下VBAコードも

'D65光源の白色点を返す
Function D65wp()

    Dim result(2) As Double '結果(配列)を返す
    
    result(0) = 95.046   'Xn
    result(1) = 100#     'Yn
    result(2) = 108.906  'Zn
    
    D65wp = result

End Function

'D65 分光分布 380-700nmまでの10nmおきのデータ
Function D65spectrum()

    Dim D65(32) As Double '結果(配列)を返す
    
    D65(0) = 49.9755  '380nm
    D65(1) = 54.6482  '390nm
    D65(2) = 82.7549  '400nm
    D65(3) = 91.486   '410nm
    D65(4) = 93.4318  '420nm
    D65(5) = 86.6823  '430nm
    D65(6) = 104.865  '440nm
    D65(7) = 117.008  '450nm
    D65(8) = 117.812  '460nm
    D65(9) = 114.861  '470nm
    D65(10) = 115.923 '480nm
    D65(11) = 108.811 '490nm
    D65(12) = 109.354 '500nm
    D65(13) = 107.802 '510nm
    D65(14) = 104.79  '520nm
    D65(15) = 107.689 '530nm
    D65(16) = 104.405 '540nm
    D65(17) = 104.046 '550nm
    D65(18) = 100#    '560nm
    D65(19) = 96.3342 '570nm
    D65(20) = 95.788  '580nm
    D65(21) = 88.6856 '590nm
    D65(22) = 90.0062 '600nm
    D65(23) = 89.5991 '610nm
    D65(24) = 87.6987 '620nm
    D65(25) = 83.2886 '630nm
    D65(26) = 83.6992 '640nm
    D65(27) = 80.0268 '650nm
    D65(28) = 80.2146 '660nm
    D65(29) = 82.2778 '670nm
    D65(30) = 78.2842 '680nm
    D65(31) = 69.7213 '690nm
    D65(32) = 71.6091 '700nm

    D65spectrum = D65

End Function

'等色関数 10度視野 380-700nmまでの10nmおきのデータ
Function xyz_func_10deg()

    Dim xyz(32, 3) As Double '結果(配列)を返す
    
    '等色関数x 10度視野
    xyz(0, 0) = 0.00016 '380nm
    xyz(1, 0) = 0.002362 '390nm
    xyz(2, 0) = 0.01911 '400nm
    xyz(3, 0) = 0.084736 '410nm
    xyz(4, 0) = 0.204492 '420nm
    xyz(5, 0) = 0.314679 '430nm
    xyz(6, 0) = 0.383734 '440nm
    xyz(7, 0) = 0.370702 '450nm
    xyz(8, 0) = 0.302273 '460nm
    xyz(9, 0) = 0.195618 '470nm
    xyz(10, 0) = 0.080507 '480nm
    xyz(11, 0) = 0.016172 '490nm
    xyz(12, 0) = 0.003816 '500nm
    xyz(13, 0) = 0.037465 '510nm
    xyz(14, 0) = 0.117749 '520nm
    xyz(15, 0) = 0.236491 '530nm
    xyz(16, 0) = 0.376772 '540nm
    xyz(17, 0) = 0.529826 '550nm
    xyz(18, 0) = 0.705224 '560nm
    xyz(19, 0) = 0.878655 '570nm
    xyz(20, 0) = 1.01416 '580nm
    xyz(21, 0) = 1.11852 '590nm
    xyz(22, 0) = 1.12399 '600nm
    xyz(23, 0) = 1.03048 '610nm
    xyz(24, 0) = 0.856297 '620nm
    xyz(25, 0) = 0.647467 '630nm
    xyz(26, 0) = 0.431567 '640nm
    xyz(27, 0) = 0.268329 '650nm
    xyz(28, 0) = 0.152568 '660nm
    xyz(29, 0) = 0.081261 '670nm
    xyz(30, 0) = 0.040851 '680nm
    xyz(31, 0) = 0.019941 '690nm
    xyz(32, 0) = 0.009577 '700nm

    '等色関数y 10度視野
    xyz(0, 1) = 0.000017 '380nm
    xyz(1, 1) = 0.000253 '390nm
    xyz(2, 1) = 0.002004 '400nm
    xyz(3, 1) = 0.008756 '410nm
    xyz(4, 1) = 0.021391 '420nm
    xyz(5, 1) = 0.038676 '430nm
    xyz(6, 1) = 0.062077 '440nm
    xyz(7, 1) = 0.089456 '450nm
    xyz(8, 1) = 0.128201 '460nm
    xyz(9, 1) = 0.18519 '470nm
    xyz(10, 1) = 0.253589 '480nm
    xyz(11, 1) = 0.339133 '490nm
    xyz(12, 1) = 0.460777 '500nm
    xyz(13, 1) = 0.606741 '510nm
    xyz(14, 1) = 0.761757 '520nm
    xyz(15, 1) = 0.875211 '530nm
    xyz(16, 1) = 0.961988 '540nm
    xyz(17, 1) = 0.991761 '550nm
    xyz(18, 1) = 0.99734 '560nm
    xyz(19, 1) = 0.955552 '570nm
    xyz(20, 1) = 0.868934 '580nm
    xyz(21, 1) = 0.777405 '590nm
    xyz(22, 1) = 0.658341 '600nm
    xyz(23, 1) = 0.527963 '610nm
    xyz(24, 1) = 0.398057 '620nm
    xyz(25, 1) = 0.283493 '630nm
    xyz(26, 1) = 0.179828 '640nm
    xyz(27, 1) = 0.107633 '650nm
    xyz(28, 1) = 0.060281 '660nm
    xyz(29, 1) = 0.0318 '670nm
    xyz(30, 1) = 0.015905 '680nm
    xyz(31, 1) = 0.007749 '690nm
    xyz(32, 1) = 0.003718 '700nm

    '等色関数z 10度視野
    xyz(0, 2) = 0.000705 '380nm
    xyz(1, 2) = 0.010482 '390nm
    xyz(2, 2) = 0.086011 '400nm
    xyz(3, 2) = 0.389366 '410nm
    xyz(4, 2) = 0.972542 '420nm
    xyz(5, 2) = 1.55348 '430nm
    xyz(6, 2) = 1.96728 '440nm
    xyz(7, 2) = 1.9948 '450nm
    xyz(8, 2) = 1.74537 '460nm
    xyz(9, 2) = 1.31756 '470nm
    xyz(10, 2) = 0.772125 '480nm
    xyz(11, 2) = 0.415254 '490nm
    xyz(12, 2) = 0.218502 '500nm
    xyz(13, 2) = 0.112044 '510nm
    xyz(14, 2) = 0.060709 '520nm
    xyz(15, 2) = 0.030451 '530nm
    xyz(16, 2) = 0.013676 '540nm
    xyz(17, 2) = 0.003988 '550nm
    xyz(18, 2) = 0 '560nm
    xyz(19, 2) = 0 '570nm
    xyz(20, 2) = 0 '580nm
    xyz(21, 2) = 0 '590nm
    xyz(22, 2) = 0 '600nm
    xyz(23, 2) = 0 '610nm
    xyz(24, 2) = 0 '620nm
    xyz(25, 2) = 0 '630nm
    xyz(26, 2) = 0 '640nm
    xyz(27, 2) = 0 '650nm
    xyz(28, 2) = 0 '660nm
    xyz(29, 2) = 0 '670nm
    xyz(30, 2) = 0 '680nm
    xyz(31, 2) = 0 '690nm
    xyz(32, 2) = 0 '700nm

    xyz_func_10deg = xyz

End Function


'物体の反射率(%で入力)または透過率からCIEXYZに変換する関数 配列で返す
Function trans2XYZ(trans As range)

    Dim result(2, 0) As Double '結果(配列)を行方向に返す
    
    '定数kを求める
    i = 0
    teisuk = 0
    For Each i In trans
        teisuk = teisuk + D65spectrum(j) * xyz_func_10deg(j, 1)
        j = j + 1
    Next
    
    teisuk = 100# / teisuk
    
    'Xを求める
    i = 0
    j = 0
    X = 0
    For Each i In trans
        X = X + i * 0.01 * D65spectrum(j) * xyz_func_10deg(j, 0)
        j = j + 1
    Next
    X = teisuk * X
    
    'Yを求める
    i = 0
    j = 0
    Y = 0
    For Each i In trans
        Y = Y + i * 0.01 * D65spectrum(j) * xyz_func_10deg(j, 1)
        j = j + 1
    Next
    Y = teisuk * Y
    
    'Zを求める
    i = 0
    j = 0
    Z = 0
    For Each i In trans
        Z = Z + i * 0.01 * D65spectrum(j) * xyz_func_10deg(j, 2)
        j = j + 1
    Next
    Z = teisuk * Z

    result(0, 0) = X
    result(1, 0) = Y
    result(2, 0) = Z
    
    trans2XYZ = result
    
End Function

'CIEXYZからCIELabに変換する関数 配列で返す
Function XYZ2Lab(X, Y, Z)

    Dim result(2, 0) As Double '結果(配列)を行方向に返す

    'D65白色点
    Xn = D65wp(0)
    Yn = D65wp(1)
    Zn = D65wp(2)
    
    X_Xn = X / Xn
    Y_Yn = Y / Yn
    Z_Zn = Z / Zn
    
    If X_Xn > ((6# / 29#) ^ 3#) Then
        fx = X_Xn ^ (1# / 3#)
    Else
        fx = ((29# / 3#) ^ 3# * X_Xn + 16#) / 116#
    End If

    If Y_Yn > ((6# / 29#) ^ 3#) Then
        fy = Y_Yn ^ (1# / 3#)
    Else
        fy = ((29# / 3#) ^ 3# * Y_Yn + 16#) / 116#
    End If

    If Z_Zn > ((6# / 29#) ^ 3#) Then
        fz = Z_Zn ^ (1# / 3#)
    Else
        fz = ((29# / 3#) ^ 3# * Z_Zn + 16#) / 116#
    End If

    L = 116# * fy - 16#
    a = 500# * (fx - fy)
    b = 200# * (fy - fz)
    
    result(0, 0) = L
    result(1, 0) = a
    result(2, 0) = b
    
    XYZ2Lab = result

End Function

'物体の反射率または透過率をLabに変換する関数 配列で返す
Function trans2lab(trans As range)

    Dim result(2, 0) As Double '結果(配列)を行方向に返す
    
    'XYZを求める
    X = trans2XYZ(trans)(0, 0)
    Y = trans2XYZ(trans)(1, 0)
    Z = trans2XYZ(trans)(2, 0)
    
    'Labを求める
    result(0, 0) = XYZ2Lab(X, Y, Z)(0, 0)
    result(1, 0) = XYZ2Lab(X, Y, Z)(1, 0)
    result(2, 0) = XYZ2Lab(X, Y, Z)(2, 0)
    
    trans2lab = result
    
End Function

'ΔEabを返す関数
Function deltaEab(L1, a1, b1, L2, a2, b2)

    deltaEab = Sqr((L2 - L1) ^ 2 + (a2 - a1) ^ 2 + (b2 - b1) ^ 2)

End Function

'ΔE00を返す関数
Function deltaE00(L1, a1, b1, L2, a2, b2)

    kL = 1#     'kL
    kC = 1#     'kC
    kH = 1#     'kH
    dLp = L2 - L1               'ΔL'
    Lbar = (L1 + L2) / 2#       'Lbar
    Cs_1 = Sqr(a1 ^ 2# + b1 ^ 2#)  'C*_1
    Cs_2 = Sqr(a2 ^ 2# + b2 ^ 2#)  'C*_2
    Cbar = (Cs_1 + Cs_2) / 2#      'Cbar
    ap_1 = a1 + (a1 / 2#) * (1# - Sqr(Cbar ^ 7# / (Cbar ^ 7# + 25# ^ 7#))) 'a'_1
    ap_2 = a2 + (a2 / 2#) * (1# - Sqr(Cbar ^ 7# / (Cbar ^ 7# + 25# ^ 7#))) 'a'_2
    Cp_1 = Sqr(ap_1 ^ 2# + b1 ^ 2#)    'C'_1
    Cp_2 = Sqr(ap_2 ^ 2# + b2 ^ 2#)    'C'_2
    Cpbar = (Cp_1 + Cp_2) / 2#         'C'bar
    dCp = Cp_2 - Cp_1                  'ΔC'
    hp_1 = WorksheetFunction.Degrees(WorksheetFunction.Atan2(ap_1, b1))       'h'_1
    hp_2 = WorksheetFunction.Degrees(WorksheetFunction.Atan2(ap_2, b2))       'h'_2
    'Δh'
    If ((Cp_1 = 0#) Or (Cp_2 = 0#)) Then
        dsmallhp = 0#
    ElseIf (Abs(hp_1 - hp_2) <= 180#) Then
        dsmallhp = hp_2 - hp_1
    ElseIf ((Abs(hp_1 - hp_2) > 180#) And (hp_2 <= hp_1)) Then
        dsmallhp = hp_2 - hp_1 + 360#
    Else
        dsmallhp = hp_2 - hp_1 - 360#
    End If
    dHp = 2# * Sqr(Cp_1 * Cp_2) * Sin(WorksheetFunction.Radians(dsmallhp / 2#))   'ΔH'
    'H'bar
    If ((Cp_1 = 0#) Or (Cp_2 = 0#)) Then
        Hpbar = hp_1 + hp_2
    ElseIf (Abs(hp_1 - hp_2) <= 180#) Then
        Hpbar = (hp_1 + hp_2) / 2#
    ElseIf ((Abs(hp_1 - hp_2) > 180#) And (hp_1 + hp_2 < 360)) Then
        Hpbar = (hp_1 + hp_2 + 360#) / 2#
    Else
        Hpbar = (hp_1 + hp_2 - 360#) / 2#
    End If
    T = 1# - 0.17 * Cos(WorksheetFunction.Radians(Hpbar - 30#)) _
           + 0.24 * Cos(WorksheetFunction.Radians(2# * Hpbar)) _
           + 0.32 * Cos(WorksheetFunction.Radians(3# * Hpbar + 6#)) _
           - 0.2 * Cos(WorksheetFunction.Radians(4# * Hpbar - 63#))       'T
    S_L = 1# + (0.015 * (Lbar - 50#) ^ 2#) / (Sqr(20# + (Lbar - 50#) ^ 2#))    'S_L
    S_C = 1# + 0.045 * Cpbar        'S_C
    S_H = 1# + 0.015 * Cpbar * T    'S_H
    R_C = 2# * Sqr(Cpbar ^ 7# / (Cpbar ^ 7# + 25# ^ 7#))      'R_C
    dtheta = 30# * Exp(-1# * ((Hpbar - 275#) / 25#) ^ 2#) 'Δθ
    R_T = -1# * Sin(WorksheetFunction.Radians(2 * dtheta)) * R_C                 'R_T
    term1 = (dLp / (kL * S_L)) ^ 2#     '第1項
    term2 = (dCp / (kC * S_C)) ^ 2#     '第2項
    term3 = (dHp / (kH * S_H)) ^ 2#     '第3項
    term4 = R_T * (dCp / (kC * S_C)) * (dHp / (kH * S_H))       '第4項
    
    deltaE00 = Sqr(term1 + term2 + term3 + term4)
            
End Function

'CIElabからCIEXYZに変換する関数 配列で返す
Function lab2XYZ(L, a, b)

    Dim result(2, 0) As Double '結果(配列)を行方向に返す
    
    '白色点D65
    Xn = D65wp(0)
    Yn = D65wp(1)
    Zn = D65wp(2)
    
    '---CIELabからCIEXYZに変換
    fy = (L + 16#) / 116#
    fx = fy + (a / 500#)
    fz = fy - (b / 200#)

    If fx > (6# / 29#) Then
        X = fx ^ 3# * Xn
    Else
        X = (3# / 29#) ^ 3# * (116# * fx - 16#) * Xn
    End If
    
    If fy > (6# / 29#) Then
        Y = fy ^ 3# * Yn
    Else
        Y = (3# / 29#) ^ 3# * (116# * fy - 16#) * Yn
    End If
    
    If fz > (6# / 29#) Then
        Z = fz ^ 3# * Zn
    Else
        Z = (3# / 29#) ^ 3# * (116# * fz - 16#) * Zn
    End If

    result(0, 0) = X
    result(1, 0) = Y
    result(2, 0) = Z
    
    lab2XYZ = result

End Function

'CIEXYZからsRGBに変換する関数 配列で返す
Function XYZ2RGB(X, Y, Z)

    Dim result(2, 0) As Double '結果(配列)を行方向に返す
    
    'XYZ→RGBの変換行列 sRGB(D65)
    M_inv00 = 3.24097
    M_inv01 = -1.537383
    M_inv02 = -0.498611
    M_inv10 = -0.969244
    M_inv11 = 1.875968
    M_inv12 = 0.041555
    M_inv20 = 0.05563
    M_inv21 = -0.203977
    M_inv22 = 1.056972

    'ガンマ補正前のRGB
    R_tmp = M_inv00 * X / 100# + M_inv01 * Y / 100# + M_inv02 * Z / 100#
    G_tmp = M_inv10 * X / 100# + M_inv11 * Y / 100# + M_inv12 * Z / 100#
    B_tmp = M_inv20 * X / 100# + M_inv21 * Y / 100# + M_inv22 * Z / 100#

    'ガンマ補正
    If R_tmp <= 0.0031308 Then
        r = (12.92 * R_tmp) * 255#
    Else
        r = (1.055 * R_tmp ^ (1# / 2.4) - 0.055) * 255#
    End If

    If G_tmp <= 0.0031308 Then
        g = (12.92 * G_tmp) * 255#
    Else
        g = (1.055 * G_tmp ^ (1# / 2.4) - 0.055) * 255#
    End If
    
    If B_tmp <= 0.0031308 Then
        b = (12.92 * B_tmp) * 255#
    Else
        b = (1.055 * B_tmp ^ (1# / 2.4) - 0.055) * 255#
    End If
       
    result(0, 0) = r
    result(1, 0) = g
    result(2, 0) = b
    
    XYZ2RGB = result

End Function

'Labの色で着色 CIELabからCIEXYZに変換後、sRGBに変換して指定セルを着色
Sub Coloring_Cell()

    i = 9
    Do While Cells(4, i) <> "" '比較対象の1から数値入力されているところまでループ
    
        '---CIELabからCIEXYZに変換
        L = Cells(4, i)
        a = Cells(5, i)
        b = Cells(6, i)
                
        X = lab2XYZ(L, a, b)(0, 0)
        Y = lab2XYZ(L, a, b)(1, 0)
        Z = lab2XYZ(L, a, b)(2, 0)
        
        '---CIEXYZからsRGBに変換
        r = XYZ2RGB(X, Y, Z)(0, 0)
        g = XYZ2RGB(X, Y, Z)(1, 0)
        b = XYZ2RGB(X, Y, Z)(2, 0)
           
        'セルを着色
        Cells(9, i).Interior.Color = RGB(r, g, b)
        
        i = i + 1
    Loop
    
End Sub

ImageMagickを使って複数のjpgをpdfにする

#include <stdio.h>
#include <stdlib.h>
#include <direct.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <windows.h>

int main(void) {
	
	/*jpgを保存するフォルダを作成し、jpgをpdfに変換、jpgを削除する*/
	char program[] = "D:\\Tools\\ImageMagick\\convert.exe";
	int ret;
	struct stat st;
	char work_dir[] = "D:\\work";
	char tmp_dir[] = "temp";
	char save_dir[256];
	char input_fn[] = "*.jpg";
	char output_fn[] = "all.pdf";
	char cmd[256];
	char output_fn_path[256];
	char del_jpg[256];
	int result;
	
	//jpgを保存するフォルダ
	sprintf(save_dir, "%s\\%s", work_dir, tmp_dir);
	
	//フォルダ存在チェック
	ret = stat(save_dir, &st);
	
	if(0 == ret){
		printf("%sは存在します。\n", save_dir);
	}
	else{
		printf("%sを作成します。\n", save_dir);
		_mkdir(save_dir);
	}
	
	//tempフォルダにあるjpgをpdfに変換
	sprintf(cmd, "%s %s\\%s %s\\%s", program, save_dir, input_fn, save_dir, output_fn);
	result = system(cmd);
	if (result == EXIT_SUCCESS){
        puts("pdfを作成しました。");
		//出力ファイルがあるフォルダを開く
//		system("explorer \"D:\\work\\temp\"");
		//出力ファイルがあるフォルダを開く ファイル選択状態にする
		sprintf(output_fn_path, "/select,%s\\%s", save_dir, output_fn);
		ShellExecute(NULL, NULL, "explorer.exe", output_fn_path, NULL, SW_SHOWNORMAL);
	}
	else{
        puts("pdf作成に失敗しました。");
    }
	
	//jpgファイルを削除する
	sprintf(del_jpg, "del %s\\%s", save_dir, input_fn);
	system(del_jpg);
	
    return 0;
	
}