Eliminación de ruido y compresión de imágenes mediante Wavelets

23 diciembre, 2017 Garikoitz Martínez Moreno Comentarios 0 comentarios

Índice

Objetivo

Aplicar análisis wavelets de eliminación de ruido y de compresión de la información a señales provenientes de experimentos de fusión nuclear por confinamiento magnético.

Introducción

Disponemos de imágenes en bruto almacenadas en forma de matriz. Cada archivo en bruto tiene un tamaño de 5631 KB y se corresponde con una imagen de 576×385 pixels.

Figura 1 – Detalle del archivo tsECH11568_image.dat

Los archivos son imágenes de diagnóstico Thomson Scattering TJ-II y se dividen en cinco grupos.

Fondo cámara CCD (BKGND).
Plasma en corte ECRH (COFF).
Plasma con calentamiento ECRH (ECH).
Plasma con calentamiento NBI (NBI).
Luz parásita sin plasma (STRAY).

Programa desarrollado en Matlab

A continuación veamos unas capturas del programa desarrollado.

Tras añadir ruido a la imagen seleccionada se presenta la figura 4.

Figura 4 – Resultado de la adición de ruido.

Tras la eliminación de ruido de la imagen seleccionada se presentan la figura 5, 6 y 7.

Figura 5 – Resultado de la eliminación de ruido.

Figura 6 – Resultado de la eliminación de ruido. Descomposición por canales RGB.

Figura 7 – Resultado de la eliminación de ruido. Descomposición por canales RGB de coeficientes.

Tras la compresión de la imagen seleccionada se presenta la figura 8.

Figura 8 – Resultado de la compresión de la imagen.

El programa se apoya en tres carpetas para su funcionamiento como indica la figura 9.

Figura 9 – Carpetas usadas por el programa realizado.

TS_files: Reservado para archivos de configuración del programa.
TS_saved_files: Es donde se almacenan automáticamente las imágenes resultantes del procesado.
TS_Signals: Todas las imágenes en bruto (.dat) que haya en esta carpeta el programa las lee automáticamente y las carga en un menú popup para facilitar el procesado. Si esta carpeta estuviera vacía, también se puede abrir una imagen en bruto desde el menú File > Open TS Image.

Añadir y eliminar ruido

Añadir ruido

Para la adición de ruido utilizaremos la función imnoise que nos permite añadir ruido fácilmente y de entre varios tipos.

J = imnoise(I,type,int)

Donde:

I: es la imagen a la que queremos añadir ruido.
type: el tipo de ruido.
Int: es la intensidad del ruido.
J: es la imagen con ruido.

Figura 10 – Tipos de ruido e intensidades disponibles.

Figura 11 – Resultado de añadir ruido speckle de intensidad 0.08.

Eliminar ruido

Para la eliminación de ruido utilizaremos la función de la transformada discreta 2D dwt2. A continuación calculamos el nivel de error y el umbral, descomponemos (wavedec2) el canal y umbralizamos para finalmente reconstruirlo (waverec2). Esto se hace con los tres canales y finalmente se unen para formar la imagen final.

[cA,cH,cV,cD] = dwt2(X,’wname’);

[C,S] = wavedec2(X,N,’wname’); XDEN = waverec2(C,S,’wname’)

Donde:

cA, cH, cV y cD: son los coeficientes de descomposición de la imagen.
X: es la imagen a tratar. XDEN: la imagen resultante.
C: es el vector de descomposición. S: es la matriz de seguimiento del vector.
wname: es el tipo de wavelet

A continuación veamos el código principal de eliminación de ruido.

%aR
        [CA,CH,CV,CD] = dwt2(redChannel,tipowave);
        subplot(4,4,5); imshow(CA/255);title('CA redChannel');
        subplot(4,4,6); imshow(CH/255);title('CH redChannel');
        subplot(4,4,7); imshow(CV/255);title('CV redChannel');
        subplot(4,4,8); imshow(CD/255);title('CD redChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(redChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(redChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),sorh,thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aR=waverec2(c_n,l,tipowave);
%aG
        [CA,CH,CV,CD] = dwt2(greenChannel,tipowave);
        subplot(4,4,9); imshow(CA/255);title('CA greenChannel');
        subplot(4,4,10); imshow(CH/255);title('CH greenChannel');
        subplot(4,4,11); imshow(CV/255);title('CV greenChannel');
        subplot(4,4,12); imshow(CD/255);title('CD greenChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(greenChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(greenChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),sorh,thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aG=waverec2(c_n,l,tipowave);
%aB
        [CA,CH,CV,CD] = dwt2(blueChannel,tipowave);
        subplot(4,4,13); imshow(CA/255);title('CA blueChannel');
        subplot(4,4,14); imshow(CH/255);title('CH blueChannel');
        subplot(4,4,15); imshow(CV/255);title('CV blueChannel');
        subplot(4,4,16); imshow(CD/255);title('CD blueChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(blueChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(blueChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),sorh,thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aB=waverec2(c_n,l,tipowave);
        XDEN = cat(3, aR, aG, aB);
        XDEN=uint8(XDEN);

Para medir en cierta medida la calidad de la reconstrucción de la imagen usaremos la relación señal a ruido de pico o en su término anglosajón Peak Signal to Noise Ratio (PSNR) y el error cuadrático medio cuyo término anglosajón es Mean Squared error (MSE). Cuanto mayor es el PSNR mejor es la codificación de la imagen.

Tabla 1 – Resultados obtenidos alterando el nivel.

Figura 12 – Comparativa eliminación de ruido con rbio1.1 y varios niveles.

En la tabla 1 podemos ver que el PSNR más alto se obtiene con el nivel más bajo (N=1) y en este caso la wavelet rbio1.1, pero si nos fijamos más en detalle en las imágenes de la figura 12, podemos comprobar que la eliminación de ruido empieza a ser satisfactoria a partir de N=4, por lo tanto, obviando el PSNR la mejor configuración sería rbio1.1 y N=4.

Tabla 2 – Resultados obtenidos alterando la wavelet y su momento de fuga.

Figura 13 – Canales RGB de la imagen SIN ruido -> sym8, N=4.

Figura 14 – Canales RGB de la imagen SIN ruido -> bior3.5, N=4.

El mejor PSNR lo obtiene la wavelet bior3.5 aunque el resultado visualmente no es del todo satisfactorio, ya que contiene demasiado ruido como se puede ver en la figura 14. Hay varias wavelets que han dado resultados satisfactorios pero en este caso la mejor ha sido la configuración sym8, N=4.

Compresión de imágenes

Utilizaremos la función wdencmp ya que nos permite mayor flexibilidad a la hora manejar los parámetros de compresión. Esta función sirve tanto para comprimir como para eliminar ruido, aunque nosotros la utilizaremos para lo primero. Su estructura es la siguiente:

[XC,CXC,LXC,PERF0,PERFL2] = wdencmp(Opt=lvd,Señal,Wavelet,N,THR,SORH)

[XC,CXC,LXC,PERF0,PERFL2] = wdencmp(Opt=gbl,Señal,Wavelet,N,THR,SORH,KEEPAPP)

Parámetros obtenidos:

XC: imagen resultante de la compresión.
CXC: parámetro de descomposición de la wavelet.
LXC: parámetro de descomposición de la wavelet.
PERF0: Coeficientes nulos.
PERFL2: Factor de compresión.

Parámetros elegidos:

Opt: Parámetro que indica el tipo de umbral (threshold) a utilizar. Puede ser global (gbl) o dependiente del nivel (lvd).
Señal: Matriz que contiene la señal a procesar.
Wavelet: Wavelet a utilizar.
N: Nivel de compresión.
THR: Umbral.
SORH:Tipo de umbral suave (s) o fuerte (h).
KEEPAPP: Si su valor es 1, los coeficientes de aproximación no se umbralizan. Solo lo utilizamos cuando opt=gbl.

Código utilizado para la compresión:

case 1
        opt = 'lvd';
        thr_h = [17 18];    % Horizontal thresholds.
        thr_d = [19 20];    % Diagonal thresholds.
        thr_v = [21 22];    % Vertical thresholds. 
        thr = [thr_h ; thr_d ; thr_v];
        redChannel = Imagen_Inicial(:, :, 1);
        greenChannel = Imagen_Inicial(:, :, 2);
        blueChannel = Imagen_Inicial(:, :, 3);
        %Al utilizar thr como matrix n=2, de lo contrario dará error.
        [xdr,cxd2r,lxd2r,perf0r,perfl2r] = wdencmp(opt,redChannel,tipowave,2,thr,sorh);
        [xdg,cxd2g,lxd2g,perf0g,perfl2g] = wdencmp(opt,greenChannel,tipowave,2,thr,sorh);
        [xdb,cxd2b,lxd2b,perf0b,perfl2b] = wdencmp(opt,blueChannel,tipowave,2,thr,sorh);
        xd = cat(3, xdr, xdg, xdb);
        perf0 = (1/3)*(perf0r+perf0g+perf0b);
        perfl2 = (1/3)*(perfl2r+perfl2g+perfl2b);
        xd=uint8(xd);
case 2 
        opt = 'gbl';
        [c s] = wavedec2(Imagen_Inicial,n,char(tipowave));
        thr=20;
        [xd,cxd,lxd,perf0,perfl2] = wdencmp(opt,c,s,tipowave,n,thr,sorh,1);
        xd=uint8(xd);

En la opción dependiente del nivel (lvd), el umbral (thr) es una matriz 3 por N por lo que en este caso utilizaremos N=2. Otra particularidad radica en que trabajamos directamente sobre la señal a procesar, de modo que para optimizar el resultado separamos los canales RGB de la imagen y los procesamos por separado.

En la opción global (gbl), se descompone la imagen a procesar mediante la función wavedec2 (descomposición 2D multinivel), a la que le pasamos la señal a procesar, el nivel de descomposición (N) y la wavelet. La función nos devuelve el vector de descomposición (C) y la matriz de seguimiento de cada componente (S) que serán los parámetros que usaremos en wdencmp.

[C,S] = wavedec2(Señal,N,’wname’)

Resultados obtenidos

De la tabla 3 podemos llegar a varias conclusiones. La primera de ellas es que un umbral suave obtiene mejor relación de compresión que uno duro, aunque el umbral duro tiene mejor calidad de imagen. La segunda hace referencia al momento de fuga (vanishing moment) de las wavelets, de modo que cuanto menor es el momento de fuga, mejores resultados de compresión se obtienen.

Ahora haremos lo mismo pero alterando el nivel de compresión. Los resultados se muestran en la tabla 4.

Tabla 4 – Resultados obtenidos alterando el nivel.

Alterando el nivel logramos una relación de compresión del 58,21%. A medida que aumenta el nivel la calidad de la imagen baja considerablemente por lo que una elección razonable sería elegir un nivel de compresión de 2 que ya nos reduce el tamaño de la imagen a la mitad.

Figura 15 – Resultado de la compresión rbio1.1 nivel 2.

Código Matlab

%{
---------------------------------------------------------------------------
EN: TS Signal Processor 1.0 developed by Garikoitz Martinez Moreno for the 
    Signal Processing course. Year 2015/2016
ES: TS Signal Processor desarrolado por Garikoitz Martinez Moreno para la
    asignatura de Procesado de Señales. Año 2015/2016
---------------------------------------------------------------------------
Compile: mcc -mve TSSP.m
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FORMULARIO PRINCIPAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = TSSP(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @TSSP_OpeningFcn, ...
                   'gui_OutputFcn',  @TSSP_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
function TSSP_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
movegui(hObject,'center');
%set(gcf, 'Position', get(0,'Screensize'));
%
if strcmp(get(hObject,'Visible'),'off')
    axes(handles.axes1);
    plot(rand(5));
    cla;
    set(gca,'fontsize',7)
    %No muestre valores ejes
    set(gca,'xtick',[]);
    set(gca,'ytick',[]);
end
%Load the TS Signals
try
    %Load images
    set(handles.popupmenu1,'String','');
    Files=dir('TS_Signals\*.dat');
    if numel(Files) > 0
      for k=1:length(Files)
         FileNames{k} = Files(k).name;
      end
      set(handles.popupmenu1,'String',FileNames);
    else
        set(handles.edit1, 'String', 'ERROR: TS_Signals\ folder is empty');
        set(handles.signaltools,'Enable','off');
    end
  set(handles.popupmenu_wavelet,'String',{'haar','db','sym','coif','bior','rbio'});
    set(handles.popupmenu_wavelet,'Value',6);
    set(handles.popupmenu_wavlevel,'String',{'1.1','1.3','1.5','2.2','2.4','2.6','2.8','3.1','3.3','3.5'});
  set(handles.popupmenu_thr,'String',{'lvd','gbl'});
    set(handles.popupmenu_thr,'Value',2);
    set(handles.popupmenu_thrsorh,'String',{'s','h'});
    set(handles.popupmenu_thrsorh,'Value',1);
    set(handles.popupmenu_noisetype,'String',{'gaussian','salt & pepper','speckle'})
    set(handles.popupmenu_noiseval,'String',{'0.01','0.02','0.05','0.08','0.1','0.5'});
    set(handles.popupmenu_noisetype,'Value',3);
    set(handles.popupmenu_noiseval,'Value',4);
    %
    a = get(handles.slider_level,'Value');
    set(handles.edit_level,'string',strcat('Level',{' '},num2str(round(a))));
    %
    set(handles.pushbutton6_noise, 'Enable', 'off');
    set(handles.pushbutton7_denoise, 'Enable', 'off');
    set(handles.pushbutton8_compress, 'Enable', 'off');
    set(handles.addnoise, 'Enable', 'off');
    set(handles.denoisesignal, 'Enable', 'off');
    set(handles.compression, 'Enable', 'off');
    %
    popupmenu_wavelet_Callback(hObject,eventdata,handles);
 catch me
    h = msgbox(me.message, 'Error','error');
end
function varargout = TSSP_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
%--------------------------------------------------------------------------
%FUNCIÓN VARIABLES GLOBALES
%--------------------------------------------------------------------------
function setGlobal(file,dir)
global filenameglobal
global dirglobal
filenameglobal = file;
dirglobal = dir;
function r = getGlobalf
global filenameglobal
r = filenameglobal;
function s = getGlobald
global dirglobal
s = dirglobal;
%--------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CREACIÓN DE OBJETOS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function popupmenu1_CreateFcn(hObject, eventdata, handles)
try
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
catch me
     h = msgbox({me.identifier me.message}, 'Error','error');
end
function edit1_CreateFcn(hObject, eventdata, handles)
try
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
catch me
    h = msgbox(me.message, 'Error','error');
end
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
try
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
catch me
    h = msgbox(me.message, 'Error','error');
end
function edit3_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit4_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit5_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit7_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit_info_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function popupmenu_wavelet_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit_level_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function popupmenu_wavlevel_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function slider_level_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider_level (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
    %
end
function edit_typethr_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function popupmenu_thr_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit9_Callback(hObject, eventdata, handles)
function edit9_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function popupmenu_thrsorh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function popupmenu_noiseval_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function popupmenu_noisetype_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit11_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function edit12_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%EVENTOS DE OBJETOS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CloseMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to CloseMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
selection = questdlg(['Close the program?'],...
                     ['Close '],...
                     'Yes','No','Yes');
if strcmp(selection,'No')
    return;
end
delete(handles.figure1)
catch me
    h = msgbox(me.message, 'Error','error');
end
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
popup_sel_index = get(handles.popupmenu1, 'Value');
popup_string = get(handles.popupmenu1, 'String');
selection = popup_string{popup_sel_index};
setGlobal(selection,'TS_Signals\');
filename = strcat('TS_Signals\', selection);
f = fopen(filename);
c = textscan(f,'%d %d %f','CollectOutput',1);
fclose(f);
Imag = accumarray(c{1},c{2}); % M (filas) y la segunda a  N (columnas)
Imag = imrotate(Imag,90);
axes(handles.axes1);
cla;
image(Imag,'CDataMapping','scaled')
set(gca,'xtick',[]);
set(gca,'ytick',[]);
msg = strcat('RAW Signal -',{' '}, selection);
set(handles.edit1, 'String', msg); 
%--------------------------------------------------------------------------
%Save the image
popup_sel_index = get(handles.popupmenu1, 'Value');
popup_string = get(handles.popupmenu1, 'String');
selection = popup_string{popup_sel_index};
[pathstr,name,ext] = fileparts(selection)
F = getframe(handles.axes1);
Image = frame2im(F);
sv = strcat('TS_saved_files\', name);
sv = strcat(sv,'.png');
imwrite(Image,sv,'png','TS_Signal processor',selection)
set(handles.pushbutton6_noise, 'Enable', 'on');
set(handles.pushbutton8_compress, 'Enable', 'on');
set(handles.addnoise, 'Enable', 'on');
set(handles.compression, 'Enable', 'on');

%--------------------------------------------------------------------------
catch me
     h = msgbox({me.identifier me.message}, 'Error','error');
end
%--------------------------------------------------------------------------
%MENÚ ABRIR SEÑAL
function opentssignal_Callback(hObject, eventdata, handles)
% hObject    handle to opentssignal (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
axes(handles.axes1);
cla;
[filename, pathname] = uigetfile('*.dat','Select the TS Signal');
if ~isequal(filename, 0)
    setGlobal(filename,pathname);
    f = fopen(strcat(pathname, filename));
    c = textscan(f,'%d %d %f','CollectOutput',1);
    fclose(f);
    Imag = accumarray(c{1},c{2}); % M (filas) y la segunda a  N (columnas)
    Imag = imrotate(Imag,90);
    axes(handles.axes1);
    cla;
    image(Imag,'CDataMapping','scaled')
    set(gca,'xtick',[]);
    set(gca,'ytick',[]);
    msg = strcat('RAW Signal -',{' '}, filename);
    set(handles.edit1, 'String', msg);
    %----------------------------------------------------------------------
    %Save the image
    [pathstr,name,ext] = fileparts(filename)
    F = getframe(handles.axes1);
    Image = frame2im(F);
    sv = strcat('TS_saved_files\', name);
    sv = strcat(sv,'.png');
    imwrite(Image,sv,'png','TS_Signal processor',filename)
    set(handles.pushbutton6_noise, 'Enable', 'on');
    set(handles.pushbutton8_compress, 'Enable', 'on');
    set(handles.addnoise, 'Enable', 'on');
    set(handles.compression, 'Enable', 'on');
    %----------------------------------------------------------------------
end
msgmnu = strcat('RAW Signal -',{' '}, filename);
set(handles.edit1, 'String', msgmnu); 
%-----------------------------------------------------
catch me
     h = msgbox({me.identifier me.message}, 'Error','error');
end
%--------------------------------------------------------------------------
%MENÚ ACERCA DE
function about_Callback(hObject, eventdata, handles)
% hObject    handle to about (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%h = msgbox('TS Signal Processor by GMM');
try
[cdata,map] = imread('TS_files\icon.png'); 
 h=msgbox({'TS Signal Processor 1.0','Developed by Garikoitz Martinez',...
            'UNED/UCM Signal Processing Course 2015/16'},'About','custom',cdata,summer);
catch me
   h = msgbox({me.identifier me.message}, 'Error','error');
end
%--------------------------------------------------------------------------
%MENÚ GUARDAR FIGURA
function savefig_Callback(hObject, eventdata, handles)
% hObject    handle to savefig (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
    selection=getGlobalf;
    pathstr=getGlobald;
    [pathstr,name,ext] = fileparts(strcat(pathstr, selection))
    F = getframe(handles.axes1);
    Image = frame2im(F);
    sv = strcat('TS_saved_files\', name);
    sv = strcat(sv,'.png');
    imwrite(Image,sv,'png','TS_Signal processor',selection)
    %--
catch me
     h = msgbox({me.identifier me.message}, 'Error','error');
end
%--------------------------------------------------------------------------
%MENÚ ELIMINAR RUIDO IMAGEN
function denoisesignal_Callback(hObject, eventdata, handles)
% hObject    handle to denoisesignal (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
    tipowaveid = get(handles.popupmenu_wavelet,'String');
    tipowave = tipowaveid{get(handles.popupmenu_wavelet,'Value')};
    levelid = get(handles.popupmenu_wavlevel,'String');
    level = levelid{get(handles.popupmenu_wavlevel,'Value')};
    switch get(handles.popupmenu_wavelet,'Value') 
      case 2 
          tipowave=strcat(tipowave,level);
      case 3 
          tipowave=strcat(tipowave,level);
      case 4 
          tipowave=strcat(tipowave,level);
      case 5 
          tipowave=strcat(tipowave,level);
      case 6 
          tipowave=strcat(tipowave,level);
      case 9 
          tipowave=strcat(tipowave,level);
      case 12 
          tipowave=strcat(tipowave,level);
      case 13 
          tipowave=strcat(tipowave,level);
      case 14 
          tipowave=strcat(tipowave,level);
      case 15 
          tipowave=strcat(tipowave,level);
      otherwise
    end 
selection=getGlobalf;
pathstr=getGlobald;
[pathstr,name,ext] = fileparts(strcat(pathstr, selection))
filename = strcat('TS_saved_files\', name);
filename = strcat(filename, '_denoised.png');
sv = strcat('TS_saved_files\', name);
sr = strcat(sv,'.png');
sv = strcat(sv,'_noisy.png');
if exist(sv, 'file')
    [SinRuido smap] = imread(sr);
    [ConRuido cmap] = imread(sv);
    %
    red = ConRuido(:, :, 1);
    green = ConRuido(:, :, 2);
    blue = ConRuido(:, :, 3);
    %
    %cmap = colormap(handles.axes1);
    n = round(get(handles.slider_level,'Value'));
    %[c l] = wavedec2(ConRuido,n,tipowave);
    thrsorhid = get(handles.popupmenu_thrsorh,'String');
    thrsorh = thrsorhid{get(handles.popupmenu_thrsorh,'Value')};
    switch get(handles.popupmenu_thrsorh,'Value') 
      case 1
          sorh = 's';
          threshold = 'Thresholding type = Soft ';
      case 2 
          sorh = 'h';
          threshold = 'Thresholding type = Hard ';
      otherwise
    end 
    thrtid = get(handles.popupmenu_thrsorh,'String');
    thrt = thrtid{get(handles.popupmenu_thrsorh,'Value')};
%DWT-----------------------------------------------------------------------
        figure('name','Noisy DWT RGB decomposition'); 
        set(gcf, 'Position', get(0,'Screensize'));
        redChannel = ConRuido(:, :, 1);
        greenChannel = ConRuido(:, :, 2);
        blueChannel = ConRuido(:, :, 3);
        subplot(4,4,1); imshow(ConRuido);title('Noisy image');
        subplot(4,4,2); imshow(redChannel);title('Red Channel');
        subplot(4,4,3); imshow(greenChannel);title('Green Channel');
        subplot(4,4,4); imshow(blueChannel);title('Blue Channel');
        %aR
        [CA,CH,CV,CD] = dwt2(redChannel,tipowave);
        subplot(4,4,5); imshow(CA/255);title('CA redChannel');
        subplot(4,4,6); imshow(CH/255);title('CH redChannel');
        subplot(4,4,7); imshow(CV/255);title('CV redChannel');
        subplot(4,4,8); imshow(CD/255);title('CD redChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(redChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(redChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),'h',thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aR=waverec2(c_n,l,tipowave);
        %aG
        [CA,CH,CV,CD] = dwt2(greenChannel,tipowave);
        subplot(4,4,9); imshow(CA/255);title('CA greenChannel');
        subplot(4,4,10); imshow(CH/255);title('CH greenChannel');
        subplot(4,4,11); imshow(CV/255);title('CV greenChannel');
        subplot(4,4,12); imshow(CD/255);title('CD greenChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(greenChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(greenChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),'h',thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aG=waverec2(c_n,l,tipowave);
        %aB
        [CA,CH,CV,CD] = dwt2(blueChannel,tipowave);
        subplot(4,4,13); imshow(CA/255);title('CA blueChannel');
        subplot(4,4,14); imshow(CH/255);title('CH blueChannel');
        subplot(4,4,15); imshow(CV/255);title('CV blueChannel');
        subplot(4,4,16); imshow(CD/255);title('CD blueChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(blueChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(blueChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),'h',thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aB=waverec2(c_n,l,tipowave);
        XDEN = cat(3, aR, aG, aB);
        XDEN=uint8(XDEN);
        tipowave=strcat(tipowave,' (DWT)');
%--------------------------------------------------------------------------
    Rs = SinRuido(:, :, 1);
    Gs = SinRuido(:, :, 2);
    Bs = SinRuido(:, :, 3);
    Rc = ConRuido(:, :, 1);
    Gc = ConRuido(:, :, 2);
    Bc = ConRuido(:, :, 3);
    Rd = XDEN(:, :, 1);
    Gd = XDEN(:, :, 2);
    Bd = XDEN(:, :, 3);
    figure('name','Denoised channels comparison'); 
    set(gcf, 'Position', get(0,'Screensize'));
    subplot(3,4,1); imshow(SinRuido);title('Original image');
    subplot(3,4,2); imshow(Rs);title('Red Channel');
    subplot(3,4,3); imshow(Gs);title('Green Channel');
    subplot(3,4,4); imshow(Bs);title('Blue Channel');
    subplot(3,4,5); imshow(SinRuido);title('Noisy image');
    subplot(3,4,6); imshow(Rc);title('Red Channel');
    subplot(3,4,7); imshow(Gc);title('Green Channel');
    subplot(3,4,8); imshow(Bc);title('Blue Channel');
    subplot(3,4,9); imshow(XDEN);title('Denoised image');
    subplot(3,4,10); imshow(Rd);title('Red Channel');
    subplot(3,4,11); imshow(Gd);title('Green Channel');
    subplot(3,4,12); imshow(Bd);title('Blue Channel');
    %----------------------------------------------------------------------
    imwrite(XDEN,filename);
    [imagen_final fmap]=imread(filename);
    msg = strcat('Signal -',{' '}, selection);
    msg = strcat(msg, ' denoised succesfully');
    set(handles.edit_info, 'String', msg); 
    %Calcs
    %http://es.mathworks.com/help/vision/ref/psnr.html
    %http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/VELDHUIZEN/node18.html
    signal = SinRuido(:);
    noisy  = ConRuido(:);
    denoised = imagen_final(:);
    %
    [rows columns ~] = size(ConRuido);
    mseRImage = (double(ConRuido(:,:,1)) - double(imagen_final(:,:,1))) .^ 2;
    mseGImage = (double(ConRuido(:,:,2)) - double(imagen_final(:,:,2))) .^ 2;
    mseBImage = (double(ConRuido(:,:,3)) - double(imagen_final(:,:,3))) .^ 2;
    mseR = sum(sum(mseRImage)) / (rows * columns);
    mseG = sum(sum(mseGImage)) / (rows * columns);
    mseB = sum(sum(mseBImage)) / (rows * columns);
    mse = (1/3)*(mseR + mseG + mseB);
    psnr = 10 * log10( 255^2 / mse);
    %----------------------------------------------------------------------
    [rows2 columns2 ~] = size(SinRuido);
    mseRImage2 = (double(SinRuido(:,:,1)) - double(ConRuido(:,:,1))) .^ 2;
    mseGImage2 = (double(SinRuido(:,:,2)) - double(ConRuido(:,:,2))) .^ 2;
    mseBImage2 = (double(SinRuido(:,:,3)) - double(ConRuido(:,:,3))) .^ 2;
    mseR2 = sum(sum(mseRImage2)) / (rows2 * columns2);
    mseG2 = sum(sum(mseGImage2)) / (rows2 * columns2);
    mseB2 = sum(sum(mseBImage2)) / (rows2 * columns2);
    mse2 = (1/3)*(mseR2 + mseG2 + mseB2);
    psnr2 = 10 * log10( 255^2 / mse2);
    %Show results
    figure('name','Results of noise removal'); 
    set(gcf, 'Position', get(0,'Screensize'));
    str=''; str2='';
    str = sprintf('Wavelet = %s\n%s \nPSNR = %.2fdB',strtrim(tipowave),strtrim(threshold),psnr);
    str2 = sprintf('PSNR = %.2fdB',psnr2);
    subplot(1,3,1);imshow(SinRuido,smap); title('Original Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
    subplot(1,3,2);imshow(ConRuido,cmap); title('Noisy Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]); xlabel(str2);
    subplot(1,3,3);imshow(imagen_final,fmap) ; title('Denoised Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]); xlabel(str);
    clear tipowave, clear snr, clear psnr, clear threshold, clear str;
else
  % File does not exist.
  warningMessage = sprintf('No existe la imagen con ruido:\n%s', sv);
  uiwait(msgbox(warningMessage));
end
catch me
    h = msgbox({me.identifier me.message}, 'Error denoising signal','Error');
end
%--------------------------------------------------------------------------
%MENÚ AÑADIR RUIDO IMAGEN
function addnoise_Callback(hObject, eventdata, handles)
try
%Añadir ruido
%https://en.wikipedia.org/wiki/Image_noise
%https://es.mathworks.com/matlabcentral/answers/76170-how-to-generate-background-noise-in-a-color-image
selection = getGlobalf;
pathstr=getGlobald;
[pathstr,name,ext] = fileparts(strcat(pathstr, selection))
filename = strcat('TS_saved_files\', name);
filename = strcat(filename, '.png');
Imagen_Inicial = imread(filename);
cmap = colormap(handles.axes1);%get colormap from fig
noiseid = get(handles.popupmenu_noisetype,'String');
noise = noiseid{get(handles.popupmenu_noisetype,'Value')};
nvalid = get(handles.popupmenu_noiseval,'String');
nval = nvalid{get(handles.popupmenu_noiseval,'Value')};
ConRuido = imnoise(Imagen_Inicial,noise,str2num(nval));
%--------------------------------------------------------------------------
msg = strcat('Noise applied to image -',{' '}, selection);
msg = strcat(msg,' succesfully');
set(handles.edit_info, 'String', msg); 
%--------------------------------------------------------------------------
%Show results
figure('name','Add noise to signal'); 
set(gcf, 'Position', get(0,'Screensize'));
subplot(1,2,1);imshow(Imagen_Inicial); title('Original Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
subplot(1,2,2);imshow(ConRuido) ; title('noisy Signal'); set(gca,'xtick',[]); set(gca,'ytick',[]);
%Save the image
[pathstr,name,ext] = fileparts(filename);
sv = strcat('TS_saved_files\', name);
sv = strcat(sv,'_noisy.png');
imwrite(ConRuido,cmap,sv);
set(handles.denoisesignal, 'Enable', 'on');
set(handles.pushbutton7_denoise, 'Enable', 'on');
%
catch me
     h = msgbox({me.identifier me.message}, 'Error adding noise','Error');
end
%--------------------------------------------------------------------------
%MENÚ COMPRIMIR IMAGEN
function compression_Callback(hObject, eventdata, handles)
% hObject    handle to compression (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%
try
    selection=getGlobalf;
    pathstr=getGlobald;
    [pathstr,name,ext] = fileparts(strcat(pathstr, selection));
    filename = strcat('TS_saved_files\', name);
    filename = strcat(filename, '.png');
    Imagen_Inicial = imread(filename);
    cmap = colormap(handles.axes1);
    tipowave='';
    level='';
    nivel='';
    tipowaveid = get(handles.popupmenu_wavelet,'String');
    tipowave = tipowaveid{get(handles.popupmenu_wavelet,'Value')};
    levelid = get(handles.popupmenu_wavlevel,'String');
    level = levelid{get(handles.popupmenu_wavlevel,'Value')};
    switch get(handles.popupmenu_wavelet,'Value') 
      case 2 
          tipowave=strcat(tipowave,level);
      case 3 
          tipowave=strcat(tipowave,level);
      case 4 
          tipowave=strcat(tipowave,level);
      case 5 
          tipowave=strcat(tipowave,level);
      case 6 
          tipowave=strcat(tipowave,level);
      case 9 
          tipowave=strcat(tipowave,level);
      case 12 
          tipowave=strcat(tipowave,level);
      case 13 
          tipowave=strcat(tipowave,level);
      case 14 
          tipowave=strcat(tipowave,level);
      case 15 
          tipowave=strcat(tipowave,level);
      otherwise
    end 
%--------------------------------------------------------------------------
% Compresión usando WDENCMP.
%http://read.pudn.com/downloads103/sourcecode/book/422976/wavelet/wavelet/wdencmp.m__.htm
%http://web.mit.edu/1.130/WebDocs/1.130/Software/Examples/example4.m
%http://es.mathworks.com/help/wavelet/examples/data-compression-using-2d-wavelet-analysis.html
%--------------------------------------------------------------------------
    n = round(get(handles.slider_level,'Value'));      % Decomposition Level
    [c l] = wavedec2(Imagen_Inicial,n,char(tipowave));  % Multilevel 2-D wavelet decomposition.
    thrsorhid = get(handles.popupmenu_thrsorh,'String');
    thrsorh = thrsorhid{get(handles.popupmenu_thrsorh,'Value')};
    switch get(handles.popupmenu_thrsorh,'Value') 
      case 1
          sorh = 's';
          threshold = 'Thresholding type = Soft ';
      case 2 
          sorh = 'h';
          threshold = 'Thresholding type = Hard ';
      otherwise
    end 
    thrtid = get(handles.popupmenu_thrsorh,'String');
    thrt = thrtid{get(handles.popupmenu_thrsorh,'Value')};
    switch get(handles.popupmenu_thr,'Value') 
      case 1
        opt = 'lvd';
        thr_h = [17 18];    % Horizontal thresholds.
        thr_d = [19 20];    % Diagonal thresholds.
        thr_v = [21 22];    % Vertical thresholds. 
        thr = [thr_h ; thr_d ; thr_v];
        redChannel = Imagen_Inicial(:, :, 1);
        greenChannel = Imagen_Inicial(:, :, 2);
        blueChannel = Imagen_Inicial(:, :, 3);
        [xdr,cxd2,lxd2,perf0,perfl2] = wdencmp(opt,redChannel,tipowave,n,thr,sorh);
        [xdg,cxd2,lxd2,perf0,perfl2] = wdencmp(opt,greenChannel,tipowave,n,thr,sorh);
        [xdb,cxd2,lxd2,perf0,perfl2] = wdencmp(opt,blueChannel,tipowave,n,thr,sorh);
        xd = cat(3, xdr, xdg, xdb);
        xd=uint8(xd);
        threshold=strcat(threshold,' and Level-dependant.');
      case 2 
          opt = 'gbl';
          thr=20;
          [xd,cxd,lxd,perf0,perfl2] = wdencmp(opt,c,l,tipowave,n,thr,sorh,1);
          xd=uint8(xd);
          threshold=strcat(threshold,' and Global.');
      otherwise
    end 
    %----------------------------------------------------------------------
    axes(handles.axes1);
    set(gca,'xtick',[]);
    set(gca,'ytick',[]);
    [pathstr,name,ext] = fileparts(filename);
    sv = strcat('TS_saved_files\', name);
    sv = strcat(sv,'_');
    sv = strcat(sv,tipowave);
    sv = strcat(sv,'_Level');
    sv = strcat(sv,num2str(n));
    sv = strcat(sv,'_thr_');
    sv = strcat(sv,sorh);
    sv = strcat(sv,'_');
    sv = strcat(sv,opt);
    sv = strcat(sv,'_compressed.png');
    sv = regexprep(sv,{'\.','png'},{'','.png'})
    imwrite(xd,cmap,sv);
    msg = strcat('Signal -',{' '}, selection);
    set(handles.edit_info, 'String', msg); 
    %Show results
    figure('name','Results of image compression'); 
    set(gcf, 'Position', get(0,'Screensize'));
    Imagen_final = imread(sv);
    cla;
    subplot(2,1,1);image(Imagen_Inicial,'CDataMapping','scaled'); title('Original Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
    subplot(2,1,2);image(Imagen_final,'CDataMapping','scaled'); title('Compressed Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
    xlabel(sprintf('Wavelet = %s || Compression level = %d \n %s \nCompression score = %2.1f%% || Null coefficients = %2.1f%%',tipowave,n,threshold,perfl2,perf0))
catch me
     h = msgbox({me.identifier me.message}, 'Error in compress button','error');
end
%--------------------------------------------------------------------------
%BOTÓN AÑADIR RUIDO IMAGEN
function pushbutton6_noise_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6_noise (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
try
%Añadir ruido
%https://en.wikipedia.org/wiki/Image_noise
%https://es.mathworks.com/matlabcentral/answers/76170-how-to-generate-background-noise-in-a-color-image
selection = getGlobalf;
pathstr=getGlobald;
[pathstr,name,ext] = fileparts(strcat(pathstr, selection))
filename = strcat('TS_saved_files\', name);
filename = strcat(filename, '.png');
Imagen_Inicial = imread(filename);
cmap = colormap(handles.axes1);%get colormap from fig
noiseid = get(handles.popupmenu_noisetype,'String');
noise = noiseid{get(handles.popupmenu_noisetype,'Value')};
nvalid = get(handles.popupmenu_noiseval,'String');
nval = nvalid{get(handles.popupmenu_noiseval,'Value')};
ConRuido = imnoise(Imagen_Inicial,noise,str2num(nval));
%--------------------------------------------------------------------------
msg = strcat('Noise applied to image -',{' '}, selection);
msg = strcat(msg,' succesfully');
set(handles.edit_info, 'String', msg); 
%--------------------------------------------------------------------------
%Show results
figure('name','Add noise to signal'); 
set(gcf, 'Position', get(0,'Screensize'));
subplot(1,2,1);imshow(Imagen_Inicial); title('Original Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
subplot(1,2,2);imshow(ConRuido) ; title('noisy Signal'); set(gca,'xtick',[]); set(gca,'ytick',[]);
%Save the image
[pathstr,name,ext] = fileparts(filename);
%F = getframe(handles.axes1);
%Image = frame2im(F);
sv = strcat('TS_saved_files\', name);
sv = strcat(sv,'_noisy.png');
imwrite(ConRuido,cmap,sv);
set(handles.pushbutton7_denoise, 'Enable', 'on');
set(handles.denoisesignal, 'Enable', 'on');
catch me
     h = msgbox({me.identifier me.message}, 'Error adding noise','Error');
end
%--------------------------------------------------------------------------
%BOTÓN ELIMINAR RUIDO IMAGEN
function pushbutton7_denoise_Callback(hObject, eventdata, handles)
try
    tipowaveid = get(handles.popupmenu_wavelet,'String');
    tipowave = tipowaveid{get(handles.popupmenu_wavelet,'Value')};
    levelid = get(handles.popupmenu_wavlevel,'String');
    level = levelid{get(handles.popupmenu_wavlevel,'Value')};
    switch get(handles.popupmenu_wavelet,'Value') 
      case 2 
          tipowave=strcat(tipowave,level);
      case 3 
          tipowave=strcat(tipowave,level);
      case 4 
          tipowave=strcat(tipowave,level);
      case 5 
          tipowave=strcat(tipowave,level);
      case 6 
          tipowave=strcat(tipowave,level);
      case 9 
          tipowave=strcat(tipowave,level);
      case 12 
          tipowave=strcat(tipowave,level);
      case 13 
          tipowave=strcat(tipowave,level);
      case 14 
          tipowave=strcat(tipowave,level);
      case 15 
          tipowave=strcat(tipowave,level);
      otherwise
    end 
selection=getGlobalf;
pathstr=getGlobald;
[pathstr,name,ext] = fileparts(strcat(pathstr, selection))
filename = strcat('TS_saved_files\', name);
filename = strcat(filename, '_denoised.png');
sv = strcat('TS_saved_files\', name);
sr = strcat(sv,'.png');
sv = strcat(sv,'_noisy.png');
if exist(sv, 'file')
    [SinRuido smap] = imread(sr);
    [ConRuido cmap] = imread(sv);
    %
    red = ConRuido(:, :, 1);
    green = ConRuido(:, :, 2);
    blue = ConRuido(:, :, 3);
    %
    n = round(get(handles.slider_level,'Value'));
    thrsorhid = get(handles.popupmenu_thrsorh,'String');
    thrsorh = thrsorhid{get(handles.popupmenu_thrsorh,'Value')};
    switch get(handles.popupmenu_thrsorh,'Value') 
      case 1
          sorh = 's';
          threshold = 'Thresholding type = Soft ';
      case 2 
          sorh = 'h';
          threshold = 'Thresholding type = Hard ';
      otherwise
    end 
    thrtid = get(handles.popupmenu_thrsorh,'String');
    thrt = thrtid{get(handles.popupmenu_thrsorh,'Value')};
%DWT-----------------------------------------------------------------------
        figure('name','Noisy DWT RGB decomposition'); 
        set(gcf, 'Position', get(0,'Screensize'));
        redChannel = ConRuido(:, :, 1);
        greenChannel = ConRuido(:, :, 2);
        blueChannel = ConRuido(:, :, 3);
        subplot(4,4,1); imshow(ConRuido);title('Noisy image');
        subplot(4,4,2); imshow(redChannel);title('Red Channel');
        subplot(4,4,3); imshow(greenChannel);title('Green Channel');
        subplot(4,4,4); imshow(blueChannel);title('Blue Channel');
        %aR
        [CA,CH,CV,CD] = dwt2(redChannel,tipowave);
        subplot(4,4,5); imshow(CA/255);title('CA redChannel');
        subplot(4,4,6); imshow(CH/255);title('CH redChannel');
        subplot(4,4,7); imshow(CV/255);title('CV redChannel');
        subplot(4,4,8); imshow(CD/255);title('CD redChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(redChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(redChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),sorh,thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aR=waverec2(c_n,l,tipowave);
        %aG
        [CA,CH,CV,CD] = dwt2(greenChannel,tipowave);
        subplot(4,4,9); imshow(CA/255);title('CA greenChannel');
        subplot(4,4,10); imshow(CH/255);title('CH greenChannel');
        subplot(4,4,11); imshow(CV/255);title('CV greenChannel');
        subplot(4,4,12); imshow(CD/255);title('CD greenChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(greenChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(greenChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),sorh,thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aG=waverec2(c_n,l,tipowave);
        %aB
        [CA,CH,CV,CD] = dwt2(blueChannel,tipowave);
        subplot(4,4,13); imshow(CA/255);title('CA blueChannel');
        subplot(4,4,14); imshow(CH/255);title('CH blueChannel');
        subplot(4,4,15); imshow(CV/255);title('CV blueChannel');
        subplot(4,4,16); imshow(CD/255);title('CD blueChannel');
        noiselev = median(abs(CD(:)))/0.6745;
        thresh = sqrt(2*log(numel(blueChannel)))*noiselev;
        clear CA, clear CH, clear CV, clear CD;
        [c,l] = wavedec2(blueChannel,n,tipowave);
        Y = wthresh(c(prod(l(1,:))+1:end),sorh,thresh);
        c_n = c;
        c_n(prod(l(1,:))+1:end) = Y;
        aB=waverec2(c_n,l,tipowave);
        XDEN = cat(3, aR, aG, aB);
        XDEN=uint8(XDEN);
        tipowave=strcat(tipowave,' (DWT)');
%--------------------------------------------------------------------------
    Rs = SinRuido(:, :, 1);
    Gs = SinRuido(:, :, 2);
    Bs = SinRuido(:, :, 3);
    Rc = ConRuido(:, :, 1);
    Gc = ConRuido(:, :, 2);
    Bc = ConRuido(:, :, 3);
    Rd = XDEN(:, :, 1);
    Gd = XDEN(:, :, 2);
    Bd = XDEN(:, :, 3);
    figure('name','Denoised channels comparison'); 
    set(gcf, 'Position', get(0,'Screensize'));
    subplot(3,4,1); imshow(SinRuido);title('Original image');
    subplot(3,4,2); imshow(Rs);title('Red Channel');
    subplot(3,4,3); imshow(Gs);title('Green Channel');
    subplot(3,4,4); imshow(Bs);title('Blue Channel');
    subplot(3,4,5); imshow(ConRuido);title('Noisy image');
    subplot(3,4,6); imshow(Rc);title('Red Channel');
    subplot(3,4,7); imshow(Gc);title('Green Channel');
    subplot(3,4,8); imshow(Bc);title('Blue Channel');
    subplot(3,4,9); imshow(XDEN);title('Denoised image');
    subplot(3,4,10); imshow(Rd);title('Red Channel');
    subplot(3,4,11); imshow(Gd);title('Green Channel');
    subplot(3,4,12); imshow(Bd);title('Blue Channel');
    %----------------------------------------------------------------------
    imwrite(XDEN,filename);
    [imagen_final fmap]=imread(filename);
    msg = strcat('Signal -',{' '}, selection);
    msg = strcat(msg, ' denoised succesfully');
    set(handles.edit_info, 'String', msg); 
    %Calcs
    %http://es.mathworks.com/help/vision/ref/psnr.html
    %http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/VELDHUIZEN/node18.html
    signal = SinRuido(:);
    noisy  = ConRuido(:);
    denoised = imagen_final(:);
    %CALCULA ERRORES
    [rows columns ~] = size(ConRuido);
    mseRImage = (double(ConRuido(:,:,1)) - double(imagen_final(:,:,1))) .^ 2;
    mseGImage = (double(ConRuido(:,:,2)) - double(imagen_final(:,:,2))) .^ 2;
    mseBImage = (double(ConRuido(:,:,3)) - double(imagen_final(:,:,3))) .^ 2;
    mseR = sum(sum(mseRImage)) / (rows * columns);
    mseG = sum(sum(mseGImage)) / (rows * columns);
    mseB = sum(sum(mseBImage)) / (rows * columns);
    mse = (1/3)*(mseR + mseG + mseB);
    psnr = 10 * log10( 255^2 / mse);
    %----------------------------------------------------------------------
    [rows2 columns2 ~] = size(SinRuido);
    mseRImage2 = (double(SinRuido(:,:,1)) - double(ConRuido(:,:,1))) .^ 2;
    mseGImage2 = (double(SinRuido(:,:,2)) - double(ConRuido(:,:,2))) .^ 2;
    mseBImage2 = (double(SinRuido(:,:,3)) - double(ConRuido(:,:,3))) .^ 2;
    mseR2 = sum(sum(mseRImage2)) / (rows2 * columns2);
    mseG2 = sum(sum(mseGImage2)) / (rows2 * columns2);
    mseB2 = sum(sum(mseBImage2)) / (rows2 * columns2);
    mse2 = (1/3)*(mseR2 + mseG2 + mseB2);
    psnr2 = 10 * log10( 255^2 / mse2);
    %Show results
    figure('name','Results of noise removal'); 
    set(gcf, 'Position', get(0,'Screensize'));
    str=''; str2='';
    str = sprintf('Wavelet = %s\n%s \nPSNR = %.2fdB || MSE = %.2f',strtrim(tipowave),strtrim(threshold),psnr,mse);
    str2 = sprintf('PSNR = %.2fdB || MSE = %.2f',psnr2,mse2);
    subplot(1,3,1);imshow(SinRuido,smap); title('Original Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
    subplot(1,3,2);imshow(ConRuido,cmap); title('Noisy Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]); xlabel(str2);
    subplot(1,3,3);imshow(imagen_final,fmap) ; title('Denoised Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]); xlabel(str);
    clear tipowave, clear snr, clear psnr, clear threshold, clear str;
else
  % File does not exist.
  warningMessage = sprintf('No existe la imagen con ruido:\n%s', sv);
  uiwait(msgbox(warningMessage));
end
catch me
    h = msgbox({me.identifier me.message}, 'Error denoising signal','Error');
end
%--------------------------------------------------------------------------
%BOTÓN COMPRIMIR IMAGEN
function pushbutton8_compress_Callback(hObject, eventdata, handles)
try
    selection=getGlobalf;
    pathstr=getGlobald;
    [pathstr,name,ext] = fileparts(strcat(pathstr, selection));
    filename = strcat('TS_saved_files\', name);
    filename = strcat(filename, '.png');
    Imagen_Inicial = imread(filename);
    cmap = colormap(handles.axes1);
    tipowave='';
    level='';
    nivel='';
    tipowaveid = get(handles.popupmenu_wavelet,'String');
    tipowave = tipowaveid{get(handles.popupmenu_wavelet,'Value')};
    levelid = get(handles.popupmenu_wavlevel,'String');
    level = levelid{get(handles.popupmenu_wavlevel,'Value')};
    switch get(handles.popupmenu_wavelet,'Value') 
      case 2 
          tipowave=strcat(tipowave,level);
      case 3 
          tipowave=strcat(tipowave,level);
      case 4 
          tipowave=strcat(tipowave,level);
      case 5 
          tipowave=strcat(tipowave,level);
      case 6 
          tipowave=strcat(tipowave,level);
      case 9 
          tipowave=strcat(tipowave,level);
      case 12 
          tipowave=strcat(tipowave,level);
      case 13 
          tipowave=strcat(tipowave,level);
      case 14 
          tipowave=strcat(tipowave,level);
      case 15 
          tipowave=strcat(tipowave,level);
      otherwise
    end 
%--------------------------------------------------------------------------
% Compresión usando WDENCMP.
%http://read.pudn.com/downloads103/sourcecode/book/422976/wavelet/wavelet/wdencmp.m__.htm
%http://web.mit.edu/1.130/WebDocs/1.130/Software/Examples/example4.m
%http://es.mathworks.com/help/wavelet/examples/data-compression-using-2d-wavelet-analysis.html
%--------------------------------------------------------------------------
    n = round(get(handles.slider_level,'Value'));      % Decomposition Level
    thrsorhid = get(handles.popupmenu_thrsorh,'String');
    thrsorh = thrsorhid{get(handles.popupmenu_thrsorh,'Value')};
    switch get(handles.popupmenu_thrsorh,'Value') 
      case 1
          sorh = 's';
          threshold = 'Thresholding type = Soft ';
      case 2 
          sorh = 'h';
          threshold = 'Thresholding type = Hard ';
      otherwise
    end 
    thrtid = get(handles.popupmenu_thrsorh,'String');
    thrt = thrtid{get(handles.popupmenu_thrsorh,'Value')};
    switch get(handles.popupmenu_thr,'Value') 
      case 1
        opt = 'lvd';
        thr_h = [17 18];    % Horizontal thresholds.
        thr_d = [19 20];    % Diagonal thresholds.
        thr_v = [21 22];    % Vertical thresholds. 
        thr = [thr_h ; thr_d ; thr_v];
        redChannel = Imagen_Inicial(:, :, 1);
        greenChannel = Imagen_Inicial(:, :, 2);
        blueChannel = Imagen_Inicial(:, :, 3);
        %Al utilizar thr como matrix n=2, de lo contrario dará error.
        [xdr,cxd2r,lxd2r,perf0r,perfl2r] = wdencmp(opt,redChannel,tipowave,2,thr,sorh);
        [xdg,cxd2g,lxd2g,perf0g,perfl2g] = wdencmp(opt,greenChannel,tipowave,2,thr,sorh);
        [xdb,cxd2b,lxd2b,perf0b,perfl2b] = wdencmp(opt,blueChannel,tipowave,2,thr,sorh);
        xd = cat(3, xdr, xdg, xdb);
        perf0 = (1/3)*(perf0r+perf0g+perf0b);
        perfl2 = (1/3)*(perfl2r+perfl2g+perfl2b);
        xd=uint8(xd);
        threshold=strcat(threshold,' and Level-dependant.');
      case 2 
        opt = 'gbl';
        [c s] = wavedec2(Imagen_Inicial,n,char(tipowave));
        thr=20;
        [xd,cxd,lxd,perf0,perfl2] = wdencmp(opt,c,s,tipowave,n,thr,sorh,1);
        xd=uint8(xd);
        threshold=strcat(threshold,' and Global.');
      otherwise
    end 
    %----------------------------------------------------------------------
    axes(handles.axes1);
    set(gca,'xtick',[]);
    set(gca,'ytick',[]);
    [pathstr,name,ext] = fileparts(filename);
    sv = strcat('TS_saved_files\', name);
    sv = strcat(sv,'_');
    sv = strcat(sv,tipowave);
    sv = strcat(sv,'_Level');
    sv = strcat(sv,num2str(n));
    sv = strcat(sv,'_thr_');
    sv = strcat(sv,sorh);
    sv = strcat(sv,'_');
    sv = strcat(sv,opt);
    sv = strcat(sv,'_compressed.png');
    sv = regexprep(sv,{'\.','png'},{'','.png'})
    imwrite(xd,cmap,sv);
    msg = strcat('Signal -',{' '}, selection);
    set(handles.edit_info, 'String', msg); 
    %Show results
    figure('name','Results of image compression'); 
    set(gcf, 'Position', get(0,'Screensize'));
    Imagen_final = imread(sv);
    cla;
    subplot(2,1,1);image(Imagen_Inicial,'CDataMapping','scaled'); title('Original Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
    subplot(2,1,2);image(Imagen_final,'CDataMapping','scaled'); title('Compressed Signal'); set(gca,'xtick',[]);set(gca,'ytick',[]);
    xlabel(sprintf('Wavelet = %s || Compression level = %d \n %s \nCompression score = %2.1f%% || Null coefficients = %2.1f%%',tipowave,n,threshold,perfl2,perf0))
catch me
     h = msgbox({me.identifier me.message}, 'Error in compress button','error');
end
function edit_info_Callback(hObject, eventdata, handles)
function popupmenu_wavelet_Callback(hObject, eventdata, handles)
try
    switch get(handles.popupmenu_wavelet,'Value')   
  case 1
      set(handles.popupmenu_wavlevel, 'Enable', 'off');
  case 2
      set(handles.popupmenu_wavlevel, 'Enable', 'on');
      set(handles.popupmenu_wavlevel,'String',{'1','2','3','4','5','6','7','8','9','10'});
  case 3
      set(handles.popupmenu_wavlevel, 'Enable', 'on');
      set(handles.popupmenu_wavlevel,'String',{'2','3','4','5','6','7','8'});
  case 4
      set(handles.popupmenu_wavlevel, 'Enable', 'on');
      set(handles.popupmenu_wavlevel,'String',{'1','2','3','4','5'});
  case 5
      set(handles.popupmenu_wavlevel, 'Enable', 'on');
      set(handles.popupmenu_wavlevel,'String',{'1.1','1.3','1.5','2.2','2.4','2.6','2.8','3.1','3.3','3.5'});
  case 6
      set(handles.popupmenu_wavlevel, 'Enable', 'on');
      set(handles.popupmenu_wavlevel,'String',{'1.1','1.3','1.5','2.2','2.4','2.6','2.8','3.1','3.3','3.5'});
%   case 7
%       set(handles.popupmenu_wavlevel, 'Enable', 'off');
%   case 8
%       set(handles.popupmenu_wavlevel, 'Enable', 'off');
%   case 9
%       set(handles.popupmenu_wavlevel, 'Enable', 'on');
%       set(handles.popupmenu_wavlevel,'String',{'1','2','3','4','5','6','7','8'});
%   case 10
%       set(handles.popupmenu_wavlevel, 'Enable', 'off');
%   case 11
%       set(handles.popupmenu_wavlevel, 'Enable', 'off');
%   case 12
%       set(handles.popupmenu_wavlevel, 'Enable', 'on');
%       set(handles.popupmenu_wavlevel,'String',{'1','2','3','4','5'});
%   case 13
%       set(handles.popupmenu_wavlevel, 'Enable', 'on');
%       set(handles.popupmenu_wavlevel,'String',{'1-1.5','1-1','1-0.5'});
%   case 14
%       set(handles.popupmenu_wavlevel, 'Enable', 'on');
%       set(handles.popupmenu_wavlevel,'String',{'1-1-1.5','1-1-1','1-1-0.5'});
%   case 15
%       set(handles.popupmenu_wavlevel, 'Enable', 'on');
%       set(handles.popupmenu_wavlevel,'String',{'1-1.5','1-1','1-0.5'});
  otherwise
end 
   catch me
    h = msgbox(me.message, 'Error','error');
end
function edit7_Callback(hObject, eventdata, handles)
function edit_level_Callback(hObject, eventdata, handles)
function popupmenu_wavlevel_Callback(hObject, eventdata, handles)
function similwave_Callback(hObject, eventdata, handles)
% hObject    handle to similwave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%http://fudaqs2.ciemat.es/TJ2WEB/SimilWave.html
url = 'http://fudaqs2.ciemat.es/TJ2WEB/SimilWave.html';
web(url,'-browser')
function ciematrepo_Callback(hObject, eventdata, handles)
% hObject    handle to ciematrepo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%http://fudaqs2.ciemat.es/TJ2WEB/
url = 'http://fudaqs2.ciemat.es/TJ2WEB/';
web(url,'-browser')
function usingwavtoolbox_Callback(hObject, eventdata, handles)
% hObject    handle to usingwavtoolbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%http://nf.nci.org.au/facilities/software/Matlab/toolbox/wavelet/ch01_i24.html
url = 'http://nf.nci.org.au/facilities/software/Matlab/toolbox/wavelet/ch01_i24.html';
web(url,'-browser')
function mwdocu_Callback(hObject, eventdata, handles)
% hObject    handle to mwdocu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%http://es.mathworks.com/help/matlab/index.html
url = 'http://es.mathworks.com/help/matlab/index.html';
web(url,'-browser')
function slider_level_Callback(hObject, eventdata, handles)
% hObject    handle to slider_level (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

a = get(handles.slider_level,'Value');
set(handles.edit_level,'string',strcat('Level',{' '},num2str(round(a))));
function edit_typethr_Callback(hObject, eventdata, handles)
function popupmenu_thr_Callback(hObject, eventdata, handles)
function popupmenu_thrsorh_Callback(hObject, eventdata, handles)
function popupmenu_noisetype_Callback(hObject, eventdata, handles)
function edit11_Callback(hObject, eventdata, handles)
function edit12_Callback(hObject, eventdata, handles)
function popupmenu_noiseval_Callback(hObject, eventdata, handles)

Archivos del proyecto

Archivos fuente (Sin compilar. Necesario Matlab >2016)
Señales TJ-II

Bibliografía

[1] C. Srisailam, P. Sharma, S. Suhane. Color Image Denoising Using Wavelet soft Thresholding. International Journal of Emerging Technology and Advanced Engineering. IJETAE Journal. [en línea]. Jabalpur 2014. Disponible en: https://goo.gl/AtT98e

[2] Dimitris G. Manilakis, Vinay K. Ingle. Applied Digital Signal Processing. London, UK: Cambridge University Press, 2013. ISBN: 978-0-521-11002-0

[3] G. P. Guerrero, H. S. Escobar. Compresión de imágenes usando wavelets. Trabajo de investigación. [en línea]. Medellín, 2007. Disponible en: https://goo.gl/en5vND

[4] L. Maliki, S. Dormido-Canto, J. Vega. Clasificador de Imágenes del Diagnóstico Thomson Scattering del TJ II Basado en Template Matching. [en línea]. España, 2008. Disponible en: https://goo.gl/8Pq97n

[5] M. Misiti, Y. Misiti, G. Oppenheim, JM. Poggi. Wavelet Toolbox user’s guide. R2016a. [en línea]. Natick (MA), EE.UU, 2016. Disponible en: https://goo.gl/T82Efm

Información básica sobre protección de datos
Responsable	Garikoitz Martínez Moreno +info...
Finalidad	Gestionar y moderar tus comentarios. +info...
Legitimación	Consentimiento del interesado. +info...
Destinatarios	Automattic Inc., EEUU para filtrar el spam. +info...
Derechos	Acceder, rectificar y cancelar los datos, así como otros derechos. +info...
Información adicional	Puedes consultar la información adicional y detallada sobre protección de datos en nuestra página de política de privacidad.