### Extract Camra Response from Random Image Data¶

The Problem`

Applications exist where is is useful to know the spectral response of an image sensor in commercial cameras. One approach to determine this is present many different images (i.e. Power Spectral Density, PSD) to the sensor, and record the R,G, and B sensors output.

Imports ..

Functions

Read No Noise Full File PSD Data

First Let's Plot the Alphas (The 'Solution')

### PCA¶

To reduce dimensionality PCA was run on the 301 feature dataset. The first 20 principle components will computed to reduce the fitness search space.

Some PCA Explanations

PCA Detail

PCA Methods

PCA Reconstruction

PCA for Dummies

#### Now Let's Try the Real Data¶

A prior PCA was computed in Haskell / R.

Here we read that in for eventual comparison to the PCA computed by SkLearn

#### Comparing R results and Python Results .. Signs are sometimes different!¶

It's Ok ..

Let's see how well behaved the first 20 PCAs are ...

This is important since this is what we are assuming is somewhat smooth. Here we plot three samples from the post PCA dataset. This is the equivalent of printing the PSD for all 301 points, but we limit it to the 20 PCA cols.

#### We're almost ready to send it off to Haskell for learning ..¶

To prepare for scaling we have to stitch together the PCAs and the 'Answer"

#### Now we can scale!¶

Write the scaled data out so Haskell can do fitness ..

#### But first let's verify the fitness code is working¶

Individual ... [0.325880913228656,0.39667643306649025,0.7473014477142477,8.434155613051841e-2,-2.6434725189356578e-2,-1.4737156547028556e-2]

Fitness Cases [FitnessCase [-0.23213254724324522,-0.6107708369235977,0.23999887959073862,0.8467276982575322,0.48775006287997447,-2.5080243207554913,-1.6827055558196122,-0.8305324734160002,-1.3677144725147194,0.8701923515463181,-3.4105364359106605,1.8618782853379916,-3.554742358123332,1.161281390129975,0.259121070757041,-1.494009766131712,-0.8479441055976709,1.5123999298371094,-0.7778076500051748,-0.2892190888978542] (-9.891512311482657e-2),FitnessCase [-0.6231630714884775,1.2510387625530222,-0.9685913744980068,-2.134637737482353,-1.1802306213456357,-0.3687864365511998,0.6374883143095859,-0.10991819345725727,0.3213625883996374,0.6128359489871866,-0.6762067819575514,1.7359594408021162,1.3894764274474027,-1.0282933311986495,1.5158998150040963,-2.258650865948252,-3.727967020197284,-5.486220238485005,1.7049540836202024,9.340896540033592] (-0.9926610004408698)]

#### Compute the polynomial¶

Note that the length of the fitness case inputs is 20, and the range of x's is
$\ \ \ \ \ 0.0 < x < 3.0$

#### Now let's see if we can get a reasonable polynomial shape ..¶

Guessing first

After playing with it a bit it looks like the search space for a reasonable shape (not just a monotinically increasing/decreasing curve) is big .. Let's try to fix this by restricting the $a's$ using some constraints ..

#### Try some parametric calculations ..¶

Try to compute $a_0 \rightarrow a_4$ by giving two constraints, \begin{align} y = 0\ \ at\ \ x = -5\ \ and\ \ x = 5\\ \end{align} and \begin{align} \frac{dy}{dx} = 0\ \ at\ x = -5 \ \ and\ \ x = 5 \end{align}

So,

\begin{align} &y = a_0 + a_1x + a_2x^2 + a_3x^3 + a4x^4 \label{y}\\ &\frac{dy}{dx} = a_1 + 2a_2x + 3a_3x^2 + 4a_4x^3\label{dydx} \end{align}

Take Eq $\eqref{y}$ and substitute 5, and -5 for x,

\begin{align} 0 &= a_0 + 5a_1 + 25a_2 + 125a_3 + 625 a_4\label{yplus}\ \ \ \ \ (x=5) \\ 0 &= a_0 - 5a_1 + 25a_2 - 125a_3 + 625 a_4\label{yminus}\ \ \ \ \ (x=-5) \end{align}

Adding Eq's $\eqref{yplus}$ and $\eqref{yminus}$, and also subtracting Eq $\eqref{yminus}$ from Eq $\eqref{yplus}$ yields,

\begin{align} 0 &= 2a_0 + 50a_2 + 1250a_4 \ \ \ \ \ (adding\ Eq\ \eqref{yplus} \ and\ Eq\ \eqref{yminus}) \\ 0 &= 10a_1 + 250a_3 \ \ \ \ \ \ \ \ \ \ \ (subtracting\ Eq\ \eqref{yminus}\ from\ Eq\ \eqref{yplus}) \end{align}

Let's do the same thing with the derivative equation $\eqref{dydx}$

\begin{align} 0 &= a_1 + 10a_2 + 75a_3 + 500a_4\label{dydxplus}\ \ \ \ \ (x=5) \\ 0 &= a_1 - 10a_2 + 75a_3 - 500a_4\label{dydxminus}\ \ \ \ \ (x=-5) \end{align}

Adding Eq's $\eqref{dydxplus}$ and $\eqref{dydxminus}$, and also subtracting Eq $\eqref{dydxminus}$ from Eq \eqref{dydxplus} yields, \begin{align} 0 &= 2a_1 + 150a_3 \ \ \ \ \ \ \ \ \ \ (adding\ Eq\ 3 \ and\ Eq\ 4) \\ 0 &= 20a_2 + 1000a_4 \ \ \ \ \ \ (subtracting\ Eq\ 4\ from\ Eq\ 3) \end{align} Collecting the parametric equations .. \begin{align} 0& = 2a_0& + &\ \ 50a_2 + &&1250a_4 \\ 0& = &10a_1\ \ + &&250a_3\label{a1a3y}\\ \\ 0& = &2a_1\ \ + &&150a_3\label{a1a3dydx} \\ 0& = &&\ \ 20a_2 + &&1000a_4 \end{align} This is not completely solvable but we can do a little work .. But wait .... There is a conflict if we look at Eq's\eqref{a1a3y}$and$\eqref{a1a3dydx}$!! $$a_1 = -25a_3\\ a_1 = -75a_3$$ Let's try both and plot them ... #### Note that one curve obeys the derivative = 0, and the other obeys the value = 0, both at x = +5, and x=-5¶ ### New Ideas¶ #### Ok .. New Ideas .. Use FFTs or a polynomial instead of alphas for the 'solution'.¶ • FFT's .. convert the PSD into an FFT sequence, evolve an FFT 'multiplier' .. Not sure if it should be a multiply or an addition. I think multiplying FFTs is the equivalent of convolution. • Polynomial .. See below .. #### Polynomial First¶ For example, $$P_{w1}*f(1) + P_{w2}*f(2) + P_{w3}*f(3) + .... P_{w301}*f(301) = S_1 \ \ \ \ \ \textbf{Eq (2)}$$ $$\vdots$$ Where f is a polynomial, $$f(x) = a_0 + a_1*x + a_2*x^2 + a_3*x^3 ... a_n*x^n\ \ \ \ \ \textbf{Eq (3)}$$ And P is the PSD input vector. And R is the scaler red sensor response Ok, Here we will just take the 20 cols of PCA with unscaled output .. and try the polynomial idea .. Later we can try scaling the output with zero centered mean and a variance of one. So we let Haskell and a basic EC algorithm work on a solution. The representation was 6 floating point numbers,$a_0 \rightarrow a_5$as shown in Eq (3). During the evolution a candate solution is chosen, then each of the 301 points is passed through the polynomial, resulting in 301 new numbers. A dot product is now performed with the input data as shown in Eq (2) yielding$S_x$. This is a potential solution which is compared with the actual solution to generate the actual fitness used for subsequent selection. For the first test three PCA columns were used as fitness cases, with a polynomial length of 4. It trained well as the original fitness was$10^{13}\$ and ended up after about 1700 generations at 250.

It's best solution was,

Let's plot it .. Note: Python Lambda, Map, and Filter

First let's plot the real solution ..

#### Now try the FFT¶

Now the inverse FFT

How can I find the point spread function ( impulse response function) of a camera? Using a 2D-Convolution technique

Maybe you could guess and use deconvlucy() until you get your fuzzballs as point-like as possible.