2008年11月06日

RGBとHSVの相互変換[R]

Rのrimageパッケージはjpeg形式の画像をRに読み込ませて取り扱うのに便利なパッケージだが、色の表現形式としてはRGBのみの対応になっている。

RGB形式は確かに便利でわかりやすいが、用途によっては色相(Hue)、彩度(Saturation・Chroma)、明度(Brightness・Lightness・Value)の3つの成分を用いて色を指定するHSVという形式が便利なことがある。

便利なことがあるゆえに、様々な言語においてRGB⇔HSV相互変換のためのプログラム、アルゴリズムが書かれ、公開されているのだが、ちょっと探した感じではR用のものは見つけられなかった。(2008/11/7追記:はてブのコメント見ててもしやと思ったら変換用の関数があった。しかも名前がrgb2hsv。ちょっと打ち込んでみてたら分かったのに恥ずかしい…。ただそのままでは使えそうにないのでもうちょっと調べてまた書く予定。hsvをrgbにする関数は本当に無いかも?とりあえず以下のスクリプトのオブジェクト名は変更。)

ので作ってみた。

まずはRGB→HSV

#3次元配列から第3次元各要素の最大値を要素とする2次元行列を生成する
max.rgb <- function(z){
	y <- z[,,1]
	y[z[,,1]<=z[,,2] & z[,,3]<=z[,,2]] <- z[,,2][z[,,1]<=z[,,2] & z[,,3]<=z[,,2]]
	y[z[,,1]<=z[,,3] & z[,,2]<=z[,,3]] <- z[,,3][z[,,1]<=z[,,3] & z[,,2]<=z[,,3]]
	y
	}

#3次元(略)最小値を要素とする(略)
min.rgb <- function(z){
	y <- z[,,1]
	y[z[,,1]>=z[,,2] & z[,,3]>=z[,,2]] <- z[,,2][z[,,1]>=z[,,2] & z[,,3]>=z[,,2]]
	y[z[,,1]>=z[,,3] & z[,,2]>=z[,,3]] <- z[,,3][z[,,1]>=z[,,3] & z[,,2]>=z[,,3]]
	y
	}

#RGB→HSV変換
my.rgb2hsv <- function(x){
	#Vの計算
	hsV <- max.rgb(x)
	#Sの計算
	hSv <- (max.rgb(x)-min.rgb(x))/max.rgb(x)
	hSv[max.rgb(x)==0] <- 0
	#Hの計算
		#入れ物の用意
			Hsv <- x[,,1]
		#Rがmaxの部分のH
			Rmax <- 60*(x[,,2]-x[,,3])/(max.rgb(x)-min.rgb(x))
			Hsv[x[,,1]==max.rgb(x)] <- Rmax[x[,,1]==max.rgb(x)]
		#Gがmaxの部分のH
			Gmax <- 60*(x[,,3]-x[,,1])/(max.rgb(x)-min.rgb(x)) +120
			Hsv[x[,,2]==max.rgb(x)] <- Gmax[x[,,2]==max.rgb(x)]
		#Bがmaxの部分のH
			Bmax <- 60*(x[,,1]-x[,,2])/(max.rgb(x)-min.rgb(x)) +240
			Hsv[x[,,3]==max.rgb(x)] <- Bmax[x[,,3]==max.rgb(x)]
		#Hが定義されない部分
			Hsv[max.rgb(x)==min.rgb(x)] <- 0
		#値がマイナスの部分
			Hsv[Hsv<0] <- Hsv[Hsv<0]+360
	array(c(Hsv, hSv, hsV),dim=c(dim(Hsv),3))  #imagematrixオブジェクトと同様にH,S,Vの配列を出力
	}

最初にRGB成分の最大値・最小値のみを要素とする行列を作成する関数を定義し、これを変換用の関数中で用いている。なんかもうちょっと綺麗な書き方をしたい・できそうな気がしたけど思いつかない。

RGB→HSV変換関数の出力はimagematrixオブジェクトと呼ばれる配列と同様のものにした。H成分、S成分、V成分がそれぞれ一枚ずつ行列に格納される。

使うにはrimageパッケージのread.jpeg関数などで作成したimagematrixオブジェクトを突っ込むだけ。

photo <- read.jpeg("写真.jpeg")
photo.hsv <- rgb2hsv(photo)

ちょっと重いのが気になる。どこが原因なんだろ。

次にHSV→RGB変換。場合分けが多いのでちょっと大変。

#HSV→RGB変換
my.hsv2rgb <- function(x){
	H <- x[,,1]
	S <- x[,,2]
	V <- x[,,3]
	R <- array(rep(0,length(H)), dim=dim(H))
	G <- array(rep(0,length(H)), dim=dim(H))
	B <- array(rep(0,length(H)), dim=dim(H))
	#Hiという数値の計算
		Hi <- floor(H/60)
	#Fの計算
		rgbF <- H/60 - Hi
	#Mの計算
		M <- V*(1-S)
	#Nの計算
		N <- V*(1-S*rgbF)
	#Kの計算
		K <- V*(1-S*(1-rgbF))
	#Hiが0のとき
		R[Hi==0] <- V[Hi==0]
		G[Hi==0] <- K[Hi==0]
		B[Hi==0] <- M[Hi==0]
	#1のとき
		R[Hi==1] <- N[Hi==1]
		G[Hi==1] <- V[Hi==1]
		B[Hi==1] <- M[Hi==1]
	#2のとき
		R[Hi==2] <- M[Hi==2]
		G[Hi==2] <- V[Hi==2]
		B[Hi==2] <- K[Hi==2]
	#3のとき
		R[Hi==3] <- M[Hi==3]
		G[Hi==3] <- N[Hi==3]
		B[Hi==3] <- V[Hi==3]
	#4のとき
		R[Hi==4] <- K[Hi==4]
		G[Hi==4] <- M[Hi==4]
		B[Hi==4] <- V[Hi==4]
	#5のとき
		R[Hi==5] <- V[Hi==5]
		G[Hi==5] <- M[Hi==5]
		B[Hi==5] <- N[Hi==5]
	#S=0のとき
		R[S=0] <- V[S=0]
		G[S=0] <- V[S=0]
		B[S=0] <- V[S=0]
	imagematrix(array(c(R, G, B),dim=c(dim(R),3)))  #imagematrixオブジェクトとして出力
	}

出力はimagematrixオブジェクトなのでそのままプロットしておk。

plot(my.hsv2rgb(photo.hsv))

ついでにhsv形式の画像をプロットする関数なども作っておくと便利かもしれない。

#HSV形式のデータをプロットする関数
hsv.plot <- function(x){
	plot(my.hsv2rgb(x))
	}

まあちょっとだけ打ち込む文字数が減るだけなんですけど。

これだけなんだけど、肝心のhsvへ変換する関数が重いのがやっぱ気になる。使えないほどじゃないんだけど、rgbへの変換より明らかに重い。もうちょっと調べて改良しないとだめか。

posted by Rion778 at 20:09 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2008年11月02日

長野へ

研究室がクーラー工事でアレだというので、日帰りで大王わさび農場と諏訪大社へ行ってまいりました。

これは大王わさび園入り口で販売されていたわさび。3本1000円のものから1本5000円のものまでピンきり。 品種は田に吊るされていた札からするに長野23号だと思われるが確かではない。ぼちぼちの辛味。

入ってすぐのわさび田では洗浄作業を行っていた。わさびは通常2年をかけて栽培する。 放っておくと泥が蓄積し透水性が低下したり病気の原因となったりするのでわさび田は洗浄が必要。 てっきり最近は高圧の放水とかで機械的にやるものだと思ってたけど手作業だった。観光用か?面積を考えるとかなりの重労働。

かの有名なわさび像であります。その他にも謎としか言いようのないオブジェクトが点在しておりました。

定職直後と思われるわさび田。左手側が上流で、手前と奥に見える水路より中央の水路が低くなっており、 それぞれの畝間を水が流れる構造となっている。 わさびは自家中毒を起こす物質を放出するために畑などでは葉が大きくなるばかりで根茎が肥大しないが(いわゆる「葉わさび」となる)、 常に流れる水によりこれが洗浄される場合には根茎が肥大する(というのが定説だが詳細はよくわかっていないらしい)。

1年程度は経過していると思われる株。中央では花が咲いている。花は佃煮などにして食される。 割と見た目が汚いし病気や虫食いも見られるが、まあこんなものが普通。 きれいな葉などなかなか取れないので佃煮やらわさび漬けやらにするわけで。また、下流側であることも影響していると思われる。 湧き水というものは年中水温が一定なので、わさびが暑い夏や寒い冬を元気に乗り越えることを可能としているのだが、 それもわさび田の下流となれば気温の影響を受けやすくなり、また養分も少なくなる。そのため、 上等田というものは通常上流の水が湧いてくる付近に存在する。

これが大王わさび農場の名の由来となっている大王神社。「大王」とはなんぞという話だが、 その昔安曇野を治めていた八面大王という怪力無双の首領がおり、その大王の「胴体」が納められているのが大王神社だとかなんとか。 なんでも昔全国統一を目指す中央政権によって大王は打ち倒されたのだが、あまりに強い大王の復活を恐れ遺体をバラバラにしたのだとか。 なんと恐ろしい伝承か。

それで一通り大王わさび農場を見た後は諏訪大社へ。諏訪大社といえば御柱祭で有名。 まー諏訪神社も御柱祭的なものも全国各地にあるので有名もへったくれもないんだけど。私の地元の近く(といっても山を越える必要があるが) にもある。樹齢800年の大杉があるのでそれなりに有名なところ。

時間の都合もあり秋宮だけ。

神楽殿。巨大な注連縄と巨大な狛犬が印象的。注連縄はさすがに出雲大社の大注連縄が日本一だが、 狛犬は諏訪大社のものが日本一だそうだ。ちなみに諏訪大社の祭神の建御名方命は出雲大社に祀られる大国主の長男だそうで。 (詳しくはWikipediaなどでどうぞ)

といっても本来の祭神はミシャグチ神をはじめとする土着神だそうで。ミシャグチ神といえば蛇の姿をしているといわれる。 そのためにとぐろを巻いた注連縄が神楽殿の左右においてあるのだろう。ちなみに宝物殿にもあった。後ろに一本だけ (実際は幣拝殿を囲むように4本ある)御柱が見えるが、これがミシャグチ神をおろす依り代なんだそうだ。

これが上の写真で右側に見えていた御柱。賽銭箱の横には「御柱に賽銭を打ち込まないでください」と注意書きが。 打ち込むような輩がいるのかしら。ちなみに一之御柱と二之御柱は見えるが、三、四之御柱は幣拝殿裏側に設置されているためあまり見えない。

とりあえずこんなもんで。

今回の総走行距離は500km以上。低反発の座布団を駆使しても若干腰が…。 まー高速ばっかだったので数値が示すほどの大変さはなかったけど。

以上ぱぱっと紹介しましたがほかにもいくらかの写真がFlickrにアップしてありますのでよろしければご覧くださいまし。

posted by Rion778 at 14:12 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2008年10月21日

最小二乗法

なにか2つの数値が一組となったようなデータがあったとしよう(体重と身長、温度と湿度、国語の成績と数学の成績...etc)。このようなデータセットは散布図によって図に表すのが適当だろう。例えばこのように(以下Rを使って説明します)。

x <- c(1, 3, 4, 5, 7, 2, 8, 9, 10, 6)
y <- c(0, 4, 3, 7, 3, 4, 7, 9, 12, 8)
plot(x, y, pch=16)

ここでは、xとyと名づけられたそれぞれのデータセット(それぞれを変数と呼ぶ)の間には、片方が大きければもう片方が大きいといったような関係、正の相関があるように見る。

ここで、yをxの関数として予想できたら便利と思うかもしれない。

例えば、次のようなモデル(特にこれは一次式による線形モデルとか呼ぶ)を考えよう。

ここでaとbというのはまだ決まっていない定数で、xiというのが個々のデータとなっている(添え字のiには1, 2, 3...と言ったデータ番号が入る)。順番はともかく、いわゆる普通の一次式だ。yの頭に傘が付いているが、これはxiに各々対応するyの値を正確には表現できない(n個のyを完全にxの式で表示したければn-1次式が必要となる)ためで、yの推定値であることを主張している(ちなみに「ワイハット」と読もう)。

さて、まずは見た感じで決めてみよう。こんな感じでどうだろうか。

a=0, b=1としてみた。ためしにこの式の線を引いてみよう。

なかなかのものではありませんか。

まあ、場合によってはこれで十分かもしれない。しかしながら、例えばa=0.01, b=0.99と言ったような値と、今の値とどっちのが「より適切」であるかはどのように判断するべきだろうか。

ここで判断の指標には今引いた「線」から「点」までの距離、特に垂直方向の距離を使う。線をグラフに書き入れよう。

この垂直方向に引かれた線、その合計が最も少なくなるような線が、最も適切な線であるとするのはそれなりに合理的だろう。というわけで、この垂直方向の線が最も短くなるような線を探索することを目標としよう。

まず、「垂直方向の線」を数式で表現する。これはy軸方向におけるモデル式とデータポイントの間の距離であり、誤差、あるいは残差と呼ばれる。残差をdという文字で表すことにする。あるyの値yiに対するdは添え字をつけてdiとする。つまり、

そうして問題はこのdiの絶対値の総和を最小にしようということだった。最小化するためにいじれるのはaとbという値だけで、それらはy[ハット]を表現するモデル式の中に含まれている。ところで「絶対値の」というのは関心があるのが差ではなく「距離」だからである。しかし、絶対値を含む数式というのは非常に計算が難しくなるため、ここではdを二乗することにより符号を正に揃えるという戦法を採ろう。

よって、目的は次の式を最小化するaとbの値を求めることに置き換えられた。で、置き換えられたついでにちょいちょいと変形してaとbを表に引きずり出してやろう。

展開して多少見づらくなったかもしれないが、ここで自由に変更できる文字はaとbだけであることを思い出すと、この式はaとbについてそれぞれ2次式になっていることに気づく。そして2次式として眺めれば、それぞれ下に凸であるということにも気づくだろう。

下に凸な2次式の最小化といえば…「微分して0とおく」というやつ。つまり、変化率が最小になるような点が極値であり、最小値となる。

「微分して」と言ってみたがここでは変数がaとbの2つある。変数が複数ある場合、偏微分という方法を使う。偏微分と聞くと構える人が居るかもしれないが、注目する変数以外を定数とみなして普通に微分するだけなので心配することはない。要するに、bを定数とみなして微分(aについて偏微分)、aを定数とみなして微分(bについて偏微分)の2回の微分をする。すると、次のように2つの式が得られる。

(ここの∂というのは偏微分であることを示す微分記号で、デルだとかラウンド・ディーだとか呼ばれる)

aを変化させたときの変化率、bを変化させたときの変化率というわけなので、こいつらを両方0と置く。

そしてこの2式を連立させ、aとbについて解けば問題は終了。要するに「解の公式」を求めるのが目的なのだ。以下は面倒なだけで特にこれと言って難しい部分もないので読み飛ばし推奨。

まず、式(1)より次式を得る。

式(2)へ代入。

bについて解く。

見やすい感じにする。

次に式(2)より次式を得る。

式(1)へ代入。

aについて解く。

見やすい感じに。

さて、これでaとbに関する「解の公式」が得られた。並べてみよう。


特にこの「形」が意味するものも無いのでさーっと眺めてもらえば問題ないのだが、とりあえず生データのみで答えが出せるというのがポイント。そして実際に計算すべき量は5つしかないのもポイント(よく観察!)。

まあとにかくこの辺は統計ソフトがちょいちょいとやってくれるので、流れだけ把握しておけばよいかなとか思う。モデルは一次式とも限らないわけだし(といっても根拠がなければ一次式にすべきといわれる)。

それで、このように残差の二乗を最小化することにより最適な(とは限らないのだが)モデルを求める手法を最小二乗法と呼ぶ。

最後に最小二乗法により得られたデータ直線を最初の方のグラフに重ねてプロットしてみよう。

まあ最初の「見た感じモデル式」もボチボチではありませんか。手元に電卓とかないときはこういうのも重要ですよ。

ちなみに、Rだと最小二乗法による回帰はこんな感じでやります。

lm(y~x)

これだけで切片と傾きが!で、もう少し詳しい情報が欲しければ

summary(lm(y~x))

と。傾きや切片の有意性、決定係数などといった情報を与えてくれます。

ちなみに最小二乗法は1805年、ルジャンドル先生が発表したわけだけれど、その10年前にはガウス先生が発見していたといういわくつき(?)の方法。ガウス先生はどうもメモに書き付けて満足して発表しないといったことがしばしばあったそうで、先生がちゃんと論文書いてたら数学の発展は何十年も早かったろうとかなんとか。

まあ、要するに論文はちゃんと書きましょうという話でした。

posted by Rion778 at 19:51 | Comment(4) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2008年10月05日

RICHO R8購入

とうとうデジカメを新調してしまった。RICHOのR8というヤツ。

ネットで探したけど安いのがなかったので、某ヤの付く電機屋で買うことに。店頭展示品で31,000円(3,000ポイント付き) ってのがあったので、他の量販店、カメラ屋を確認してから購入。値切ろうとしたけど全く駄目だった。展示品は無理とか何とか。残念。

でもまあポイントを当日使えるようにしてもらえたので、4GのSDと5年保障つけて現金出費合計3万3千円。 もう今月はモヤシばかり食べるしかないな。あと本も買わない。

とりあえずナシが手元にあったのでためし撮り。

ナシ

ついでに刺さってる爪楊枝で接写のテスト。

つまようじ

アルコールとカフェインの中毒症状で震えのとまらない(嘘です)僕の手で適当に撮影してもこんな感じ。 三脚を使えばもっと綺麗に撮れる?

posted by Rion778 at 17:16 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2008年09月16日

えびふりゃー

天ぷらなべを買ったので揚げ物練習中。

今日はばーちゃんに電話で作り方を聞いたエビフライ。殻のむいたやつを買おうかと思ったけど、 からつきのやつor衣までついたやつという選択肢しかなかったため、結局殻つきブラックタイガーを15匹300円で買ってきて調理。 手がエビ臭い…

思った程度にはうまくできた。思った程度には失敗した。

腹に切り込みを入れて伸ばさないと曲がるといわれてたのでそれっぽい処理をしたはずなんだけど結局曲がった。 もっと伸ばさないとだめなのか。

エビフライ

ばーちゃん曰く小麦粉を2回つけるのがポイントらしい。理由はなんとなく。

小麦粉→卵→小麦粉→卵→パン粉

このフローをばーちゃん語から汲み取るのにどれだけ苦労したか…

さっきばーちゃんから「えび天はうまくできましたか」とメールがきた。

本気なのか冗談なのかわからないのがポイント!

posted by Rion778 at 23:35 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする
2008年09月07日

ふりっかー

存在を知ってからも、アカウントを取ってからも相当に放置していたFlickrに写真をアップロードなどしてみる。案の定面白いのでもっと写真撮ろう。もうちっと良いデジカメが欲しいなー。

オニヤンマ

こいつは実家の玄関の前にとまっていたオニヤンマ。Flickrにアップしてから貼り付けた。便利です。僕の携帯(SO905iCS)で適当に撮影したんですが、このオニヤンマをめぐってばーちゃんと母上がこんなやりとりを。

ばば「うーん。こう、後ろから撮らないと…」

はは「欲しいの?トンボの写真欲しいの?撮ってあげようか?」

ばば「別に欲しいわけじゃ…」

はは「欲しいんでしょ?デジカメとってくるから」

ばば「べ、べつにいらないったら!でもどうしても撮りたいっていうんならもらってあげなくもないんだからね!」(一部補正)

デジカメを持ってくる母親。そして飛び立つオニヤンマ。

実家はいたるところに虫が落ちている(例えば廊下にカブトムシが落ちてたり)ようなところ。どうやらでかいハチの巣が家のどこかにあるそうで、「ハチが大量に出入りしてるよ!」って隣の家に言われたらしい。それに関してははとばばがこんなやりとりを。

ばば「お金もないし別に危なくないだろうから『もしどうしても嫌だったら言ってね』って言っておいた。」

はは「まあスズメバチじゃなければ大丈夫じゃない」

ばば「ハチの巣といえばこの前裏庭にすごく大きなハチの巣があって、おじいさんが生きていたころだけど」

はは「…それって『この前』じゃないよね」

ばば「…20年くらい前?」

冗談なのか本気なのか分からないのがポイント!

本当は朝ご飯だけ食べたら帰ろうかと思ったけど、居心地がいいので晩ご飯まで食べたら帰ろうかと考え中。

posted by Rion778 at 13:50 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2008年08月26日

部分と全体を比較する

全体の比率を期待値にし、部分部分の比率がどれだけ期待値から食い違っているのかを観察するのが独立性の検定。

部分内のバラつき(グループ内の変動)を、部分を一つのデータと見たときの全体のバラつき(グループ間の変動)と比較し、「部分の塊のバラつき」が「部分内のバラつき」よりどれだけ大きいかを観察するのが分散分析。

似ているようだけど考え方の順序は逆だ。片方は全体から部分へ、もう片方は部分から全体へ。

「全体」から外れているからそれぞれは異なる。

「部分」よりもばらついているからそれぞれはばらついている。

「全体」と「部分」を同一のものさしで測っていいのか、ということは少し気になる。

posted by Rion778 at 19:27 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2008年07月27日

三角関数などの根を小さい方から順に求める[Maxima]

前回RでやったのをMaximaで。

root:[0,0,0,0,0,0,0,0,0,0];
f(x):=cos(x)*(x*3);
myfindroot(f):=block([j:0,a:0.1,k:0],
 loop,k:k+1,
 for i:1 while f(j)*f(j+a)>0 do(a:a+0.1),
 root[k]:find_root(f(x),x,j,j+a),
 j:root[k]+0.1,
 if k=10 then return(root),
 go(loop));

前回は関数のエラーを利用していてスマートじゃなかったのでちゃんと符号を見て関数を使うようにした。 まーでも範囲を0.1ずつ増やしているので、0.1の範囲にいくつも根があると見逃してしまうだろう。もっと上手い方法はないものかしら。

しかしMaximaは使い方がよくわからないなー。一冊くらい本を買って勉強しないとだめかな。 Rみたいに任意の長さのベクトルをパッと作る方法は無いんだろうか。あとエディタからパッと命令を送る方法は無いんだろうか。 エディタをいじらないと駄目な気がしているけどどうなんだろう。やっぱRが一番使いやすいなー。

posted by Rion778 at 21:50 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。