2013年8月23日金曜日

Fourier Slice Theorem (retry)

The other day, I tried implementing Fourier Slice Theorem but couldn't reconstruct the projecttion line from the slice of the Fourier domain.  I thought the slicing operation was not accurate enought, and I tried some interpolation when slicing this time.

to be continue...

2013年8月19日月曜日

Refocus by integral projection

リフォーカスのgif画像を作ってみた。

リフォーカスのおおざっぱな原理はこちら。自作レンダリングソフトを使ってSubaperture Imageを取得してそれを足し合わせる。


フーリエ変換でやる方は挫折。挫折した理由はこちらをご覧ください。

Fourier Slice Theorem (2D)

Light Field Cameraの勉強をしているとFourier Slice Theoremというのが出てくる。Lytro創始者Ren NgのFourier Slice PhotographyではN次元まで拡張されている。4D light fieldの2D Sliceはうまく行かなかったので、まずは2次元でちゃんと出来るか確認してみた。
図はRen Ng氏のドクター論文より。



まずはProjectionを行う。図の左側だ。基底変換をした後に積分をすればよい。
これを実装しようと思ったが、アホな頭には思った以上に難しい。例によってツールを使う。どうやらこのような変換はRadon Transformとして知られているおりCTなどのトモグラフィに応用されているようだ。Pythonで使う場合はskit-imageを使う。このライブラリで変換を行うと2次元画像が出力される。が、今回は1次元(1つの角度)だけ取り出して描画する。これでProjection演算は完了。図は元画像を10°傾けてプロジェクションしたもの。
元画像は0と255の2値からなる画像なので、実はこの変換の結果は角度に関係ない。


続いて、図をフーリエ変換したあとにSliceしてそれを逆フーリエ変換する。図では→、↓、←、という変換経路。FFT→Rotation→Slice→iFFTというふうに実装してみた。これの変換をすれば先ほどと同じ出力が得られる、というのがFourier Slice Theoremである。スライスする部分は自分で実装。

まず下図は0°でスライスした結果。FFTの結果はパワー(絶対値の2乗)を画像として載せている。てっぺんが少し平らになったがなかなか良く再現出来ている。



続いて10°の結果を載せる。だいぶ形がおかしい。


別の現実的な画像でやってみる。まず0°。スケールが違うので少し違って見えなくもないが、再現出来ていると言えよう。


続いて10°。


んー。似てない。
多分Sliceするときに精度が落ちてしまったのだろう。4D Light Fieldでも同じ原因で出来なかったのだろう。しかし、こんなに変わってしまうとは。なにか工夫すればうまくいくのか…

追記:
どうやら調べてみると、離散的になデータをグリッドで扱うところに誤差があるっぽい。
同じようなことをやろうとした人たちがいるようだ。
https://groups.google.com/forum/#!topic/comp.soft-sys.matlab/Fe9Wf8Slhzs

最後の方では補完を試みたけど、全然ダメだったとか、
実際はこの方法使われてないとか書いてある。


画像の回転

画像をFFTした結果の回転を行うプログラムを組む。

画像をスライスした時の値(複素数)が欲しかったので組み始めた。回転なんてちょろいと思ったら、意外と組めなかった。
Scikit-imageにもrotateという関数が準備されていたが、回転後の画像が切れてしまう。

今回は以下のページを参考にした。
http://homepage2.nifty.com/tsugu/sotuken/rotation/

Pythonでそのまま組むと速度がかなり遅い。いままで何も考えてなかったけど回転って計算負荷あるんやなー。今回は実験なので速度のことは(そんなに)気にしない。




また使うか分からんけど、一応ソースも貼っておこう。

import numpy as np

def rotate_img(src, a):
 #src: numpy array
 #rotation angle in degree
 #returns rotated image

 a = np.pi/180. * a 
 vcos = np.cos(a)
 vsin = np.sin(a)

 #calc dist image size
 wsrc = src.shape[1]
 hsrc = src.shape[0]
 wdist = int(np.abs(wsrc*vcos) + np.abs(hsrc*vsin) + 0.5)
 hdist = int(np.abs(wsrc*vsin) + np.abs(hsrc*vcos) + 0.5)

 #allocate memory
 dist = np.zeros((hdist, wdist), dtype=src.dtype)

 #center of the images
 cxsrc = wsrc / 2
 cysrc = hsrc / 2
 cxdist = wdist / 2
 cydist = hdist / 2

 for y2 in xrange(hdist):
  for x2 in xrange(wdist):
   x1 = ((x2-cxdist)*vcos - (y2-cydist)*vsin) + cxsrc
   y1 = ((x2-cxdist)*vsin + (y2-cydist)*vcos) + cysrc

   if x1 >= 0 and x1 < wsrc and y1 >= 0 and y1 < hsrc:
    dist[y2,x2] = src[y1,x1]
 return dist

def slice_img(src, a):
 #src: numpy array
 #rotation angle in degree
 #return sliced 1d array

 a = np.pi/180. * a 
 vcos = np.cos(a)
 vsin = np.sin(a)

 #calc dist image size
 wsrc = src.shape[1]
 hsrc = src.shape[0]
 wdist = int(np.abs(wsrc*vcos) + np.abs(hsrc*vsin) + 0.5)
 hdist = int(np.abs(wsrc*vsin) + np.abs(hsrc*vcos) + 0.5)

 #allocate memory
 dist = np.zeros((wdist), dtype=src.dtype)

 #center of the images
 cxsrc = wsrc / 2
 cysrc = hsrc / 2
 cxdist = wdist / 2
 cydist = hdist / 2

 for x2 in xrange(wdist):
  x1 = ((x2-cxdist)*vcos) + cxsrc #y2 = cydist
  y1 = ((x2-cxdist)*vsin) + cysrc

  if x1 >= 0 and x1 < wsrc and y1 >= 0 and y1 < hsrc:
   dist[x2] = src[y1,x1]
 return dist


if __name__ == "__main__":
 import Image
 import matplotlib.pyplot as plt
 img = Image.open("_image0-0.png").convert('L')

 a = 10.

 src = np.array(img)
 dist = rotate_img(src, a)
 sliced = slice_img(src, a)

 #dist = Image.fromarray(np.uint8(dist))
 
 plt.figure()
 plt.subplot(131)
 plt.title("original")
 plt.imshow(src, cmap=plt.cm.Greys_r)
 
 plt.subplot(132)
 plt.title("rotated")
 plt.imshow(dist, cmap=plt.cm.Greys_r)

 plt.subplot(133)
 plt.plot(np.arange(len(sliced)), sliced)
 plt.title("sliced at y = N/2")
 plt.show()



情熱大陸がおもしろくなくなった。

以前は毎週情熱大陸を楽しみにしていた。
それがいつからかは分からんけど、徐々に見なくなった。
さっき葉加瀬太郎のバイオリンが流れてきて、今日は情熱大陸の日だと知った。(でも見んかった)

どうも最近は芸能人が多い気がする。マイナーでもいいので第一線で活躍している人を紹介してほしい。それじゃ視聴率が取れないということなんやろうけど、アサヒがスポンサーならわかってくれるんじゃないか。昔は情熱大陸を見たあとはモチベーションが上がって俺もがんばるぜモードに入っていた。

過去の放送を見ると芸能人以外も意外に取り上げている。どうも面白くなくなった原因はそこにあるのではないらしい。プロデューサーが変わったのか、自分の趣向が変わったのか。

2013年8月18日日曜日

boostのublasのベクトル計算メモ

Pythonではnumpyを使わせて頂いてだいぶ作業が楽だ。
C++でベクトル計算するときは、boostにublasという行列ライブラリがあるのでそれを使って計算ができる。vector宣言がSTLと被るので注意が必要。

ベクトルや行列同士の積は関数で計算するようだ。
以下使い方のメモ

#include <iostream>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace std;

int main() {
 using namespace boost::numeric;

 ublas::vector<double> v1(3);
 ublas::vector<double> v2(3);
 ublas::matrix<double> m(3,3);

 for(unsigned i = 0; i < v1.size(); i++) {
  v1[i] = i;
  v2[i] = i+3;
 }
 
 for (unsigned i = 0; i < m.size1(); i++)
        for (unsigned j = 0; j < m.size2(); j++)
            m(i, j) = 3 * i + j;
       

 cout << "v1: " << v1 << endl;
 cout << "v2: " << v2 << endl;

 cout << "add:   "  << v1 + v2  << endl;
 cout << "sub:   "  << v1 - v2  << endl;
 cout << "multi: " << 5.0 * v1 << endl;
 cout << "div:   "  << v1 / 5.0 << endl;

 //product operation
 cout << "inner: "  << ublas::inner_prod(v1,v2) << endl;
 cout << "outer: "  << ublas::outer_prod(v1,v2) << endl;

 //with matrix
 cout << "prod: " << ublas::prod(m,v1) << endl;

}


ヘビでもわかるライトフィールドカメラの原理 その3.2 データと分解能

その2.4でLytroのカメラの構造について説明したが、今回はどのようなデータが実際に取得されているかについて簡単に説明する。




左の図はライトフィールドカメラの構造を2次元で表したもので、センサーは1次元で表されているが、実際はもちろん2次元。そこに4次元のデータ(4D Light Field)を記録する訳である。それを示したのが右の図である。センサーは(x, y)の座標を持ち、その中に更に(u, v)の座標がある。

(x, y)はマイクロレンズの位置によって決まり、(u, v)はその下のピクセルの座標であり、絞りのどの位置から来た光かが記録される。これらは当然ながら離散的であり、(x, y)の分割数はマイクロレンズの数、(u, v)の分割数はマイクロレンズ下のピクセルの数によって決まる。
つまり、空間分解能はマイクロレンズの数、方向分解能はマイクロレンズ下のピクセル数によって決まる。

これらのデータはデジタルデータなので自由に事後処理することができる。たとえば(x, y)と(u, v)の関係を入れ替えてみる。



すると(u, v)を指定したとき、その座標に形成される画像は前回説明したsubaperture imageになるということが分かる。(絞りの場所を指定した時の画像なので当然と言えば当然だが…)


事後処理についてはたぶんいろんなアルゴリズムがあってここでは説明できないので割愛する。



ヘビでもわかるライトフィールドカメラの原理 その3.1 フォーカス再考

前回までにライトフィールドカメラの原理と構成を説明して、ライトフィールドの取得方法がわかった。取得したライトフィールドをデジタルで処理する訳であるが、それを理解するためにリフォーカスするということについてわかりやすく概念を図示しておこう。




よりレンズに近い1点から出た赤い光線はよりレンズから遠いところで、
よりレンズから遠い1点から出た青い光線はよりレンズから近いところへ結像する。図ではどちらも撮像面で結像していない。つまり、赤い光線はも青い光線もピンぼけしている状態である。

ライトフィールドカメラで光線が記録できれば絞りのある1点からくる光線のみを取り出すことができる。図のように絞りの上側から来る光線と下側から来る光線を取り出してみる。これらの画像は「ボケていない視差のある画像」と考えることができる。Ngの論文ではsubaperture imageなどと読んでいる。




例えば赤い点にフォーカスしたいとき、Image1を上へ、Image2を下へ赤い点が一致するように動かせばよい。このとき青い点は結像していないのでボケた像になる。逆に青い点にフォーカスしたいときは、反対方向へ動かせば良い。実際には画像は2つだけでなく複数のSubaperture imageを取得することができるので、それらも足し合わせることになる。

つまり、フォーカスするということは視差のある画像をずらしながら足していく行為だといえる。

2013年8月14日水曜日

numpy de 2d fft

numpyで画像のFFTを行う。

まずグレースケールで基本を抑える。
import numpy as np
import Image

if __name__ == "__main__":
 img = Image.open("pngs/image0-0.png").convert('L') #open and convert to grayscale
 img = np.array(img)  #convert to numpy array
 
 ft = np.fft.fft2(img)  #fft
 Pow = np.abs(ft)**2  #get power spectrum
 Pow = np.log10(Pow)  #adjust scale to log

 Pmax = np.max(Pow)  #ajust to 8-bit scale
 Pow = Pow / Pmax * 255 

 pow_img = Image.fromarray(np.uint8(np.fft.fftshift(Pow))) #convert to Image
       #fftshift makes DC componets to the origin

 pow_img.show()   #display








2013年8月13日火曜日

Sublime tex2 でPythonのBuild設定(備忘録)

Sublime textでは"cmd+shft+b"でBuild実行してくれて、画面したに実行結果がでる。
でもPythonを使うとき、初期設定ではimport Errorが起きる。

そんなときはPackages/Python/Python.sublime_buildに
 "path": "/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin",

を加えるとうまく実行してくれた。めでたしめでたし。

2013年8月12日月曜日

Perspective Sihft like Lytro

I change the last code and made a gif of perspective shift.

The rendering algorithm is simple. Just trace rays coming from the same aperture position.
The aperture is devided in to 4x4. Here is a nice fig of explanation.



The resulting gif is a bit ackward, but would be smooth as the number of division gets bigger.




2013年8月9日金曜日

PILで画素を操作、保存するとき。

画像を加工するときにOpenCVとかライブラリを使えばいろいろな変換ができるけど、画素単位の操作やったり、1つひとつ動作を確認しながらやりたいとき用のメモ。
PythonのPILとNumpyを使う。


おまじない
import numpy as np
import Image

画像を開く
img = Image.open("tmp.png")
グレースケールで開く
img = Image.open("tmp.png").convert('L')
新しい画像の生成。
img = Image.new('RGB', (5,2)) #(w,h)

Imageからnumpy.arrayへの変換
a = np.array(img)

numpy.arrayからImageへの変換
img2 = Image.fromarray(a)
intに変換するときは
img2 = Image.fromarray(np.uint8(a))


ImageをPILのビューアーを使って表示
img.show()


保存
img.save("tmp.png")





2013年8月6日火曜日

Pattern Recognition and Machine Learning for Dummy section 11.1

PRML 11章をやろうとしたが、11.1の最初の数式でいきなりつまずいた。
電車の中でどういう意味か考えた。

zのサンプル密度はp(z)によって決まる。ほんで今の場合は一様なのでp(z)=1となる。
この時y=f(z)は図のようになるが、zが一様だとするとよく出る値とそうでない値がある。要するに傾きによってyの密度が決まるので式(11.5)のようになる。

なーんだ。

九寨溝、黄龍の行き方(旅の日程の巻)

7月末に4泊5日で九寨溝、松藩、黄龍の旅行にいってきた。
もちろん今回もすべて個人で手配した。
前回のマチュピチュでは高山病でゲロゲロだったので、
今回(黄龍で海抜3700mmくらい)は余裕をもってゆったりとした旅を志した。
しかし、下調べをしてなさ過ぎて色々と苦労した。
その苦労した体験を誰かと共有できれば幸いである。

旅の日程はおおまかに

  1. 成田(9時発)→上海→重慶→九寨溝空港(夜9時くらい)
    • ホテルへ(タクシー)
  2. 九寨溝堪能
  3. 松藩へ移動(バス。7:30〜10時くらい)
    • 松藩見物
    • 黄龍へ移動(タクシー)
  4. 黄龍国家級風景名勝区、見学
  5. 九寨溝空港(7:50)→北京→成田(18:00)
観光のほとんどは九寨溝と黄龍の見学に費やして、あとは移動。
トラブルもいろいろあったけど、とりあえず美しい風景をお楽しみくだいさ。


①成田(9時発)→上海→重慶→九寨溝空港(夜9時くらい)

 まずは飛行機で上海へ。乗り換え組は別コースで、パネルを持った人に付いて行く。英語はほぼ通じない。重慶に到着するとチェックインできない人たちが居座っていた。時間があったのでそこで食事をとる。メニューの写真を見ながら48元のものを注文。

なぜか余計に2元取られた。抗議すると、どうやら箸代が2元らしい(写真下の袋が箸)。ご飯を食べてまったりしていると、どうやら国内飛行機は別のところでチェックインするとのこと。そこまで徒歩で行く。





九寨溝空港につくとあたりはもう暗い。ここからタクシー2時間のって九寨溝の近くにとってあるホテルへ。


地球の歩き方によると夜間なら200元で九寨溝まで連れてってくれるらしい。時間が9時くらいだったのでバスは諦めてタクシーを探す。タクシーもそんなに見当たらない。客引きがよってきて、看板の一部を指さす。看板には九寨溝まで300元と書いてあった。これは公式のものだろうか、ちょっと疑わしかったけど、夜も更けてきたので300元で乗ることにした。片言の英語だったがメーターだからダイジョブといっている。そして、夜の道を走りだした…ここで中国を体験することになる。


片言の英語で「日本はサッカー強いねー」なんていってるが、一向にメーターをつける気配がないので、「メーター、メーター!」というとやっとメーターをつけた。メーターを早く付けないと自分たちが不利なのに、なぜ付けない。と思っていた。
しばらくして、車は泊まった。まだ15分くらいしか走ってない。そこにはもう1台車が止まっていて、運転手はその車の主と思われる人と会話を始めた。
全然再発進する気配がないので、

「どうした、早く行けや。」

と促すと、友達が来るからそいつの車に乗れと…

こっちも強く出たいが、周り真っ暗でしかもここはどこな状態。置き去りにされたら高地で野宿になってしまうので強くは出れない。

しばらくすると友達がきた。車のグレードは明らかに落ちていた。そいつに乗り込んで走っていく。ほんとに連れてってくれるのだろうか、とずっと不安に思っていたが、2時間くらいでホテルへ連れて行ってくれた。どうやら下請けの運転手のようだった。
空港で捕まえたタクシーは始めから下請けに渡すつもりだったのでメーターを上げなかったのだ。


つづく








九寨溝、黄龍の行き方(反省の編)

前の記事で旅の内容をダラダラと書いた。
今回は中国旅行において気になった点をまとめておく。


  1. 現金は余裕をもって換金しておく。
  2. 予約はExpediaとかを通して行う。
  3. 基本的にみんな悪いやつじゃない。
  4. ご飯はあまり美味しくない。
1. 現金は余裕をもって換金しておく。
 海外旅行の基本中の基本かもしれない、中国ではVISAが使えないのでなるべく現金を用意していく。ホテルやATMに「VISA」なんて大々的に書いてあっても使えない。ここにもそう書いてある。今回は優しい日本人の方と偶然遭遇したので難を逃れることができた。


2. 予約はExpediaとかを通して行う。
 今回は1つはBooking.comを通して、もう1つはChinaHotelHot.comを通して行った。
後者では、地球の歩き方に載っている松潘華龍山荘というところに泊まった。
どちらもカウンターでは英語は通じなかったが、前者の方がかなりまともなホテルだった。Expediaとかに登録されているホテルなので、ホテルの人が慣れている。後者は英語で予約できたが、中国人向け。また、世界中の人が選んでいるので情報収集の観点でも役に立つ。例えば、黄龍の近にホテルを取ろうと思ったが、Expediaとかで探しても出てこない。それは、そこに世界の需要がないからである。みんな近傍の町からそこへアクセスするし、町には旅行者向けにそういう移動ができるシステムが存在する。

3. 基本的にみんな悪いやつじゃない。
 ぼったくろうとする奴もおるけど、本気じゃない。インドの場合は本気。インドのやばさを体験しておいてよかった。
 英語は通じなかったが、基本的にみんな暇そうなんで助けようとしてくれた。

4. ご飯はあまり美味しくない。

 上海ではその辺の店でもすごく美味しかったが、今回はそうではない。カップラーメンもっていったが、すぐに食べてしまった。あ、カップラーメンもってくときはフォークも忘れずに。ホテルに箸やフォークはありません(特によ。







2013年8月4日日曜日

path tracing でレンダリング(c++) その2DOF (Depth Of Field)

先日作ったレンダラにレンズを加えてDOF(Depth Of Field)を再現してみた。
センサーとレンズの位置関係を決めるとピクセルと共役な点(集光する点)が求まる。
はじめにその共役点を求めて、絞りの各座標からそこへ向けて光線を飛ばす。




センサーの分割だけでなく絞りも分割するので、ピンホールカメラに対して計算時間がかなり増えてしまう。MacbookAirだとなかなか時間がかかる(
並列化もしてないし。)ファンも結構回る。






上の図は瞳分割4x4でに50sample/direction。下の図は瞳分割なしの50*16 sample/pixel

ところで家のデスクトップはCore i5の2.67GHzだが、こっちのほうが遅い。
理由は多分32bitのGCCを使ったせいだろう。こんなにも速度が違うのか、初めて知った。

/*
pathtrace_basic: basic path tracing from smallpt
pathtrace_dof: introducing a camera model.
this subaperture
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
//#define drand48() (double(rand()) / RAND_MAX)

class Vec {
public:
 double x, y, z;

 // constructor
 Vec(double x_ = 0.0, double y_ = 0.0, double z_ = 0.0) {x = x_; y = y_; z = z_;} 


 //operator overloading
 Vec operator+(const Vec &b) { return Vec(x+b.x, y+b.y, z+b.z); }
 Vec operator-(const Vec &b) { return Vec(x-b.x, y-b.y, z-b.z); }
 Vec operator*(double b) { return Vec(b*x, b*y, b*z); }
 Vec operator%(const Vec &b) { return Vec(y*b.z - z*b.y, z*b.x-x*b.z, x*b.y-y*b.x);} // cross product
 Vec mult(const Vec &b) { return Vec(x*b.x, y*b.y, z*b.z); }
 Vec& norm() { return *this = *this * (1/sqrt(x*x + y*y + z*z)); }
 double dot(const Vec &b) { return x*b.x + y*b.y + z*b.z; }
};

class Ray {
public:
 Vec o, d;
 Ray(Vec o_, Vec d_) : o(o_), d(d_) {} // constructor
};

class Plane2D {
public:
 double x, y;
 int xres, yres;
 Plane2D(double x_, double y_, int xres_, int yres_): 
 x(x_), y(y_), xres(xres_), yres(yres_) {}
};

class Sensor : public Plane2D {
public:
 Vec *c; //color array
 Sensor(double x_, double y_, int xres_, int yres_): Plane2D(x_, y_, xres_, yres_) {
  c = new Vec[xres*yres];
 }
 ~Sensor(){
  delete c;
 }
};

class Lens : public Plane2D {
public:
 double f;
 Lens(double f_, double x_, double y_, int xres_, int yres_): Plane2D(x_, y_, xres_, yres_) {
  f = f_;
 } 
};

class Camera {
public:
 Sensor* snsr;
 Lens* lens;
 double a;   //distance from lens to sensor, suppose lens and snsr are parallel.
 double beta;  //magnification
 Vec pos;   //center of the lens
 Vec dir;
 Vec up;
 Vec hrz;

 //constructor
 Camera(Sensor& s, Lens& l, double a_, Vec pos_, Vec dir_, Vec up_) {
  snsr = &s; lens = &l;
  a = a_; pos = pos_; dir = dir_; up = up_;
  hrz = (up%dir).norm();
  beta = l.f / (a-l.f);
 }

 Vec conj(int xi, int yi);
 Vec lens_coord(int ui, int vi);
};

Vec Camera::conj(int xi, int yi) {
 //input coordinate on the sensor
 //returns coordinate of the conjugate position
 Vec d = hrz * (0.5 - (0.5+xi)/snsr->xres) * snsr->x + up * (0.5 - (0.5+yi)/snsr->yres) * (-snsr->y) + dir * a; //from pix to the lens' center
 return pos + d * beta;
}
Vec Camera::lens_coord(int ui, int vi) { 
 return pos + hrz*((0.5+ui)/lens->xres - 0.5)*lens->x + up*((0.5+vi)/lens->yres - 0.5) * (-lens->y);
}



enum Refl_type { DIFF, SPEC, REFR};

class Sphere {
public:
 double rad;
 Vec pos, emi, color;
 Refl_type refl;

 //constructor
 Sphere(double rad_, Vec pos_, Vec emi_, Vec color_, Refl_type refl_):
  rad(rad_), pos(pos_), emi(emi_), color(color_), refl(refl_) {}

 double intersect(const Ray &r) {
  //returns distance, 0 if missed
  Vec op = pos - r.o;
  double t, eps=1e-4;
  double b = op.dot(r.d);
  double det = b*b - op.dot(op) + rad*rad;
  if (det < 0) {
   return 0;
  } else { 
   det = sqrt(det);
  }

  return (t = b-det) > eps ? t : ((t = b + det) > eps ? t : 0);
 }
};

Sphere spheres[] = {//Scene: 
 //     radius, position,               emission,    color,     material 
   Sphere(1e5,    Vec(-1e5-50, 0.0, 0.0), Vec(),     Vec(.75,.25,.25), DIFF),//Left 
  Sphere(1e5,    Vec( 1e5+50, 0.0, 0.0), Vec(),     Vec(.25,.25,.75), DIFF),//Rght 
  Sphere(1e5,    Vec(0.0, 0.0, -1e5-50), Vec(),     Vec(.75,.75,.75), DIFF),//Back 
  //Sphere(1e5,    Vec(0.0, 0.0,  1e5+1000),Vec(),     Vec(),            DIFF),//Frnt 
  Sphere(1e5,    Vec(0.0, -1e5-50, 0.0), Vec(),     Vec(.75,.75,.75), DIFF),//Botm 
  Sphere(1e5,    Vec(0.0,  1e5+50, 0.0), Vec(),     Vec(.75,.75,.0), DIFF),//Top 
  Sphere(16.5,   Vec(-20, -50+16.5, -50), Vec(),     Vec(.99,0.0,.0),  DIFF),//
  Sphere(16.5,   Vec(  0, -50+16.5, 0), Vec(),     Vec(1,1,1)*.999,  REFR),//Glas 
 Sphere(16.5,   Vec( 20, -50+16.5, 50), Vec(),     Vec(.99, .99, .0),  DIFF),//
  Sphere(300,    Vec(0.0, 349., 50.0),   Vec(6,6,6), Vec(),    DIFF) //Lite 
 }; 

inline double clamp(double x) {
 return x < 0 ? 0 : x > 1 ? 1 : x;
}

inline int toInt(double x) {
 return int(pow(clamp(x), 1/2.2)*255+.5);
}

inline bool intersect(const Ray &r, double &t, int &id) {
 //loop over the spheres,
 //and store the closest distance and the id.
 int n = sizeof(spheres) / sizeof(Sphere);
 double d, inf = t = 1e20;

 for (int i = 0; i < n; i++) {
  if ((d=spheres[i].intersect(r))>1e-4 && d < t) {
   t = d;
   id = i;
  }
 }
 return t<inf;
}

Vec radiance(const Ray &r_, int depth_) {
 double t; //distance to intersection
 int id = 0;
 Ray r = r_;
 int depth = depth_;

 Vec color(0,0,0); //accumulated color
 Vec reflect(1,1,1); //accumulated reflectance

 while(1) {
  if (!intersect(r, t, id)) //if missed
   return color;

  const Sphere &obj = spheres[id];  // hit object
  Vec x = r.o + r.d * t;     // intersection point
  Vec n = (x-obj.pos).norm();   //surface normal
  Vec nl = n.dot(r.d) < 0 ? n : n*(-1); //dot becomes negative if ray is goin into the sphere.
  Vec oc = obj.color;     //object color (reflectance).

  color = color + reflect.mult(obj.emi);  //accumulate color if any emission
  //printf("%f, %f, %f", color.x, color.y, color.z);

  double p = oc.x > oc.y && oc.x > oc.z ? oc.x : oc.y > oc.z ? oc.y : oc.z; // max object reflectance
  
  //russian roulet
  if (drand48() < p) {
   oc = oc * (1/p); //normalize by maximum reflectnce
  } else {
   return color;
  }

  reflect = reflect.mult(oc); //accumulate reflectance

  if (obj.refl == DIFF) { //lambert diffusion
   double r1 = 2 * M_PI * drand48();
   double r2 = drand48(), r2s = sqrt(r2);

   // new axis @ a intersection point
   Vec w = nl;
   Vec u = ((fabs(w.x) > 0.1 ? Vec(0,1,0): Vec(1,0,0))%w).norm();
   Vec v = w % u;

   //random sample in a semisphere
   Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();

   r = Ray(x, d);
   continue;


  } else if (obj.refl == SPEC) { //specular reflection (like mirror)
   r = Ray(x, r.d - n*2*n.dot(r.d));
   continue;
  }

  //dielectric Refraction
  Ray reflRay(x, r.d - n*2*n.dot(r.d)); //reflection ray
  bool into = n.dot(nl) > 0;     // true if ray is goin into the object
  double ni = 1.0, nr = 1.5;
  double nir = into ? ni/nr : nr/ni;
  double ddn = r.d.dot(nl);
  double cos2t = 1 - nir*nir*(1-ddn*ddn);

  if (cos2t < 0) { //total reflection
   r = reflRay;
   continue;
  }
  Vec refr_dir = (r.d * nir - n*((into?1:-1) * (ddn*nir + sqrt(cos2t)))).norm(); //refracted ray
  double R0 = (nr-ni)*(nr-ni)/((nr+ni)*(nr+ni));  //normal reflectance
  double co = 1 - (into ? -ddn : refr_dir.dot(n)); 

  double Re = R0 + (1 - R0) * co * co * co * co * co; //sheered reflectance
  double Tr = 1 - Re;

  if (drand48() < Re) { 
   r = reflRay;  //reflected
  } else {
   r = Ray(x, refr_dir); //refracted
  }
  continue;
 }
}

int main(int argc, char *argv[]){
 srand(time(NULL));
 clock_t sclock, eclock;

 int samps = argc==2 ? atoi(argv[1]): 1; // # samples
   int xres = 256, yres = xres;
   int ures = 1, vres = ures;
   Sensor snsr(100., 100., xres, yres);
   double fl = 100;
   Lens lens(fl, fl/1.4, fl/1.4, ures, vres);
   Camera cam(snsr, lens, 2*fl, Vec(0.0, -0.1, 50+fl*2), Vec(0., 0., -1.), Vec(0.,1.,0.));
   int uv = cam.lens->xres * cam.lens->yres;
   Vec r;

//#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP
   sclock = clock(); 
   for (int y=0; y < cam.snsr->yres; y++){                       // Loop over image rows
      fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps,100.*y/(cam.snsr->yres-1));
     for (int x=0; x < cam.snsr->xres; x++) {   // Loop over img cols
      Vec pconj = cam.conj(x,y);
      r = Vec(); //radiance

      for (int v=0; v<cam.lens->yres; v++) {
       for (int u=0; u<cam.lens->xres; u++) {
        Vec p = cam.lens_coord(u,v);
        Vec d = (pconj - p).norm();//lens to the conj point
        
        for(int s=0; s<samps; s++) {
         r = r + radiance(Ray(p,d), 0) * (1./(samps*uv));
        }
       }
      }
      int i = (cam.snsr->xres-1-y)*cam.snsr->xres + x;
           cam.snsr->c[i] = Vec(clamp(r.x),clamp(r.y),clamp(r.z));
        }
   }
   eclock = clock();
   printf("time: %.2f1[s]\n",(double)(eclock-sclock)/CLOCKS_PER_SEC);

   FILE *f = fopen("image.ppm", "w");         // Write image to PPM file.
   fprintf(f, "P3\n%d %d\n%d\n", cam.snsr->xres, cam.snsr->yres, 255);
   for (int i=0; i< cam.snsr->xres * cam.snsr->yres; i++) {
     fprintf(f,"%d %d %d ", toInt(cam.snsr->c[i].x), toInt(cam.snsr->c[i].y), toInt(cam.snsr->c[i].z));
    }
}