function [objetooverlapping,restaurada,lambdasoptimas]=regu_perceptual(Imobservada,b,lambdas,N,M,H,k1,k2,umbral,criterio_lambda,pinta); Dada una imagen observada que ha pasado por un proceso de degradación representado por K y a la que se ha añadido ruido, se intenta regularizar usando el operador de penalización L y un parámetro de regularización lambda. N es el tamaño del bloque y M es el tamaño del bloque más un borde, esto se hace para que los efectos de borde no se vean en la solución. pinta es un valor 0/1 para indicar si se muestra una figura de la evolución de la solución. EJEMPLO: N=16; M=32; [H,k1,k2]=constantes2(fs,M); Nlambdas=70; lambdas=logspace(-2,2,Nlambdas); figure,[IrestPinv,lamboptPinv]=regu_perceptual(BlurredNoisy,PSF,lambdas ,N,M,H,k1,k2,umbral,2,1);
0001 % function [objetooverlapping,restaurada,lambdasoptimas]=regu_perceptual(Imobservada,b,lambdas,N,M,H,k1,k2,umbral,criterio_lambda,pinta); 0002 % 0003 % Dada una imagen observada que ha pasado por un proceso de degradación representado 0004 % por K y a la que se ha añadido ruido, se intenta regularizar usando el operador 0005 % de penalización L y un parámetro de regularización lambda. N es el tamaño del bloque 0006 % y M es el tamaño del bloque más un borde, esto se hace para que los efectos de borde 0007 % no se vean en la solución. 0008 % pinta es un valor 0/1 para indicar si se muestra una figura de la evolución de la solución. 0009 % 0010 % 0011 % EJEMPLO: 0012 % 0013 % N=16; 0014 % M=32; 0015 % 0016 % [H,k1,k2]=constantes2(fs,M); 0017 % 0018 % Nlambdas=70; 0019 % 0020 % lambdas=logspace(-2,2,Nlambdas); 0021 % 0022 % figure,[IrestPinv,lamboptPinv]=regu_perceptual(BlurredNoisy,PSF,lambdas ,N,M,H,k1,k2,umbral,2,1); 0023 0024 0025 0026 0027 0028 function [restaurada,lambdamedio]=regu_perceptual(Imobservada,b,lambdas,N,M,H,k1,k2,umbral,criterio_lambda,pinta); 0029 0030 0031 [fil,col]=size(Imobservada); 0032 0033 0034 nver = fil/N; 0035 nhor = col/N; 0036 0037 objetoparaoverlapping = rand(nver*M,nhor*M); 0038 0039 0040 restaurada=rand([fil,col]); 0041 0042 a=(M-N)/2; 0043 0044 Im2=ampliaconborde(Imobservada,a); 0045 0046 [fil2,col2]=size(Im2); 0047 0048 0049 nlamb=length(lambdas); 0050 0051 0052 B=freqz2(b,M,M); 0053 0054 B=ifftshift(B); 0055 0056 BB = abs(B).^2; 0057 Bconj = conj(B); 0058 0059 energB = sum(sum(abs(B).^2)); 0060 0061 0062 contadorfilas=0; 0063 0064 valorlambda=0; 0065 0066 lam=[]; 0067 0068 %r=zeros(M); 0069 % 0070 %contador = 1; 0071 %for i=N/2:N:fil 0072 % for j=N/2:N:col 0073 % r = r+ respuestafourier(Im2(i+a-M/2+1:i+a+M/2,j+a-M/2+1:j+a+M/2),M,H,k1,k2); 0074 % contador = contador +1; 0075 % end 0076 %end 0077 0078 0079 %% Estimación de la respuesta media en varios bloques 0080 0081 0082 r=zeros(M); 0083 0084 0085 Nbloques = 30; 0086 0087 filas = round(N/2 + rand(Nbloques,1)*(fil-N)); 0088 columnas = round(N/2 + rand(Nbloques,1)*(col-N)); 0089 0090 lamb = 0; 0091 0092 %R = []; 0093 0094 for i=1:Nbloques 0095 0096 observado=Im2(filas(i)+a-M/2+1:filas(i)+a+M/2,columnas(i)+a-M/2+1:columnas(i)+a+M/2); 0097 0098 r = r+ respuestafourier(observado,M,H,k1,k2); 0099 0100 %R = [R r]; 0101 end 0102 0103 0104 respuestamedia = r/Nbloques; 0105 0106 maximo = max(max(respuestamedia)); 0107 %media = mean(mean(R)); 0108 0109 0110 %% Fin estimación de la respuesta media 0111 0112 0113 0114 0115 %%% Estimación del parametro de regularizacion 0116 % 0117 %Nbloques = 25; 0118 % 0119 %filas = round(N/2 + rand(Nbloques,1)*(fil-N)); 0120 %columnas = round(N/2 + rand(Nbloques,1)*(col-N)); 0121 % 0122 %lamb = 0; 0123 % 0124 %for i=1:Nbloques 0125 % 0126 % observado=Im2(filas(i)+a-M/2+1:filas(i)+a+M/2,columnas(i)+a-M/2+1:columnas(i)+a+M/2); 0127 % 0128 % F=costeinvrespfour(observado,M,H,k1,k2,maximo/umbral); 0129 % 0130 % [sol,lambda]=restaura_bloque(observado,B,F,M,N,lambdas,2); 0131 % 0132 % lamb = lamb + lambda; 0133 % 0134 %end 0135 % 0136 %lamb = lamb / Nbloques 0137 % 0138 %%% Fin estimación del parametro de regularizacion 0139 % 0140 0141 0142 0143 for i=N/2:N:fil 0144 %fprintf('\nFila %d : ',i) 0145 contadorfilas = contadorfilas+1; 0146 contadorcolumnas = 0; 0147 0148 for j=N/2:N:col 0149 %fprintf('%d, ',j) 0150 contadorcolumnas = contadorcolumnas + 1; 0151 0152 observado=Im2(i+a-M/2+1:i+a+M/2,j+a-M/2+1:j+a+M/2); 0153 0154 %umb = umbral*mean2(Im2(i-N/2+1:i+N/2,j-N/2+1:j+N/2))/256 0155 0156 % [dummy,varianzaruido] = wiener2(observado,[2 2]); 0157 0158 % if (mean2(observado)<128) 0159 % umbralito = umbral/2; 0160 % else 0161 % umbralito = umbral; 0162 % end 0163 0164 F=costeinvrespfour(observado,M,H,k1,k2,maximo/umbral); 0165 0166 [sol,lambd]=restaura_bloque(observado,B,F,M,N,lambdas,criterio_lambda); 0167 0168 lam=[lam lambd]; 0169 0170 restaurada(i-N/2+1:i+N/2,j-N/2+1:j+N/2)=sol(a+1:a+N,a+1:a+N); 0171 0172 0173 if (pinta == true) 0174 figure(500),imshow(restaurada,[0 256],'notruesize'),pause(0.2); 0175 end 0176 0177 objetoparaoverlapping((contadorfilas-1)*M+1:contadorfilas*M,(contadorcolumnas-1)*M+1:contadorcolumnas*M) = sol; 0178 0179 0180 end 0181 end 0182 0183 lambdamedio = mean(lam); 0184 0185 restaurada = do_overlap(objetoparaoverlapping,M,N,8); 0186 0187 0188 0189 0190 0191 0192 0193 0194 0195 0196