Analisi lunari mediante strumenti informatici

Di più
3 Anni 10 Mesi fa - 3 Anni 10 Mesi fa #39722 da doktorenko
Questa è una discussione sperimentale: nel senso che i partecipanti sono invitati a ripetere e a riprodurre le esperienze proposte; tutto il codice necessario per generare le immagini verrà pubblicato.

Il principio per la formazione delle immagini è quello della camera ottica, ma capovolto per ragioni di efficienza: invece di intercettare un determinato fascio di raggi luminosi riflessi dal soggetto, si inviano dal sensore verso il soggetto tanti raggi quanti sono i punti sensibili della nostra pellicola virtuale; raggi quindi che hanno origine ognuno in un punto distinto del sensore, che passano tutti per uno stesso punto (il cosiddetto foro stenopeico), e che da questo foro di separano per incontrare -eventualmente- la superficie di un soggetto illuminato.

La prima parte del codice è relativa appunto alla formazione di questi raggi; come linugaggio di programmazione uso, per iniziare, python e il modulo numpy.

Creiamo quindi la nostra camera ottica (la convezione 3d usata è: asse y altezza, z profondità) :
Code:
import numpy as np prodottoS=lambda a,b: np.sum(a*b,axis=0) normalizza=lambda a: a/np.linalg.norm(a,axis=0) class Camera: def __init__(self,centro,risoluzione): self.c=np.array(centro).reshape(3,1) F=.61 mm=.55 xv,yv=np.meshgrid(np.linspace(-.5,.5,risoluzione), np.linspace(-.5,.5,risoluzione)) zv=np.zeros((risoluzione,risoluzione))-F o=np.stack((xv*mm,yv*mm,zv)) self.r=-normalizza(o).reshape(3,-1) self.o=(o+self.c[:,np.newaxis]).reshape(3,-1) risoluzione=128 cam=Camera([0,.35,-3],risoluzione)

La camera "cam" così costruita:

1) ha un sensore quadrato formato da 128*128 punti ("risoluzione")
2) ha centro (inteso come il punto in comune del fascio di raggi proiettato) nel punto (0,.35,-3)
3) ha una focale fissa di .61 e una dimensione del sensore di .55 per lato
4) non si può orientare (punta direttamente l'asse z)

l'origine di ciascun raggio è memorizzata in "cam.o", la direzione in "cam.r"

fine della prima parte
Ultima Modifica 3 Anni 10 Mesi fa da doktorenko. Motivo: modifica codice per renderlo piu` aderente alla descrizione

Si prega Accesso a partecipare alla conversazione.

Di più
3 Anni 10 Mesi fa - 3 Anni 10 Mesi fa #39761 da doktorenko
Nella seconda parte, il codice che si occupa di calcolare le intersezioni tra le forme geometriche (per il momento solamente sfere) e i raggi che partono dalla camera:
Code:
class Sfera: C=0 R=[] n=1 def __init__(self,centro,raggio,colore): if (self.n==1): Sfera.C=np.array(centro).reshape(3,1) self.c=np.array(centro).reshape(3,1) self.r=raggio Sfera.C=np.concatenate((Sfera.C,self.c)) Sfera.R.append(raggio) self.r2=raggio*raggio self.col=np.array(colore) self.i=Sfera.n Sfera.n+=1 def intercetta(self,origine,direzione,d0,ind): L=self.c-origine np.seterr(all='ignore') tca=prodottoS(L,direzione) d2= prodottoS(L,L)-tca*tca thc=np.sqrt(self.r2-d2) t0=tca-thc t1=tca+thc t0=np.where(t1<t0,t1,t0) t0=np.where(t0<0,np.nan,t0) m=np.where(t0<d0) d0[m]=t0[m] ind[m]=self.i np.seterr(all='warn')

Per Il calcolo dell`intersezione tra raggio e un punto della sfera si procede nella maniera seguente:



1) si traccia L (segmento che unisce il punto di partenza O del raggio al centro C della sfera)
2) si calcola la lunghezza di tca (segmento OB) come prodotto scalare di L per il versore direzione raggio r
3) si ricava la lunghezza di d (cateto BC del triangolo rettangolo OBC)
4) si ricava la lunghezza thc (cateto del triangolo rettangolo CBP, con CP raggio della sfera)
5) si ricavano le due lunghezze t0 e t1 (distanza OP0 e OP1) rispettivamente sottraendo e sommando tca e thc

(il calcolo dei punti tridimensionali P0 e P1 viene eseguito in un secondo momento)

fine seconda parte (qui i dettagli della procedura)
Ultima Modifica 3 Anni 10 Mesi fa da doktorenko.

Si prega Accesso a partecipare alla conversazione.

Di più
3 Anni 10 Mesi fa #39774 da doktorenko
Terza parte.

Aggiungiamo una funzione che restituisce la distanza tra l`origine di ciascun raggio e l`eventuale oggetto incontrato (contestualmente viene restituito anche l`indice dello stesso oggetto). Se il raggio non incontra ostacoli nel suo cammino, il risultato convenzionale restituito dalla funzione è; indice=0, distanza=1e20.
Code:
def distanze(o,r,scena): t=np.zeros(o.shape[1:])+1e20 ind=np.zeros(t.shape).astype(int) for i in scena: i.intercetta(o,r,t,ind) return t,ind

Come prima prova aggiungiamo una funzione semplificata di disegno della scena:
Code:
def disegna(o_camera,r_camera,scena,luce): t,ind=distanze(o_camera,r_camera,scena) img=np.zeros((t.shape)) visibile=(ind>0) img[visibile]=ind[visibile] return img

Il codice principale, con la creazione della scena (4 sfere), e la visuzalizzazione dell'immagine prodotta (in questo caso uso il modulo matplotlib, ma si possono usare anche altri moduli):
Code:
luce=normalizza(np.array([0, 5, 10]).reshape(3,1)) scena=[Sfera([.75,.1,1.],.6,[1.0,1.0,1.0]), Sfera([-.75,.1,2.25],.6,[1.0,1.0,1.0]), Sfera([-2.75,.1,3.5],.6,[1.0,1.0,1.0]), Sfera([0,-99999.5,0],99999.,[1.0,1.0,1.0])] Sfera.C=Sfera.C.reshape(-1,3) Sfera.R=np.array(Sfera.R) img=disegna(cam.o,cam.r,scena,luce) from matplotlib import pyplot as plt plt.imshow(img.reshape((risoluzione,risoluzione))) plt.show()

risultato dell'elaborazione;



fine terza parte

Si prega Accesso a partecipare alla conversazione.

Di più
3 Anni 10 Mesi fa - 3 Anni 10 Mesi fa #39783 da doktorenko
Nella quarta parte scriviamo la funzione che calcola la luminosità del suolo lunare (il cuore del programma). Per semplicità, usiamo i valori precalcolati da Scevcenco; eccoli in forma di immagine 16 bit, formato png:



(se vedete un'immagine completamente nera non preoccupatevi)

La codice è questo:
Code:
def scev(N,L,R): NL=prodottoS(N,L) NR=prodottoS(N,-R) phi_l=normalizza(L-NL*N) phi_r=normalizza(-R-NR*N) cos_phi=prodottoS(phi_l,phi_r) phi=(np.degrees(np.arccos(cos_phi))/10).astype(int) thr=(np.degrees(np.arccos(NR))/10).astype(int) thl=(np.degrees(np.arccos(NL))/10).astype(int) lum=TAB[thl,thr,phi] return lum

Lo schema di funzionamento:



in entrata la funzione richiede 3 versori:

N) normale al piano p colpito dalla luce nel punto P
L) versore con origine P che punta la fonte di luce
R) versore con origine P che punta verso una determinata cella del sensore

vengono in seguito calcolati:

angolo phi := angolo tra le proiezioni di L e R sul piano p
angolo th_r := angolo tra R e N
angolo th_l := angolo tra L e N

e finalmente si recuperano i valori di luminosità per questi stessi angoli dalla tabella elaborata da Sc.

Prima di terminare, controlliamo se la tabella è stata ricavata con successo dall'immagine png allegata in precedenza; scriviamo questo piccolo programma di verifica:
Code:
from matplotlib import pyplot as plt b=plt.imread("scevcenco.png")*65.535 plt.imshow(b) plt.show()

risultato:

[img


fine quarta parte (nella quinta parte useremo la funzione su oggetti tridimensionali, per simulare il suolo lunare)
Ultima Modifica 3 Anni 10 Mesi fa da doktorenko. Motivo: codice più compatto; corrette didascalie

Si prega Accesso a partecipare alla conversazione.

Di più
3 Anni 10 Mesi fa - 3 Anni 10 Mesi fa #39810 da doktorenko
Nella quinta parte modifichiamo la funzione che disegna gli oggetti in modo da aggiungere ombre e riflessi; innanzitutto inseriamo una funzione per il calcolo delle coordinate dei punti raggiunti dai raggi uscenti dalla camera:
Code:
def disegna(O,R,scena,l): Pxyz = lambda i: O[:,i]+t[i]*R[:,i]# xyz_punto = xyz_origine_raggio + distanza * direzione_raggio

poi una funzione che calcola la riflessione prodotta da una superficie lucidata a specchio:
Code:
riflesso = lambda i : normalizza(R[:,i]-N[:,i]*2*prodottoS(R[:,i],N[:,i]))

raggio_uscente = normalizza(raggio_incidente - normale_superficie * 2 (raggio_incidente x normale_superficie) )

due pseudo funzioni per sintetizzare il codice:
Code:
distanza = lambda i : distanze(O[:,i]+N[:,i]*0.00001,R[:,i],scena) normale = lambda i : normalizza(O[:,i]-Sfera.C[ind[i]].T)

Calcoliamo ora le distanze raggi_camera - oggetti visibili:
Code:
img=np.zeros((R.shape[1:])) N=np.zeros((R.shape)) t,ind=distanze(O,R,scena)

creiamo tre indici:

1) oggetti := tutti i raggi che colpiscono un oggetto nell`inquadratura
2) regolite := sfera #4, materiale Scev. (sottoinsieme di 1)
3) riflettenti := altre sfere (sottoinsieme di 1)
Code:
oggetti=np.where(ind!=0)[0] regolite=oggetti[(ind[oggetti]==4)] riflettenti=oggetti[ind[oggetti]!=4]

calcoliamo ora le coordinate 3d di ogni punto raggiunto (indice 1), e la normale alla superficie in quello stesso punto:
Code:
O[:,oggetti]=Pxyz(oggetti) N[:,oggetti]=normalizza(O[:,oggetti]-Sfera.C[ind[oggetti]].T)

cambiamo la direzione ai raggi che raggiungono superfici riflettenti (indice #3):
Code:
R[:,riflettenti]=riflesso(riflettenti)

calcoliamo la distanza tra l'origine di questi nuovi raggi (punto di rimbalzo) e le eventuali superfici toccate:
Code:
t[riflettenti],ind[riflettenti]=distanza(riflettenti)

di questi nuovi raggi, raggruppiamo in un nuovo indice tutti quelli che toccano la superficie lunare (indice #4 := regolite_riflessa):
Code:
regolite_riflessa=riflettenti[ind[riflettenti]==4]

calcoliamo le coordinate (e le normali) dei punti toccati:
Code:
O[:,regolite_riflessa]=Pxyz(regolite_riflessa) N[:,regolite_riflessa]=normale(regolite_riflessa)

raggruppiamo gli indici #4 e #2 in un unico indice (tutta la regolite inquadrata direttamente e di riflesso):
Code:
regolite=np.concatenate((regolite,regolite_riflessa))

da questo ultimo gruppo di punti, facciamo partire dei raggi in direzione della fonte di luce:
Code:
t[regolite],ind[regolite]=distanze(O[:,regolite]+N[:,regolite]*0.00001,l,scena)

eliminiamo dal gruppo tutti i punti in ombra (sono i punti dai quali partono i raggi che incontrano un ostacolo)
Code:
regolite=regolite[ind[regolite]==0]

calcoliamo finalmente la luminosità dei punti rimanenti secondo la nostra tabella precalcolata:
Code:
img[regolite]=scev(N[:,regolite],l,R[:,regolite]) return img

riscriviamo il codice principale:
Code:
from matplotlib import pyplot as plt risoluzione=512*2 cam=Camera([0,.35,-3],risoluzione) luce=normalizza(np.array([0, 5, 10]).reshape(3,1)) scena=[Sfera([.75,.1,1.],.6,[1.0,1.0,1.0]), Sfera([-.75,.1,2.25],.6,[1.0,1.0,1.0]), Sfera([-2.75,.1,3.5],.6,[1.0,1.0,1.0]), Sfera([0,-99999.5,0],99999.,[1.0,1.0,1.0])] Sfera.C=Sfera.C.reshape(-1,3) Sfera.R=np.array(Sfera.R) TAB=(plt.imread("scevcenco.png")*65.535).reshape((9,10,19)) img=disegna(cam.o,cam.r,scena,luce) plt.imshow(img.reshape((risoluzione,risoluzione))) plt.show()

risultato (in falsi colori):



fine quinta parte (nel prossimo messagiio invio il listato unito)
Ultima Modifica 3 Anni 10 Mesi fa da doktorenko. Motivo: correzione codice

Si prega Accesso a partecipare alla conversazione.

Di più
3 Anni 10 Mesi fa #39871 da doktorenko
Codice riunito (e ottimizzato):
Code:
import numpy as np def prodottoS(a,b): if (b.ndim==1): return a.dot(b) else: return np.einsum("ij,ij->i",a,b) def normalizza(a): if (a.ndim==1): return a/np.linalg.norm(a,axis=0) else: return a/np.expand_dims(np.sqrt(prodottoS(a,a)),axis=1) class Camera: def __init__(self,centro,risoluzione): global o self.c=np.array(centro) F=.61 mm=.55 xv,yv=np.meshgrid(np.linspace(-.5,.5,risoluzione), np.linspace(-.5,.5,risoluzione)) zv=np.zeros((risoluzione,risoluzione))-F o=np.dstack((xv*mm,yv*mm,zv)).reshape(-1,3) self.r=-normalizza(o) self.o=(o+self.c) class Sfera: C=0 R=[] n=1 def __init__(self,centro,raggio,colore): if (self.n==1): Sfera.C=np.array(centro) self.c=np.array(centro) self.r=raggio Sfera.C=np.concatenate((Sfera.C,self.c)) Sfera.R.append(raggio) self.r2=raggio*raggio self.col=np.array(colore) self.i=Sfera.n Sfera.n+=1 def intercetta(self,origine,direzione,d0,ind): L=self.c-origine np.seterr(all='ignore') tca=prodottoS(L,direzione) d2= prodottoS(L,L)-tca*tca thc=np.sqrt(self.r2-d2) t0=tca-thc t1=tca+thc t0=np.where(t1<t0,t1,t0) t0=np.where(t0<0,np.nan,t0) m=np.where(t0<d0) d0[m]=t0[m] ind[m]=self.i np.seterr(all='warn') def distanze(o,r,scena): t=np.zeros(o.shape[:-1])+1e20 ind=np.zeros(t.shape).astype(int) for i in scena: i.intercetta(o,r,t,ind) return t,ind def disegna(O,R,scena,l): Pxyz = lambda i: O[i]+R[i]*t[i,np.newaxis] riflesso = lambda i : normalizza( R[i]-(N[i]*2*prodottoS(R[i],N[i])[:,np.newaxis])) distanza = lambda i : distanze(O[i]+N[i]*0.00001,R[i],scena) normale = lambda i : normalizza(O[i]-Sfera.C[ind[i]]) img=np.zeros((R.shape[:-1])) N=np.zeros((R.shape)) t,ind=distanze(O,R,scena) oggetti=np.where(ind!=0)[0] regolite=oggetti[(ind[oggetti]==4)] riflettenti=oggetti[ind[oggetti]!=4] O[oggetti]=Pxyz(oggetti) N[oggetti]=normalizza(O[oggetti]-Sfera.C[ind[oggetti]]) R[riflettenti]=riflesso(riflettenti) t[riflettenti],ind[riflettenti]=distanza(riflettenti) regolite_riflessa=riflettenti[ind[riflettenti]==4] O[regolite_riflessa]=Pxyz(regolite_riflessa) N[regolite_riflessa]=normale(regolite_riflessa) regolite=np.concatenate((regolite,regolite_riflessa)) t[regolite],ind[regolite]=distanze(O[regolite]+N[regolite]*0.00001,l,scena) regolite=regolite[ind[regolite]==0] img[regolite]=scev(N[regolite],l,R[regolite]) return img def scev(N,L,R): NL=prodottoS(N,L) NR=prodottoS(N,-R) phi_l=normalizza(L-N*NL[:,np.newaxis]) phi_r=normalizza(-R-N*NR[:,np.newaxis]) cos_phi=prodottoS(phi_l,phi_r) phi=(np.degrees(np.arccos(cos_phi))/10).astype(int) thr=(np.degrees(np.arccos(NR))/10).astype(int) thl=(np.degrees(np.arccos(NL))/10).astype(int) return TAB[thl,thr,phi] from matplotlib import pyplot as plt import time risoluzione=512*2 cam=Camera([0,.35,-3],risoluzione) luce=np.array([0, 5, 10]) luce=normalizza(luce) scena=[Sfera([.75,.1,1.],.6,[1.0,1.0,1.0]), Sfera([-.75,.1,2.25],.6,[1.0,1.0,1.0]), Sfera([-2.75,.1,3.5],.6,[1.0,1.0,1.0]), Sfera([0,-99999.5,0],99999.,[1.0,1.0,1.0])] Sfera.C=Sfera.C.reshape(-1,3) Sfera.R=np.array(Sfera.R) TAB=(plt.imread("scevcenco.png")*65.535).reshape((9,10,19)) t0=time.time() img=disegna(cam.o,cam.r,scena,luce) print (time.time()-t0) plt.imshow(img.reshape((risoluzione,risoluzione))) plt.show()

Si prega Accesso a partecipare alla conversazione.

Tempo creazione pagina: 0.359 secondi
Powered by Forum Kunena