0

I have a 1D function of N values. I want to make a 2D histogram so that I get an image of nx*ny pixels and then cut the image summing along 1 dimension. The function before and after should be the same. I tried with a gaussian but I am missing a factor sqrt. Please see the code. Am I missing something in my fucntion?

I draw r from a random number

import numpy as np 
import random
import matplotlib.pyplot as plt
sigma=0.5
N=100000

r=np.random.normal(0, sigma, N)
ind = np.where( r>=0 )
r   =r[ind]
N=len(r)
phi=2*np.pi*np.random.rand(N)
x=r*np.cos(phi)
y=r*np.sin(phi)
Ir=np.zeros(N)
Ir[:]=1

Now I want the see the distribution Ir=f(x, y), as a 2D image.

def getImage2D(x, y, fun, nx=100, ny=100, xmin=-1, xmax=1, ymin=-1, ymax=1):

dx=(xmax-xmin)*(1.0/nx);
dy=(ymax-ymin)*(1.0/ny);
image  = np.ndarray(shape=(nx, ny), dtype=float); image.fill(0.0)

for i in range(len(fun)):
    mr = (x[i]-xmin)/dx;
    nr = (y[i]-ymin)/dy;
    m, n=int(mr), int(nr)
    image[m, n]=image[m, n]+fun[i];

return image


P=getImage2D(x, y, Ir, nx=101, ny=101, xmin=-3, xmax=3, ymin=-3, ymax=3)
#P=getImage2D(x, y, r**0.5*Ir, nx=101, ny=101, xmin=-3, xmax=3, ymin=-3, ymax=3)
Px=np.sum(P, axis=1)
Px=Px/np.max(Px)
plt.figure()
plt.imshow(P)
plt.show(block=False)

If I plot the cut Px (after summing along y), I am not getting the same gaussian of width sigma!!

v=np.linspace(np.min(r), np.max(r), len(r))
v=v/np.max(v)
plt.figure()
plt.plot(np.linspace(-3, 3, 101), Px)
plt.plot(np.sort(r)[::-1], v, 'g')
plt.show(block=False)

Why the width is not same in both cases? If I put a weight r**0.5 with Ir then the width is same(sigma=0.5). Is there any mistake that I am doing in getImage2D function?

SUV
  • 53
  • 6
  • Can't you just use 2d `histogram` from `numpy`? – Sleepyhead Apr 16 '14 at 00:02
  • nuppy histogram give the same result as my getImage2D function but here also I need a factor r**0.5 in weights but why it is?? VV=np.histogram2d(x, y, bins=101, weights=Ir) #VV=np.histogram2d(x, y, bins=101, weights=r**0.5*Ir) H, xedges, yedges =VV H=np.sum(H, axis=1) plt.figure() plt.plot(np.linspace(xedges[0], xedges[-1], 101), H/np.max(H)) plt.plot(np.sort(r)[::-1], v, 'r') plt.show(block=False) – SUV Apr 16 '14 at 06:40

0 Answers0