So I have an array (it's large - 2048x2048), and I would like to do some element wise operations dependent on where they are. I'm very confused how to do this (I was told not to use for loops, and when I tried that my IDE froze and it was going really slow).

Onto the question:

h = aperatureimage
h[:,:] = 0
indices = np.where(aperatureimage>1)

for True in h:
    h[index] = np.exp(1j*k*z)*np.exp(1j*k*(x**2+y**2)/(2*z))/(1j*wave*z)

So I have an index, which is (I'm assuming here) essentially a 'cropped' version of my larger aperatureimage array. *Note: Aperature image is a grayscale image converted to an array, it has a shape or text on it, and I would like to find all the 'white' regions of the aperature and perform my operation.

How can I access the individual x/y values of index which will allow me to perform my exponential operation? When I try index[:,None], leads to the program spitting out 'ValueError: broadcast dimensions too large'. I also get array is not broadcastable to correct shape. Any help would be appreciated!

One more clarification: x and y are the only values I would like to change (essentially the points in my array where there is white, z, k, and whatever else are defined previously).

EDIT:

I'm not sure the code I posted above is correct, it returns two empty arrays. When I do this though index = (aperatureimage==1) print len(index)

Actually, nothing I've done so far works correctly. I have a 2048x2048 image with a 128x128 white square in the middle of it. I would like to convert this image to an array, look through all the values and determine the index values (x,y) where the array is not black (I only have white/black, bilevel image didn't work for me). I would then like to take all the values (x,y) where the array is not 0, and multiply them by the h[index] value listed above.

I can post more information if necessary. If you can't tell, I'm stuck.

EDIT2: Here's some code that might help - I think I have the problem above solved (I can now access members of the array and perform operations on them). But - for some reason the Fx values in my for loop never increase, it loops Fy forever....

import sys, os
from scipy.signal import *
import numpy as np
import Image, ImageDraw, ImageFont, ImageOps, ImageEnhance, ImageColor



def createImage(aperature, type):
    imsize = aperature*8
    middle = imsize/2
    im = Image.new("L", (imsize,imsize))
    draw = ImageDraw.Draw(im)
    box = ((middle-aperature/2, middle-aperature/2), (middle+aperature/2, middle+aperature/2))
import sys, os
from scipy.signal import *
import numpy as np
import Image, ImageDraw, ImageFont, ImageOps, ImageEnhance, ImageColor



def createImage(aperature, type):

    imsize = aperature*8 #Add 0 padding to make it nice
    middle = imsize/2 # The middle (physical 0) of our image will be the imagesize/2
    im = Image.new("L", (imsize,imsize)) #Make a grayscale image with imsize*imsize pixels
    draw = ImageDraw.Draw(im) #Create a new draw method
    box = ((middle-aperature/2, middle-aperature/2), (middle+aperature/2, middle+aperature/2)) #Bounding box for aperature
    if type == 'Rectangle':
        draw.rectangle(box, fill = 'white') #Draw rectangle in the box and color it white
    del draw
    return im, middle


def Diffraction(aperaturediameter = 1, type = 'Rectangle', z = 2000000, wave = .001):

    # Constants
    deltaF = 1/8 # Image will be 8mm wide
    z = 1/3.
    wave = 0.001
    k = 2*pi/wave

    # Now let's get to work
    aperature = aperaturediameter * 128 # Aperaturediameter (in mm) to some pixels
    im, middle = createImage(aperature, type) #Create an image depending on type of aperature 
    aperaturearray = np.array(im) # Turn image into numpy array

    # Fourier Transform of Aperature
    Ta = np.fft.fftshift(np.fft.fft2(aperaturearray))/(len(aperaturearray)) 

    # Transforming and calculating of Transfer Function Method
    H = aperaturearray.copy() # Copy image so H (transfer function) has the same dimensions as aperaturearray
    H[:,:] = 0 # Set H to 0
    U = aperaturearray.copy()
    U[:,:] = 0
    index = np.nonzero(aperaturearray) # Find nonzero elements of aperaturearray
    H[index[0],index[1]] = np.exp(1j*k*z)*np.exp(-1j*k*wave*z*((index[0]开发者_开发百科-middle)**2+(index[1]-middle)**2)) # Free space transfer for ap array
    Utfm = abs(np.fft.fftshift(np.fft.ifft2(Ta*H))) # Compute intensity at distance z

    # Fourier Integral Method
    apindex = np.nonzero(aperaturearray)
    U[index[0],index[1]] = aperaturearray[index[0],index[1]] * np.exp(1j*k*((index[0]-middle)**2+(index[1]-middle)**2)/(2*z))
    Ufim = abs(np.fft.fftshift(np.fft.fft2(U))/len(U))

    # Save image
    fim = Image.fromarray(np.uint8(Ufim))
    fim.save("PATH\Fim.jpg")
    ftfm = Image.fromarray(np.uint8(Utfm))
    ftfm.save("PATH\FTFM.jpg")
    print "that may have worked..."
    return

if __name__ == '__main__':
    Diffraction()

You'll need numpy, scipy, and PIL to work with this code.

When I run this, it goes through the code, but there is no data in them (everything is black). Now I have a real problem here as I don't entirely understand the math I'm doing (this is for HW), and I don't have a firm grasp on Python.

U[index[0],index[1]] = aperaturearray[index[0],index[1]] * np.exp(1j*k*((index[0]-middle)**2+(index[1]-middle)**2)/(2*z))

Should that line work for performing elementwise calculations on my array?

Could you perhaps post a minimal, yet complete, example? One that we can copy/paste and run ourselves?

In the meantime, in the first two lines of your current example:

h = aperatureimage
h[:,:] = 0

you set both 'aperatureimage' and 'h' to 0. That's probably not what you intended. You might want to consider:

h = aperatureimage.copy()

This generates a copy of aperatureimage while your code simply points h to the same array as aperatureimage. So changing one changes the other. Be aware, copying very large arrays might cost you more memory then you would prefer.

What I think you are trying to do is this:

import numpy as np

N = 2048
M = 64

a = np.zeros((N, N))
a[N/2-M:N/2+M,N/2-M:N/2+M]=1

x,y = np.meshgrid(np.linspace(0, 1, N), np.linspace(0, 1, N))

b = a.copy()

indices = np.where(a>0)

b[indices] = np.exp(x[indices]**2+y[indices]**2)

Or something similar. This, in any case, sets some values in 'b' based on the x/y coordinates where 'a' is bigger than 0. Try visualizing it with imshow. Good luck!

Concerning the edit

You should normalize your output so it fits in the 8 bit integer. Currently, one of your arrays has a maximum value much larger than 255 and one has a maximum much smaller. Try this instead:

fim = Image.fromarray(np.uint8(255*Ufim/np.amax(Ufim)))
fim.save("PATH\Fim.jpg")
ftfm = Image.fromarray(np.uint8(255*Utfm/np.amax(Utfm)))
ftfm.save("PATH\FTFM.jpg")

Also consider np.zeros_like() instead of copying and clearing H and U.

Finally, I personally very much like working with ipython when developing something like this. If you put the code in your Diffraction function in the top level of your script (in place of 'if __ name __ &c.'), then you can access the variables directly from ipython. A quick command like np.amax(Utfm) would show you that there are indeed values!=0. imshow() is always nice to look at matrices.