La trasformata DCT applicata alle immagini
description
Transcript of La trasformata DCT applicata alle immagini
Metodi del Calcolo Scientifico
Progetto 2: La trasformata DCT applicata alle
immagini
Federico Valentini, 726584Sergio Cattaneo, 726786
Appello del 16 Giugno 2014
1 Consegna
Questo progetto ha lo scopo di studiare la trasformata DCT applicata alleimmagini in toni di grigio in formato bitmap.Si richiede la scrittura di un software in un linguaggio di programmazione ascelta che esegua i seguenti compiti:
• leggere una immagine in formato BMP in toni di grigio (1 byte perpixel)
• calcolarne la DCT2 (DCT bidimensionale) usando obbligatoriamentela libreria FFTW
• disegnare un istogramma dei coefficienti di Fourier (con la possibilitadi scegliere il numero dei bins)
• chiedere all’utente un valore soglia e mettere a zero le ampiezze dellefrequenze al di sotto (in valore assoluto) di questa soglia
• calcolare la IDCT2 dei nuovi coefficienti e disegnare sia l’immagineoriginale sia l’immagine ottenuta dopo il taglio dei coefficienti in mododa poterle confrontare
1
2 Teoria matematica
2.1 DCT - Trasformata discreta del coseno
Strettamente legata alla trasformata discreta di Fourier, la DCT e unatrasformazione separabile linearmente che per un’immagine A in input resti-tuisce in output:
Bpq = αpαq∑M−1
m=0
∑N−1n=0 Amncos
π(2m+1)p2M cosπ(2n+1)q
2N
dove:
0 ≤ p ≤M − 10 ≤ q ≤ N − 1
αp =
{ 1√M
p = 0√2M 1 ≤ p ≤M − 1
αq =
{ 1√N
q = 0√2N 1 ≤ q ≤ N − 1
2.2 IDCT - Antitrasformata discreta del coseno
Amn =∑M−1
p=0
∑N−1q=0 αpαqBpqcos
π(2m+1)p2M cosπ(2n+1)q
2N
dove:
0 ≤ m ≤M − 10 ≤ n ≤ N − 1
2
3 Implementazione
3.1 FFTW
FFTW[1] e una libreria nativa in C, open-source, per il calcolo della trasfor-mata discreta di Fourier.La funzione principale che abbiamo utilizzato e stata fftw plan r2r. Questafunzione riceve in input la dimensione della matrice e il tipo di trasformatada applicare, per esempio la DCT e rappresentata dalla stringa FFTW REDFT10mentre la IDCT da FFTW REDFT01.
Da sottolineare che una volta calcolate DCT ed IDCT, sono state effettuatele seguenti normalizzazioni:
• Divisione di tutta la matrice risultate da DCT per un fattore di
4√
ROWS2
√COLS
2
• Ulteriore divisione sulla prima riga e prima colonna per un fattore di√2
• Dopo la manipolazione dei coefficienti, queste due normalizzazioni ven-gono annullate
• Infine, dopo aver applicato la IDCT, divisione per un coefficiente di4 x ROWS x COLS
3.2 OpenCV
Per gestire le immagini in input e generare le finestre con i risultati, ci siamoserviti di OpenCv[2], libreria open-source per la computer vision implemen-tata in C, C++ e Python.Questa libreria e fortemente supportata da una community di piu di 47000utenti, costatemente aggiornata (ultima versione del 25 Aprile 2014) e disponi-bile su piu piattaforme (Windows, Linux/Mac, Android e iOS)
3.3 GSL: GNU Scientific Library
GSL[3] e una libreria open-source di calcolo numerico per C/C++, che com-prende fra le altre cose funzioni per numeri complessi, algebra lineare, gen-erazione di numeri casuali; tuttavia per questo programma ci siamo servitisolo delle primitive relative al calcolo dell’istogramma.
3
4 Sviluppo
4.1 Dettagli macchina
Per la realizzazioni di questi test ci siamo affidati ad una macchina LenovoIdeaPad S510P con processore Intel Core i7-4500U (1.80 / 3.00 GHZ, 2 core,4 MB CACHE L3) con 4 GB di RAM.Il sistema operativo utilizzato e Manjaro 0.8.9 (una fork di ArchLinux).
4.2 Strutture dati e funzioni utilizzate
Per il salvataggio e l’elaborazione dei dati, ci siamo serviti di due vettori:
in = (double*) fftw_malloc(sizeof(double) * ROWS * COLS);
out = (double*) fftw_malloc(sizeof(double) * ROWS * COLS);
Nel primo vengono salvati i pixels che compongono l’immagine, mentrenel secondo i coefficienti della trasformata discreta coseno.Per la stampa a video delle immagini (prima e dopo la trasformazione) ab-biamo utilizzato la struttura dati IplImage presente nella libreria OpenCV.Dopo aver letto l’immagine il primo passo computazionale da fare e lacreazione di un piano
fftw_r2r_kind kind[RANK] = {FFTW_REDFT10, FFTW_REDFT10};
int dimentions[2] = {ROWS, COLS};
p = fftw_plan_r2r(RANK, dimentions, in, out, kind, FFTW_MEASURE);
dove vengono specificati dimensione della matrice, le strutture dati nec-essarie al calcolo e il tipo di trasformata da calcolare per ogni direzione,in particolare, come gia accennato, FFTW REDFT10 per la DCT mentreFFTW REDFT01 per la IDCT.
4
5 Risultati ottenuti
Abbiamo deciso di riportare le nostre osservazioni su diverse tipologie diimmagini fornite.
5.1 Immagine di grandi dimensioni
Lo scopo principale di questa esecuzione e stato semplicemente una misuradel tempo di computazione nel caso limite di immagini molto grandi.
Dimensione: 7216 x 5412 pixels. Tempo di calcolo: circa 25 secondi
5
5.2 Immagine artificiale
Nel secondo test abbiamo utilizzato un’immagine generata al computer, pertestare il comportamento di Fourier (e derivate) in immagini di questa tipolo-gia.
Valore soglia: 30. Con questo valore si iniziano a notare delle prime distor-sioni anche se l’immagine risulta ancora ben definita.
Valore soglia: 100. Lo sfondo inizia a non essere piu uniforme come l’immagineoriginale mentre gli oggetti in primo piano cominciano a diventare sempremeno definiti.
6
Valore soglia: 300. L’immagine inizia a non essere piu riconoscibile e pre-senta una qualita estremamente degradata rispetto all’immagine originale.
7
5.3 Immagine reale scura
Nel terzo test abbiamo analizzato un’immagine composta per la maggiorparte da un’area omogenea molto scura attorno ad una zona luminosa e piucontrastata.
Valore soglia: 40. Solo la zona del fumo prodotto dai fuochi d’artificio risentepiu visibilmente del taglio effettuato.
8
Valore soglia: 100. L’immagine inizia a degradarsi in maniera piu uniforme,pur mantenendo il soggetto principale della foto abbastanza distinguibileancora.
5.4 Immagine reale chiara e contrastata
Come ultimo test abbiamo scelto un’immagine reale senza ampie zone omo-genee scure.
Valore soglia: 50. L’immagine risulta molto ancora molto simile all’immagineoriginale.
9
Valore soglia: 100. La distorsione risulta essere piu evidente anche sel’immagine ancora risulta essere ben definita
Valore soglia: 500. Immagine degradata e non distinguibile.
10
A Appendice: codice sorgente applicazione
#include <math.h>
#include <stdio.h>
#include <fftw3.h>
#include <opencv2/core/core.hpp>
#include <opencv/highgui.h>
#include <opencv/cv.h>
#include <gsl/gsl_histogram.h>
#define RANK 2 // 2 dimensioni
int main();
double greatestCoeff(double* out, int size){
int i = 1;
double max = out[0];
for (i; i<size; i++){
if (out[i] > max) max = out[i];
}
return max;
}
double smallestCoeff(double* out, int size){
int i = 1;
double min = out[0];
for (i; i<size; i++){
if (out[i] < min) min = out[i];
}
return min;
}
IplImage* myHist(double* coeff, int size, int numBins, float
scaleX, float scaleY){
int i;
double min = smallestCoeff(coeff, size);
double max = greatestCoeff(coeff, size) + 1;
gsl_histogram* hist = gsl_histogram_alloc(numBins);
gsl_histogram_set_ranges_uniform(hist, min, max);
for (i = 0; i < size; i++){
gsl_histogram_increment(hist, coeff[i]);
}
11
double M = 256*scaleX;
double N = 128*scaleY;
double base = M / numBins;
IplImage* imgHist = cvCreateImage(cvSize(M, N), 8, 3);
cvZero(imgHist);
double histMax = max;
for(i = 0; i < numBins; i++){
cvRectangle(imgHist, cvPoint(i*base, N),
cvPoint((i+1)*base-1, N -
gsl_histogram_get(hist,
i)/N),
cvScalar(255, 255, 255, 255),
CV_FILLED, 8, 0);
}
printf("\n");
printf("Istogramma in forma testuale: \n\n");
gsl_histogram_fprintf (stdout, hist, "%g", "%g");
return imgHist;
}
int main(int argc, char** argv){
printf("\nINIZIO\n\n");
IplImage* img = cvLoadImage(argv[1], 0);
printf("Immagine in input: %s\n", argv[1]);
int ROWS = img->height;
int COLS = img->width;
long N = ROWS*COLS;
int a, b;
// Allocazione delle risorse
double *in, *out;
in = (double*) fftw_malloc(sizeof(double) * ROWS * COLS);
out = (double*) fftw_malloc(sizeof(double) * ROWS * COLS);
unsigned long i, j, n;
printf("Dimensioni immagine: ");
printf("%dx%d pixels\n\n", ROWS, COLS);
//Creazione piano dct2
fftw_plan p;
fftw_r2r_kind kind[RANK] = {FFTW_REDFT10, FFTW_REDFT10};
int dimentions[2] = {ROWS, COLS}; //dimentions
12
p = fftw_plan_r2r(RANK, dimentions, in, out, kind,
FFTW_MEASURE);
//Popolazione vettore da matrice
unsigned char value;
for (i=0; i<ROWS; i++){
for (j=0; j<COLS; j++){
value = CV_IMAGE_ELEM(img, unsigned char, i,
j);
in[(i*COLS) + j] = value;
}
}
printf("Calcolo dct2\n\n");
fftw_execute(p);
for (n = 0; n < (ROWS * COLS); n++){
out[n] /= 4 * sqrt(ROWS/2.0) * sqrt(COLS/2.0);
}
for (n = 0; n < COLS; n++){
out[n] /= sqrt(2.0);
}
for (n = 0; n < (ROWS * COLS); n = n + COLS){
out[n] /= sqrt(2.0);
}
int numBins;
printf("Inserire il numero di bins in cui si vuole
suddividere \nl’istogramma dei coefficienti: ");
scanf("%d", &numBins);
IplImage* imgHist = myHist(out, ROWS*COLS, numBins, 2.0,
2.0);
cvNamedWindow("Istogramma dei
coefficienti",CV_WINDOW_AUTOSIZE);
cvShowImage("Istogramma dei coefficienti", imgHist);
printf("\nPremere un tasto per continuare\n\n");
cvWaitKey(0);
double min = smallestCoeff(out, ROWS*COLS);
printf("Coefficiente pi piccolo: %f\n", min);
13
double max = greatestCoeff(out, ROWS*COLS);
printf("Coefficiente pi grande: %f\n\n", max);
double cut;
printf("Inserire il valore soglia per il filtro passa alto
sui coefficienti: ");
scanf("%lf", &cut);
for (n = 0; n < (ROWS * COLS); n++){
if (fabs(out[n]) < cut){
out[n] = 0.0;
}
}
for (n = 0; n < (ROWS * COLS); n = n + COLS){
out[n] *= sqrt(2.0);
}
for (n = 0; n < COLS; n++){
out[n] *= sqrt(2.0);
}
for (n = 0; n < (ROWS * COLS); n++){
out[n] *= 4 * sqrt(ROWS/2.0) * sqrt(COLS/2.0) ;
}
// Creazione piano idct2
kind[0] = FFTW_REDFT01;
kind[1] = FFTW_REDFT01;
p = fftw_plan_r2r(RANK, dimentions, out, in, kind,
FFTW_MEASURE);
printf("\nCalcolo idct2\n\n");
fftw_execute(p);
int Z = 4*ROWS*COLS;
for (n = 0; n < (ROWS * COLS); n++){
in[n] /= (int) Z;
if(in[n] < 0) in[n] = 0;
if(in[n] > 255) in[n] = 255;
}
IplImage* img2 =
cvCreateImage(cvSize(COLS,ROWS),IPL_DEPTH_8U,1);
14
for (i=0; i<ROWS; i++){
for (j=0; j<COLS; j++){
value = in[(i*COLS) + j];
CV_IMAGE_ELEM(img2, unsigned char, i, j) =
(unsigned char) value;
}
}
printf("Riselezionare l’istogramma e premere un tasto
per\nvisualizzare le immagini.\n");
cvWaitKey(0);
cvNamedWindow("ORIGINALE",CV_WINDOW_AUTOSIZE);
cvShowImage("ORIGINALE", img);
cvNamedWindow("MODIFICATA",CV_WINDOW_AUTOSIZE);
cvShowImage("MODIFICATA", img2);
printf("Premere un tasto per terminare\n");
cvWaitKey(0);
// Disallocazione delle risorse
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
cvReleaseImage(&img);
cvReleaseImage(&img2);
printf("\nFINE\n\n");
return 0;
}
15
References
[1] FFTW: http://www.fftw.org/
[2] OpenCV: http://www.opencv.org/
[3] GSL: http://www.gnu.org/s/gsl/
16