lunedì 30 maggio 2016

Foto aeree con Google Carboard

Visto che l'esperimento con gli anaglifi era riuscito ma non nel modo ottimale ho provato a migliorare, questa volta usando un vero stereoscopio,ovvero Google Cardboard, anche se dal costo decisamente economico (uno stereoscopio da tavolo per foto aeree costa sopra il migliaio di euro)


Sul telefono ho montato una semplice pagina web sul telefono Android che scalava e centrava i due fotogrammi precedentemente registrati a mano (in versioni successive sarebbe comodo avere la possibilita' di spostare uno dei fotogrammi via telecomando bluetooh per centrare le immagini guardando nel visore)

Con questo sistema la visualizzazione 3D e' decisamente buona e molto piu' simile a quella del vero stereovisore per foto aeree rispetto a quella anaglifica

index.html
-------------------------
<head>
<link rel="stylesheet" type="text/css" href="style.css">
</head>
<body>
<div class ="stereo">
    <div id="sinistra">
      <img alt="Sinistra" title="Sinitra" src="594_2.png" style="width: 100%; max-height: 100%"/>
   </div>
    <div id="destra">
      <img alt="Sinistra" title="Sinitra" src="595_2.png" style="width: 100%; max-height: 100%"/>

   </div>
</div>
</body>
-------------------------

style.css
-------------------------
.stereo{
        width: 100%;
        overflow: hidden;
    }

    #sinistra {
        float: left;
        width: 50%;
    }
   #sinistra img
   {

      margin: auto;
      display: block;
   }
    #destra {
        float: left;
        width: 50%;
    }
   #destra img
   {
      margin: auto;
      display: block;
   }
-------------------------



Foto aeree con anaglifi

A lavoro mi capita spesso la necessita' di visualizzare foto aeree ma non ho un visore stereoscopico. Ho scoperto un sistema semplice (e lecito) per scaricare le stereocoppie dal sito della Regione Toscana (evitando un generoso esborso) e cosi' mi sono messo alla ricerca di un sistema per visualizzare la terza dimensione


Il primo tentativo e' stato quello di utilizzare gli anaglifi mediante l'uso degli occhialini rossi/ciano (qualcuno si ricorda del film lo Squalo 3D??)

Il sistema piu' semplice, una volta scelti i due fotogrammi, e' quello di usare il programma AnaBuilder che permette di coregistrare in automatico le due immagini sia come traslazione che come rotazione dal menu Actions/AutoFit with borders detection


Il risultato non e' male, anche considerando che non c'e' intervento umano, ma il programma riesce a lavorare solo su piccole porzioni della stampa. Introducendo una stampa completa della foto aerea, normalmente in formato 24x24 cm scannerizzata poi a 600 dpi, il programma si blocca.



Il risultato finale e' una immagine blandamente in 3D perche' manda tutta l'esagerazione verticale che fornisce un vero stereoscopio

ps. per scaricare la stereocoppie si puo' andare su Geoscopio, selezionare l'area di interesse, cliccare su Fotogrammi per esempio 2013 (non tutti gli anni hanno la copertura completa del territorio!!) , poi con lo strumento Freccia +, si seleziona uno squadro di un fotogramma e si clicca su visualizza fotogramma. Si ripete il tutto con il fotogramma vicino ed abbiamo una stereocoppia giocando con  poi con i livelli di zoom per avere il dettaglio desiderato. Per rimuovere la finestra della legenda si puo' premere il tasto h (il visualizzatore di immagini e' IIPImage)


Fauna caldinese - Nido

Nido, forse di passero, caduto a terra








giovedì 26 maggio 2016

Bossac: extra argument found su IDE Arduino Due

E' un po ' di tempo che non usavo Arduino Due ed la proma novita' e che per compilare i programmi nella IDE 1.6.8 si devono scaricare le definizioni della scheda da Strumenti/Scheda/Gestore Scheda selezionando Arduino SAM Boards 32 bits ARM Cortex-M3


Una volta modificato in questo modo l'ambiente di sviluppo ho provato a compilare uno sketch ma sia su Windows che Linux compariva lo stesso messaggio

bossac : extra argument found

il problema risiede nello switch w e si risolve andando nel file

/home/linnocenti/.arduino15/packages/arduino/hardware/sam/1.6.8/platform.txt

si va quindi alla riga 106 e si elimina dopo -w il codice {upload.verify}. La riga corretta e' la seguente

tools.bossac.upload.pattern="{path}/{cmd}" {upload.verbose} --port={serial.port.file} -U {upload.native_usb} -e -w -b "{build.path}/{build.project_name}.bin" -R
in questo modo la compilazione arriva al termine correttamente

mercoledì 25 maggio 2016

Troubleshooting Mysql

Questa volta sono stato contattato perche' un server, che in passato ho montato ed amministrato direttamente per un annetto, aveva iniziato a fare le bizze.
Una rapida occhiata con htop ha mostrato che la macchina era in picco con il processore tra 85% e 100% e con tutti i processi principali occupati da Mysql



Htop


Una volta capito il responsabile c'era da capire il motivo. Il primo sospetto, visto il tempo di esecuzione di ogni processo, era che le query fossero talmente lunghe e frequenti da sovrapporsi nel tempo da creare un effetto a cascata che rendeva dopo un po' il tempo il server non utilizzabile

Per definire il problema ho usato il comando

mysqladmin -u root -p -i 1 processlist
in questo modo si ottiene una lista dei processi in corso sul server MySql 
evidenziando anche la query ed il tempo di esecuzione di ciascuna query.
In questo modo si ottiene una  lista dei processi in corso sul server Mysql evidenziando anche la query in corso ed il tempo di esecuzione. Impostando il parametro i ad 1 si ottiene i dati relativi ad 1 secondo di tempo. Nel parziale listato si e' visto che in un secondo venivano eseguite oltre 120 query complesse con il tempo in secondo impiegato da ciascun  thread
------------------------------------------------------------------------
+-------+----------------------+--------------------------+--------+---------+------+----------------------+------------------------------------------------------------------------------------------------------+----------+
| Id    | User                 | Host                     | db     | Command | Time | State                | Info                                                                                                 | Progress |
+-------+----------------------+--------------------------+--------+---------+------+----------------------+------------------------------------------------------------------------------------------------------+----------+
| 12296 | utente               | localhost                | db | Query   | 2    | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.f     
                                        WHERE NET = (SELECT NET.NET_ID
         | 0.000    |
| 12299 | utente               | localhost                | db | Query   | 16   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.f     
                                        WHERE NET = (SELECT NET.NET_ID
         | 0.000    |
| 12303 | utente               | localhost                | db | Query   | 4    | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.f     
                                        WHERE NET = (SELECT NET.NET_ID
         | 0.000    |
| 12308 | utente               | localhost                | db | Query   | 33   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.f     
                                        WHERE NET = (SELECT NET.NET_ID
         | 0.000    |
| 13392 | utente               | localhost                | db | Query   | 35   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.E     
                                        WHERE NET = (SELECT NET.NET_ID
                     | 0.000    |
| 13396 | utente               | localhost                | db | Query   | 11   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.E     
                                        WHERE NET = (SELECT NET.NET_ID
                     | 0.000    |
| 13399 | utente               | localhost                | db | Query   | 5    | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.E     
                                        WHERE NET = (SELECT NET.NET_ID
                     | 0.000    |
| 13403 | utente               | localhost                | db | Query   | 26   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.E     
                                        WHERE NET = (SELECT NET.NET_ID
                     | 0.000    |
| 13411 | utente               | localhost                | db | Query   | 40   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.E     
                                        WHERE NET = (SELECT NET.NET_ID
                     | 0.000    |
| 14407 | utente               | localhost                | db | Query   | 45   | Sending data         | SELECT AVG(VAL) AS VAL
                                        FROM db.E     
                                        WHERE NET = (SELECT NET.NET_ID
                     | 0.000    |
| 14411 | utente               | localhost                | db | Query   | 54   | Sending data         | SELECT AVG(VAL) AS VAL 

------------------------------------------------------------------------
ovviamente un traffico che ha messo in ginocchio il server.

Un metodo alternativo e' quello di installare mytop (un clone del comando top specifico per MySql)
Il comando si puo' lanciare passando i parametri da file di configurazione oppure tramite switch sulla linea di comando. Ho usato la seconda strada

mytop -u root -d databasename --prompt
(e' necessario inserire il nome del database, non accetta il monitoraggio di tutti i db)

Mytop del server quando sono state interrotte le query lunghe

In questa modalita' vengono mostrate le query per secondo e le query lente con il traffico di rete











lunedì 23 maggio 2016

Shodan

Shodan.io e' un sito web che propone un database di indirizzi di dispositivi IOT che sono esposti su IP pubblici. La ricerca puo' essere fatta in modo geografico (per esempio city:Florence, country:"IT", con le coordinate geografiche mediante geo:42.9693,-74.1224 ) oppure con dati derivanti per servizi esposti  (port:443 per esempio) e sistema operativo (os:linux) o network (net:150.xxx.xxx.xxx/16)
Le query si possono eseguire anche senza da URL con per esempio http://www.shodanhq.com/search?q=port%3A5632 ma per le query avanzate e' necessario essere loggati



Non e' niente che non si possa fare anche tramite Google ma e' sicuramente un sistema piu' semplice e speditivo. Con uno strumento del genere ovviamente si puo' giocare oppure fare danni (mirando delle macchine con vulnerabilita'). La cosa ancora piu' strana e' che ci sia qualcuno che espone su Internet con indirizzi pubblici dei dispositivi privi di password o con credenziali di default
Giusto per vedere se la cosa funziona ho provato a cercare delle webcam.

Schermata di amministrazione

Webcam
Nell'ultimo caso sono incappato in una camera con il sistema di amministrazione ed i comandi PTZ perfettamente accessibili, probabilmente una camera Maygion o simile. Ovviamente sono uscito senza fare danni ma per esperienza,avendo anche la possibilita' di modificare il firmware, avrei potuto rendere inservibile il sistema.
 

ps. non mi e'stato possibile rintracciare il proprietario della camera per avvisarlo.







venerdì 20 maggio 2016

Track Object con OpenCv ed IPCamera PTZ in realtime

L'idea alla base di questo post e' quella di usare OpenCV per riconoscere un oggetto in stream video da una IP camera dotata di movimento PTZ in modo che l'oggetto rimanga sempre al centro dell'inquadratura

Come camera ho ripreso la Maygion (gia' vista qui, un modello piuttosto vecchiotto da 0.3MPx) e l'oggetto da tracciare era la scatola a sinistra

Tramite OpenCV viene catturato lo stream video della camera IP, l'immagine viene elaborato alla ricerca del colore verde calcolando il centroide, a questo punto si calcola la posizione del centroide rispetto al centro dell'immagine e si pilota la camera per posizionare il centroide piu' vicino possibile al centro dell'immagine.


Questo e' il video della prova

Per rendere le cose piu' semplici, in questo primo esempio viene mostrato come catturato il flusso derivante dalla camera Ip mediante OpenCV. Di fatto il sistema e' praticamente a quello che si applica ad una webcam con la sola differenza che si deve passare l'url dello stream (variabile da produttore a produttore e da modello a modello)
---------------------------------------------
import numpy as np
import cv2

vcap = cv2.VideoCapture("http://192.168.43.207:81/videostream.cgi?stream=0&usr=admin&pwd=admin")


if vcap.isOpened():

rval,frame = vcap.read()
print "acquisito"
else:
rval = False
print "non acquisito"

while rval:

cv2.imshow("preview",frame)
rval,frame = vcap.read()
key = cv2.waitKey(20)
if key == 27:
break
cv2.destroyWindow("preview")
---------------------------------------------

Questo e' invece il codice per pilotare il movimento della camera (in pratica si mandano comandi via stringa http)
---------------------------------------------
import time
try:
    from urllib.request import urlopen
except ImportError:
    from urllib2 import urlopen

ur = "http://192.168.43.207:81/"


direzione = "up"



html=urlopen(ur+"moveptz.xml?dir="+direzione+"&user=admin&password=admin")

time.sleep(0.1)
html=urlopen(ur+"moveptz.xml?dir=stop&user=admin&password=admin")
---------------------------------------------

Per tracciare il movimento della sagoma verde ho ripreso, modificandolo leggermente, questo esempio di Conan Zhao and Simon D. Levy
---------------------------------------------
import cv2
import numpy as np
import time

try:

    from urllib.request import urlopen
except ImportError:
    from urllib2 import urlopen

ur = "http://192.168.43.207:81/"



# For OpenCV2 image display

WINDOW_NAME = 'GreenBallTracker' 

def track(image):

    # Blur the image to reduce noise
    blur = cv2.GaussianBlur(image, (5,5),0)

    # Convert BGR to HSV

    hsv = cv2.cvtColor(blur, cv2.COLOR_BGR2HSV)

    # Threshold the HSV image for only green colors

    lower_green = np.array([40,70,70])
    upper_green = np.array([80,200,200])

    # Threshold the HSV image to get only green colors

    mask = cv2.inRange(hsv, lower_green, upper_green)
    
    # Blur the mask
    bmask = cv2.GaussianBlur(mask, (5,5),0)

    # Take the moments to get the centroid

    moments = cv2.moments(bmask)
    m00 = moments['m00']
    centroid_x, centroid_y = None, None
    if m00 != 0:
        centroid_x = int(moments['m10']/m00)
        centroid_y = int(moments['m01']/m00)

    # Assume no centroid

    ctr = (-1,-1)

    # Use centroid if it exists

    if centroid_x != None and centroid_y != None:

        ctr = (centroid_x, centroid_y)


        # Put black circle in at centroid in image

        cv2.circle(image, ctr, 4, (255,0,0))
x_schermo = centroid_x-320
y_schermo = centroid_x-240
print "X="+str(x_schermo) + " ; Y=" + str(y_schermo)
tempo = 0.02
if (x_schermo > 0): 
html=urlopen(ur+"moveptz.xml?dir=left&user=admin&password=admin")
time.sleep(tempo)
html=urlopen(ur+"moveptz.xml?dir=stop&user=admin&password=admin")
else: 
html=urlopen(ur+"moveptz.xml?dir=right&user=admin&password=admin")
time.sleep(tempo)
html=urlopen(ur+"moveptz.xml?dir=stop&user=admin&password=admin")


    # Display full-color image

    cv2.imshow(WINDOW_NAME, image)

    # Force image display, setting centroid to None on ESC key input

    if cv2.waitKey(1) & 0xFF == 27:
        ctr = None
    
    # Return coordinates of centroid
    return ctr

# Test with input from camera

if __name__ == '__main__':

    capture = cv2.VideoCapture(ur+"videostream.cgi?stream=0&usr=admin&pwd=admin")


    while True:


        okay, image = capture.read()


        if okay:


            if not track(image):

                break
          
            if cv2.waitKey(1) & 0xFF == 27:
                break

        else:


           print('Capture failed')

           break

---------------------------------------------

giovedì 19 maggio 2016

LCD 3.2 " Waveshare su Raspberry

Ho provato uno schermo LCD da 3.2 pollici per Raspberry da montare direttamente sul connettore a 40 pin.



Se si prova ad usare una Raspbian pura il risultato e' il seguente.. nessuna immagine tranne lo schermo bianco



Il problema e' che lo schermo comunica con la Raspberry tramite SPI e su Raspbian di default non e' montato il modulo per gestire lo schermo. Su Waveshare sono presenti tutte le istruzioni passo passo per compilare il modulo e montarlo ma la soluzione piu' diretta e' scaricarsi una distribuzione Raspbian gia' modificata



Lo schermo ha una risoluzione di 320x240 ed una volta avviato X si vede che le finestre sono tutte piu' grandi di questa dimensioni. La  soluzione per spostare l'area di visualizzazione e' quella di cliccare sinistro tenendo premuto il tasto Alt. In questo modo si puo' trascinare la finestra


Riducendo poi il font di default di X e di XTerm si riesce "quasi" a lavorarci

Con la Raspbian gia' modificata e' gia' montato anche il modulo per il touchscreen per cui e' possibile utilizzare il Raspberry con il pennino (in dotazione) od anche con le dita.In alto e' disponibile una tastiera virtuale ma e' veramente difficile usarla

Sulla sinistra dello schermo sono presenti tre pulsanti che sono connessi ai pin GPIO 12,16 e 18

Per completezza il dispositivo dello schermo si trova su /dev/fb1

Ripulire perdita di liquido da batterie / pile

Capita ogni tanto di scordarsi le batterie all'interno di un dispositivo e poi di trovare la sopresa della perdita del liquido interno con la formazione di un materiale biancastro che avvolge i contatti e la batteria stessa.
La cosa si fa tremenda quando il materiale impedisce l'uscita della batteria stessa, come nel caso in foto



Una soluzione molto semplice e  funzionante e' quella di usare l'aceto da cucina. Si sviluppa una reazione tra l'acido acetico e la patina che genera una effervescenza e che di fatto libera la pila e pulisce i contatti




Dal punto di vista strettamente chimico non saprei bene come funziona la cosa. Il reagente contenenuto nell'aceto e' sicuramente l'acido acetico. Su alcuni siti indicano che il materiale bianco e' idrossido di potassio (per pile alcaline) ma la reazione (ripresa da qui) produce acqua e non gas per cui non e' quello che ho visto io

                                     NaOH(aq) + CH3COOH(aq) → CH3COONa(aq) + H2O(l)



Insomma, non so perche' ma decisamente funziona

Ho provato lo stesso sistema su una vecchia calcolatrice Casio di fine anni 70. In questo caso l'ossido si presenta di colore verde rame (non so se e' il colore del contatto ossidato od il colore di quanto fuoriuscito dalla batteria). I contatti si presentano fortemente corrosi


In questo caso il sistema dell'aceto non funziona. Credo che in quel periodo si usassero pile con composizione differente

lunedì 16 maggio 2016

Opzioni di compilazione gcc in Eclipse e Netbeans

Un promemoria su come configurare gli switch di compilazione di gcc con Eclipse e Netbeans
In entrambi i casi si deve cliccare destro sul nome del progetto ed aprire le proprieta'

Eclipse

Optimizzazione e debug -O3 -Wall

 Include switch -I
Per lo switch -lm per esempio si scrive sono m omettendo -l
Library switch -l -L


Netbeans

Inclusione di librerie ed headers


Per aggiungere lo switch -l in Netbeans si apre Add Option/Other Option e si digita lo switch completo (al contrario di Eclipse)


collect2.exe: ld returned error 5 exit status

L'errore ld returned error 5 sul file collect2.exe  dovrebbe essere limitato all'IDE di Arduino 1.6.8 su piattaforma Windows XP.

In pratica compilando sketch perfettamente funzionanti puo' verificarsi che non venga generato l'eseguibile


L'errore e' visibile solo abilitando l'output dettagliato in fase di compilazione (da File/Impostazioni/Mostra output dettagliato) e si tratta dell'impossibilita' di rinominare un file temporaneo


L'errore non e' presente nella versione per Linux od in versioni precedenti dell'IDE di Arduino







Fractint per Arduino

Questo post prende spunto da Fractint, un software ormai molto datato (e' uno dei programmi open source che da piu' tempo e' sempre sviluppato) che permetteva il calcolo "veloce" di frattali nei calcolatori che erano privi di coprocessore matematico utilizzando gli interi per calcoli in virgola mobile.



Voluto vedere se riuscivo a fare qualcosa di simile con Arduino. Per fare cio' ho utilizzato la libreria AVRFix , libreria appositamente studiata per processori ATMEL

Il programma senza nessun tipo di ottimizzazione ha mostrato un tempo di esecuzione pari a 7581 ms. A parita' di condizioni di schermo la libreria AVRFix con il tipo di variabile short ha quasi dimezzato il tempo di calcolo per un totale di 4434 ms

Ovviamente questa non e' una soluzione generale. Basta infatti cambiare il tipo di varibiale lasciando intatto il codice per usa AVRFix che le prestazioni decadono rapidamente a 10590 ms per il tipo di varibile l a 18464 ms per il tipo di variabile ll

No AVRFix (7581 ms)
-----------------------------------------
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4
Adafruit_SSD1306 display(OLED_RESET);



#define SCREEN_WIDTH 32
#define SCREEN_HEIGHT 128

#if (SSD1306_LCDHEIGHT != 32)
#error("Height incorrect, please fix Adafruit_SSD1306.h!");
#endif

//float re_min = -2.019531;
float a = -2.003906;
float im_min = -1.275;
float re_factor = 0.019531;
float im_factor = 0.075;

unsigned long tempo;

float b;
float x,y,x_new,y_new;

int k,j,i;

void setup()   {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();

  tempo = millis();
  
  for (i=0;i<SCREEN_HEIGHT;i++)
     {
     a = a + re_factor;
     b = im_min;
     for (j=0;j<SCREEN_WIDTH;j++)
      {
        b=b+im_factor;
        x = 0;
        y = 0;
      
        for (k=0;k<=64;k++)
            {
            x_new = (x*x)-(y*y)+a;
            y_new = (2*x*y)+b;
            if (((x_new*x_new)+(y_new*y_new))>4)
                {
                display.drawPixel(i, j, WHITE);
                //
                break;
                }
       x = x_new;
       y = y_new;
       } 
      }
     }
Serial.println(millis()-tempo);
display.display();
}


void loop() {
}
-----------------------------------------

AVRFix (tipo short 4728 ms)
-----------------------------------------
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4

Adafruit_SSD1306 display(OLED_RESET);
extern "C"{
  #include <avrfix.h>
}

#define WIDTH 32
#define HEIGHT 128
lfix_t Im = ftosk(-1.275);
lfix_t DRe = ftosk(0.02);
lfix_t DIm = ftosk(0.075);
lfix_t A = ftosk(-2.0);
lfix_t B = ftosk(0.0);
lfix_t DUE = ftosk(2.0);
lfix_t X,Y;
lfix_t XN,YN,YNN;
lfix_t tt;
int j,k,i;
float test;

unsigned long tempo;

void setup() {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();
  tempo = millis();
  
  for (i=0;i<HEIGHT;i++)  
    {
      A=A+DRe;
      B=Im;
      for (j=0;j<WIDTH;j++)
        {
          B=B+DIm;
          X = ftosk(0.0); 
          Y = ftosk(0.0);
          for (k=0;k<=64;k++)
              {
               XN = smulsk(X, X)-smulsk(Y, Y)+A;
               YN = smulsk(X, Y);
               YNN = smulsk(X, Y);
               YN = smulsk(DUE, YNN) + B;
               tt = smulsk(XN,XN)+smulsk(YN,YN);
               test = sktof(tt);
               if (test > 4.0)
                    {
                        display.drawPixel(i, j, WHITE);
                        break;
                    }
               X=XN;
               Y=YN;
              }
        }
    }
  

  Serial.println(millis()-tempo);
  display.display();
  
}

void loop() {
}
-----------------------------------------

AVRFix (tipo l 10590 ms)
-----------------------------------------
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4

Adafruit_SSD1306 display(OLED_RESET);
extern "C"{
  #include <avrfix.h>
}

#define WIDTH 32
#define HEIGHT 128
//lfix_t Re = ftolk(-2.00390625); 
lfix_t Im = ftosk(-1.275);
lfix_t DRe = ftok(0.019531);
lfix_t DIm = ftok(0.075);
lfix_t A = ftok(-2.00390625);
lfix_t B = ftok(0.0);
lfix_t DUE = ftok(2.0);
lfix_t X,Y;
lfix_t XN,YN,YNN;
lfix_t tt;
int j,k,i;
double test;

unsigned long tempo;

void setup() {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();
  tempo = millis();
  
  for (i=0;i<HEIGHT;i++)  
    {
      A=A+DRe;
      B=Im;
      for (j=0;j<WIDTH;j++)
        {
          B=B+DIm;
          X = ftok(0.0); 
          Y = ftok(0.0);
          for (k=0;k<=64;k++)
              {
               XN = mulk(X, X)-mulk(Y, Y)+A;
               YN = mulk(X, Y);
               YNN = mulk(X, Y);
               YN = mulk(DUE, YNN) + B;
               tt = mulk(XN,XN)+mulk(YN,YN);
               test = ktod(tt);
               if (test > 4.0)
                    {
                        display.drawPixel(i, j, WHITE);
                        break;
                    }
               X=XN;
               Y=YN;
              }
        }
    }
  

  Serial.println(millis()-tempo);
  display.display();
  
}

void loop() {

}
-----------------------------------------

AVRFix (tipo ll 18464 ms)
-----------------------------------------
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

#define OLED_RESET 4

Adafruit_SSD1306 display(OLED_RESET);
extern "C"{
  #include <avrfix.h>
}

#define WIDTH 32
#define HEIGHT 128
//lfix_t Re = ftolk(-2.00390625); 
lfix_t Im = ftolk(-1.275);
lfix_t DRe = ftolk(0.019531);
lfix_t DIm = ftolk(0.075);
lfix_t A = ftolk(-2.00390625);
lfix_t B = ftolk(0.0);
lfix_t DUE = ftolk(2.0);
lfix_t X,Y;
lfix_t XN,YN,YNN;
lfix_t tt;
int j,k,i;
double test;

unsigned long tempo;

void setup() {
  Serial.begin(115200);
  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);  
  display.display();
  delay(100);
  display.clearDisplay();
  tempo = millis();
  
  for (i=0;i<HEIGHT;i++)  
    {
      A=A+DRe;
      B=Im;
      for (j=0;j<WIDTH;j++)
        {
          B=B+DIm;
          X = ftolk(0.0); 
          Y = ftolk(0.0);
          for (k=0;k<=64;k++)
              {
               XN = lmullk(X, X)-lmullk(Y, Y)+A;
               YN = lmullk(X, Y);
               YNN = lmullk(X, Y);
               YN = lmullk(DUE, YNN) + B;
               tt = lmullk(XN,XN)+lmullk(YN,YN);
               test = lktod(tt);
               if (test > 4.0)
                    {
                        display.drawPixel(i, j, WHITE);
                        
                        break;
                    }
               X=XN;
               Y=YN;
              }
        }
    }
  

  Serial.println(millis()-tempo);
  display.display();
  
}

void loop() {

}
-----------------------------------------


Disclaimer: il titolo di questo post non e' da intendersi come un plagio verso il nome del piu' noto software Fractint ma piuttosto come un tributo ad un programma che usavo sui 386 nei primi anni 90



Roccia lunare 14310

----------------------------------------------
Aggiornamento importante a questo post
----------------------------------------------
Ho avuto la fortuna di poter osservare da vicino (e come geologo e' una esperienza particolare) un frammento di roccia lunare ritagliato da campione piu' grande prelevato dalla missione Apollo 14


La roccia deriva dal punto G della seconda escursione lunare (scheda completa qui)



Le quattro facce (fotografabili) del campione. Di particolare interesse e' la foto sottostante perche' questo lato era quello esposto (gli altri lati derivano dal taglio effettuato a terra) e presenti segni dei micro meteoriti che hanno "lavorato" la superficie






Si tratta di una roccia decisamente giovane (circa 4 milioni di anni)  per lo standard lunare.
L'esame petrografico (originale a questo link) mostra un basalto feldspatico a grana fine e viene definito come composizione KREEP (Potassio, Terre Rare e Fosforo), ultimo resto del frazionamento del magma che ha formato la crosta lunare e che essendo composto da elementi incompatibili per entrare nella struttura cristallina rimangono isolati all'interno del magma stesso (fonte wikipedia, wikipedia inglese). I magmi ricchi in KREEP si sarebbero segregati tra la crosta ed il mantello lunare. Sono presenti fenocristalli di plagioclasi e pirosseni negli interstizi

Questo era l'aspetto della roccia lunare prima di essere sezionata per le analisi, piu' o come come la hanno vista Shepard e Mitchell (in origine era mezza sepolta dalla regolite dalla regolite). Non ci sono foto realizzate sulla Luna di questo esemplare


Come detto non esistono foto dirette del punto di campionamento. Questa dovrebbe essere nella zona compresa tra la stazione C1 ed H (nel mezzo tra le due c'e' la stazione G) ripresa da Project Apollo Archive

AS14-68-9471HR.jpg




venerdì 13 maggio 2016

Mandelbrot con GMP

Avevo gia' utilizzato la libreria GMP in questo post per il calcolo di PiGreco
Questa volta la utilizzero' per effettuare zoom estremi dell'insieme di Mandelbrot grazie alla sua caratteristica di poter utilizzare una precisione arbitraria ovvero di poter settare il numero di bit che rappresentano una quantita' in virgola mobile. Il tipo float cambia a seconda da compilatore a compilatore e da processore a processore ma diciamo che in generale puo' essere considerato a 32 bit, il che va bene quando il numero di cifre significative non e' molto esteso

Per mostrare la differenza nell'uso di GMP rispetto al solo tipo float ho scritto due programmi ed ho effettuato due zoom dell'insieme di Mandelbrot

Zoom 1:
In questo caso il risultato e' identico
GMP

Float


Zoom 2
Riducendo la finestra di visualizzazione di centesimo rispetto allo zoom precedente GMP continua a lavorare in modo corretto mentre il programma scritto in C++ puro "esplode"
GMP


Float

Naturalmente tutto cio' ha un prezzo in termini di memoria consumata (la precisione e' vero che e' arbitraria ma deve comunque fare i conti con la ram disponibile) ed il tempo di calcolo (ovviamente piu' lento con GMP).Curiosamente provando a lanciare questo programma in macchina virtuale VMWare, il programma mi diceva che non aveva abbastanza memoria disponibile e sono passato ad un calcolatore reale

esempio con l'utilizzo di GMP
----------------------------------------------------------

#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <iostream>

using namespace std;

#include <pngwriter.h>

#define SCREEN_WIDTH 1280
#define SCREEN_HEIGHT 960


int main(void){

mpf_t A,B, Im, DRe, DIm,test;
mpf_t X,Y,XN,YN,XX,YY,Xs,Ys;
float test2;

int j,k,i;

long int precision = 256;
double r;




mpf_set_default_prec(precision);

mpf_init2(A,precision);
mpf_init2(B,precision);
mpf_init2(Im,precision);
mpf_init2(DRe,precision);
mpf_init2(DIm,precision);
mpf_init2(test,precision);
mpf_init2(X,precision);
mpf_init2(Y,precision);
mpf_init2(XN,precision);
mpf_init2(YN,precision);
mpf_init2(XX,precision);
mpf_init2(YY,precision);
mpf_init2(Xs,precision);
mpf_init2(Ys,precision);

// punto vertice
// Re = -0.74364085
// Im = 0.13182733
// Diametro = 0.00012068
// Zoom = 25497

pngwriter png(SCREEN_WIDTH,SCREEN_HEIGHT,0,"mandelbrot2.png");
double iterazioni = 1024;

mpf_set_d(A,-0.74364085);
//mpf_set_d(A,-0.74370119);
mpf_set_d(B,0.13182733);
mpf_set_d(Im,0.13182733);
mpf_set_d(DRe,0.000000009428125);
mpf_set_d(DIm,0.00000001257083);



mpf_set_d(X,0.0);
mpf_set_d(Y,0.0);
mpf_set_d(Xs,0.0);
mpf_set_d(Ys,0.0);
mpf_set_d(XX,0.0);
mpf_set_d(YY,0.0);
mpf_set_d(XN,0.0);
mpf_set_d(YN,0.0);

cout << "Inizio" << endl;

for (i=0;i<SCREEN_HEIGHT;i++)
    {
    mpf_add(A,A,DRe);
    mpf_set(B,Im);
    for (j=0;j<SCREEN_WIDTH;j++)
    {
    mpf_add(B,B,DIm);
    mpf_set_d(X,0.0);
    mpf_set_d(Y,0.0);

    for (k=0;k<=iterazioni;k++)
        {
        mpf_pow_ui(Xs,X,2);
        mpf_pow_ui(Ys,Y,2);
        mpf_sub(XN,Xs,Ys);
        mpf_add(XN,XN,A); //XN = (x*x)-(y*y) +a
        mpf_mul(YN,X,Y);
        mpf_mul_ui(YN,YN,2);
        mpf_add(YN,YN,B); //YN = 2*X*Y

        mpf_pow_ui(XX,XN,2);
        mpf_pow_ui(YY,YN,2);
        mpf_add(test,XX,YY);// test = (XN*XN)+(YN*YN)
        test2 = mpf_get_d(test);
        if (test2 > 4.0)
            {
            r = k%2;
            png.plot(j,i,r,r,r);
            break;
            }
        mpf_set(X,XN);
        mpf_set(Y,YN);
        }
    }
    }

png.close();

mpf_clear(A);
mpf_clear(B);
mpf_clear(A);
mpf_clear(Im);
mpf_clear(DRe);
mpf_clear(DIm);
mpf_clear(X);
mpf_clear(Y);
mpf_clear(Xs);
mpf_clear(Ys);
mpf_clear(XX);
mpf_clear(YY);
mpf_clear(XN);
mpf_clear(YN);
mpf_clear(test);

cout << "Fine" << endl;

return 0;
}


----------------------------------------------------------
esempio in C++ puro
----------------------------------------------------------
#include <pngwriter.h>

#define SCREEN_WIDTH 1280
#define SCREEN_HEIGHT 960



// punto vertice
// Re = -0.74364085
// Im = 0.13182733
// Diametro = 0.00012068
// Zoom = 25497

double iterazioni = 1024;

float a = -0.74364085;
float b = 0.13182733;
float im = 0.13182733;

float dre = 0.000000009428125;
float dim = 0.00000001257083;


double r;

float x,y,x_new,y_new;

int test;

int k,j,i;

int keypress = 0;


int main() {

pngwriter png(SCREEN_WIDTH, SCREEN_HEIGHT,0,"mand_norm2.png");


for (i=0;i<SCREEN_HEIGHT;i++)
 {
 a=a+dre;
 b=im;
 for (j=0;j<SCREEN_WIDTH;j++)
  {
  b = b+dim;

  x = 0;
  y = 0;
  test = 0;

  for (k=0;k<=iterazioni;k++)
   {
   x_new = (x*x)-(y*y)+a;
   y_new = (2*x*y)+b;
   if (((x_new*x_new)+(y_new*y_new))>4)
    {
    r = k%2;
    png.plot(j,i, r, r, r);
    break;
    }
   x = x_new;
   y = y_new;
   }

  }

 }
png.close();
return(0);

}

----------------------------------------------------------