2015-06-06 13:31:00 weixin_33853827 阅读数 95

局部标准差在图像处理邻域具有广泛的应用,但是直接计算非常耗时,本文利用积分图像对局部标准差的计算进行加速。

局部标准差:

标准差定义如下(采用统计学中的定义,分母为):


其中

为了计算图像的局部标准差,首先设定局部区域的大小为 ,则局部区域的像素点个数

对标准差的公式进行化简:


,故:


可以看出,局部标准差计算中需要对图像进行局部求和操作,即

我们可以先分别计算出图像的积分图像,这样就能在常量时间计算出上述的局部和。

时间复杂度:

图像中共个像素点,每个局部区域共有个像素点,直接计算局部标准差的时间复杂度 ,利用积分图像计算局部标准差的时间复杂度为。(计算积分图像的时间为 、计算一个像素点的局部标准差需要常量时间,共个像素点,故最后的计算复杂度仍为 


Matlab程序:

Localstd.m:直接计算图像的局部标准差

function Std=Localstd(Img,d)
[m,n]=size(Img);
Var=zeros(m,n);
Img = padarray(Img,[d d],'symmetric');%边界填充
N=(2*d+1)^2;
for i=1:m
    for j=1:n
        i1=i+d;
        j1=j+d;
        B=Img(i1-d:i1+d,j1-d:j1+d);%取局部区域
        Var(i,j)=sum(sum((B-mean(B(:))).^2))/(N-1);%标准差公式
    end
end
Std=sqrt(Var);

fastLocalstd.m:利用积分图像进行快速计算

function Std=fastLocalstd(Img,d)
[m,n]=size(Img);
Img = padarray(Img,[d+1 d+1],'symmetric');%边界填充
Img2=Img.^2;
Int=integral(Img);%计算Img的积分图像
Int2=integral(Img2);%计算Img^2的积分图像
Var=zeros(m,n);
N=(2*d+1)^2;%局部像素点个数
for i=1:m
    for j=1:n
        i1=i+d+1;
        j1=j+d+1;
        %利用积分图像求局部和
        sumi2=Int2(i1+d,j1+d)+Int2(i1-d-1,j1-d-1)-Int2(i1+d,j1-d-1)-Int2(i1-d-1,j1+d);
        sumi=Int(i1+d,j1+d)+Int(i1-d-1,j1-d-1)-Int(i1+d,j1-d-1)-Int(i1-d-1,j1+d);
        Var(i,j)=(sumi2-sumi^2/N)/(N-1);
    end
end
Std=sqrt(Var);


function I=Integral(Img) 
Img=double(Img);
[m,n]=size(Img);
I=zeros(m,n);
for i=1:m
    for j=1:n
        if i==1 && j==1             %积分图像左上角
            I(i,j)=Img(i,j);
        elseif i==1 && j~=1         %积分图像第一行
            I(i,j)=I(i,j-1)+Img(i,j);
        elseif i~=1 && j==1         %积分图像第一列
            I(i,j)=I(i-1,j)+Img(i,j);
        else                        %积分图像其它像素
            I(i,j)=Img(i,j)+I(i-1,j)+I(i,j-1)-I(i-1,j-1);  
        end
    end
end

main.m:

clear all;
close all;
clc;
Img=double(imread('lena.tif'));
d=1;
tic
J1=Localstd(Img,d);
toc
tic
J2=fastLocalstd(Img,d);
toc
figure;
imshow([Img/max(Img(:)),J1/max(J1(:)),J2/max(J2(:))]);





lena.tif大小为256*256,Localstd运行时间为 0.984451s,而fastLocalstd运行时为0.025528s。

当然,对于局部标准差,Matlab本身提供了函数stdfilt,下面是它的核心代码:

function Std=Stdfilt(I,d)
I=double(I);
h=ones(2*d+1);
n=(2*d+1)^2;
conv1 = imfilter(I.^2,h,'symmetric') / (n-1); 
conv2 = imfilter(I,h,'symmetric').^2 / (n*(n-1));
Std = sqrt(conv1-conv2);

可以看出,它按照化简后的公式直接对图像进行卷积操作,因为imfilter函数是使用c实现的,且内部应该进行了优化,故速度较快。对上面的lena.tif运行时间为:0.064522s。


版权声明:本文为博主原创文章,未经博主允许不得转载。

转载于:https://www.cnblogs.com/luo-peng/p/4646216.html

2019-12-25 13:32:21 qq_39622795 阅读数 26

基于局部统计可变阈值的图像分割

  1. 灰度图上进行
  2. 每个像素点处设一个阈值
  3. 领域
  4. 一般使用领域内标准差,标准差表示对比度
  5. 全局平均还是领域平均看情况自己选择
  6. 求标准差、均值前一定要图像转换为float类型,tofloat不是matlab内置函数,需要自己添加
  7. 分割结果图均为二值图像

使用默认8领域

%默认8领域
I=imread('F:\20191214162428.jpg');
I1=rgb2gray(I);
I1=tofloat(I1);
I2=stdfilt(I1);
figure,imshow(I2);
m=imfilter(I1,1/9*ones(3),'replicate');
m2=mean2(I1);
I3=(I1>15*I2)&(I1>0.4*m);
figure,imshow(I3);

function [out,revertclass] = tofloat(inputimage)

% 类型转换
% 单纯的取 x 用的匿名函数句柄(玩的有点花)
identify = @(x) x;
% 将输入转换为单精度的函数句柄
tosingle = @im2single;

table = {'uint8',tosingle,@im2uint8 

         'uint16',tosingle,@im2uint16 

         'logical',tosingle,@logical

         'double',identify,identify

         'single',identify,identify};
% strcmp(s1,s2),输入字符串的数组可以说任意类型组合
% find 找出其中的真值
classIndex = find(strcmp(class(inputimage),table(:,1)));

if isempty(classIndex)

    error('不支持的图像类型');

end
% 找到对应的转换类型
out = table{classIndex,2}(inputimage);
% 记录对应逆转换类型
revertclass = table{classIndex,3};
end

原图:
在这里插入图片描述
局部标准差图像:
在这里插入图片描述
使用局部平均值的分割结果图:
在这里插入图片描述
使用全局平均值的分割结果图:
在这里插入图片描述
可见使用全局平均值的效果更好

使用自定义领域

%自定义领域
fil=[1,0,1,0,1;
    0,0,1,0,1;
    1,1,0,1,0;
    0,0,0,0,1;
    1,1,1,1,0];
I=imread('F:\20191214162428.jpg');
I1=rgb2gray(I);
I1=tofloat(I1);
I2=stdfilt(I1,fil);
figure,imshow(I2);
fil=fil/sum(fil(:));
m=imfilter(I1,fil,'replicate');
m2=mean2(I1);
I3=(I1>11*I2)&(I1>0.9*m);
figure,imshow(I3);

function [out,revertclass] = tofloat(inputimage)

% 类型转换
% 单纯的取 x 用的匿名函数句柄(玩的有点花)
identify = @(x) x;
% 将输入转换为单精度的函数句柄
tosingle = @im2single;

table = {'uint8',tosingle,@im2uint8 

         'uint16',tosingle,@im2uint16 

         'logical',tosingle,@logical

         'double',identify,identify

         'single',identify,identify};
% strcmp(s1,s2),输入字符串的数组可以说任意类型组合
% find 找出其中的真值
classIndex = find(strcmp(class(inputimage),table(:,1)));

if isempty(classIndex)

    error('不支持的图像类型');

end
% 找到对应的转换类型
out = table{classIndex,2}(inputimage);
% 记录对应逆转换类型
revertclass = table{classIndex,3};
end

分割结果图
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

2016-01-27 17:14:38 xiaotie1005 阅读数 4858


使用局部标准差实现图像的局部对比度增强算法。

Posted on 2013-09-16 15:01 Imageshop 阅读(426) 评论(1编辑 收藏

      图像的对比度增强算法在很多场合都有着重要的应用,特别是在医学图像上,这是因为在众多疾病的诊断中,医学图像的视觉检查时很有必要的。而医学图像由于本身及成像条件的限制,图像的对比度很低。因此,在这个方面已经开展了很多的研究。这种增强算法一般都遵循一定的视觉原则。众所周知,人眼对高频信号(边缘处等)比较敏感。虽然细节信息往往是高频信号,但是他们时常嵌入在大量的低频背景信号中,从而使得其视觉可见性降低。因此适当的提高高频部分能够提高视觉效果并有利于诊断。

     在这一方面,传统的线性对比度拉升以及直方图均衡化是使用的最为广泛的全局图像增强方法。对比度拉升线性的调整了图像的动态范围,而直方图均衡化栖利用累计直方图分布概率重新映射图像的数据。这些方法虽然简单,但是都没有考虑到局部的信息。并且,全局直方图均衡化(GHE)还会产生使得一些噪音过度加强。

     在局部对比度增强方面,有两种方式是最为有名的,一种是自适应直方图均衡化(AHE),这个算法可以参考我的博文限制对比度自适应直方图均衡化算法原理、实现及效果。还有一种就是自适应对比度增强(ACE)。AHE算法使用局部的直方图的相关信息对数据进行映射。这改变了图像的对比度,但是需要大量的计算。后来有人利用了双线性差值技术克服了这个问题,首先将图像分块,然后分别计算这些快内部的映射关系。为了增强某一个像素点的值,映射关系通过与这个像素所在块相邻的四个块的映射关系差值获得。在这个算法中,仅仅需要一个块大小的参数(在我的博文中还对参数进行了扩展)。

      ACE算法采用了反锐化掩模技术,我们对此过程解释如下:首先图像被分成两个部分。一是低频的反锐化掩模(unsharp mask)部分,可以通过图像的低通滤波(平滑,模糊技术)获得。二是高频成分,可以过原图减去反锐化掩模获取。然后高频部分被放大(放大系数集为对比度增益CG)并加入到反锐化掩模中去,最后得到增强的图像。ACE算法的核心就是如何计算CG,这里将介绍两种简单的CG计算方法。

    在继续说明之前,先贴一幅这个算法处理的图给大家看看效果,免得都是通篇的文字。

   

                  原图                                     处理后的图

      在谈及CG之前,我们还是需要看看一些基础的工作。首先上面谈到了unsharp mask,他对应于图像的低频成分。对于具体的像素,一般可以通过计算以该像素为中心的局部区域的像素平均值来实现。我们假定x(i,j)是图像中某点的灰度值,局部区域的定义为:以(i,j)为中心,窗口大小为(2n+1)*(2n+1)的区域,其中n为一个整数。当然这个窗口区域也不一定就要是正方形。局部的平均值,也就是低频部分,可以用下式计算:

                    

      而局部方差为:

        

      上式中σx(i,j)就是所谓的局部标准差(LSD)。定义f(i,j)表示x(i,j)对应的增强后的像素值。则ACE算法可以表示如下:

                 

     其中的函数G(i,j)就是上文所讲的CG。一般情况下CG总是大于1的,这样高频成分就能得到增强。

     关于CG的取值,一种最简单的方式就是令其为常量,假定为C,一般C>1,这样式3就变为:

                    

     这种情况下,图像中所有的高频部分都被同等放大,可能有些高频部分会出现过增强的现象的。

   

                    原图                                                                      n=50,c=2                                                              n=50,c=3

    在上图中,分别使用了C=2及C=3的的情况,在C=3时,可见星球边缘的部分被过增强,出现成片的白色。

    一种解决的方案就是使用不同的增益。Lee等人提出了如下的方案:

            

     上式中,D是个常数,这样,CG是空间自适应的,并且和局部均方差成反比,在图像的边缘或者其他变化剧烈的地方,局部均方差比较大,因此CG的值就比较小,这样就不会产生振铃效应。然而,在平滑的区域,局部均方差就会很小,这样CG的值比较大,从而引起了噪音的放大,所以需要对CG的最大值做一定的限制才能获得更好的效果。

     D这个常数的取值有些文章介绍说用图像的全局平均值,我编程的时候是用的图像的全局均方差,并且增加一个参数Amount,用来再次控制高频增强的程度。

     对式我们进行了编程,并测试一些图像,当进行一些参数的调整时,都能获得比较理想的增强效果,举例如下:

  

                              原图                                                     n=20,Amount=100,MaxCG=3                                 n=20,Amount=100,MaxCG=10

  

                              原图                                                   n=25,Amount=100,MaxCG=2                                n=50,Amount=150,MaxCG=2

  

                             原图                                                   n=50,Amount=150,MaxCG=3                                n=50,Amount=150,MaxCG=10

     不仅是对医学图像,对于我日常生活中的一些图示很理想的图片,经过上述ACE算法在大部分情况下也可获得较好的效果:

 

 

  为什么要对CG的最大值进行限制,我们举下面这个图作为例子:

 

                                原图                                                                                                  n=50,Amount=150,CG不限制

   在上图中,图像右上角部位属于平滑区域,但是仔细观察有着微弱的层次感,如果不对CG进行限制,则由于LSD的值很小,导致CG的值特别大,从而产生像素饱和的现象(即像素的值超过了255)。如右图所示,解决这个问题一种方法就是增加n的值,即增加取样半径,通过取样半径的增加,使得以平滑区域为中心的窗口能够覆盖一些非平滑的区域,从而产生较大的LSD值。另外一种就是采用CG限幅的方式,效果如下:

    

                                 n=200,Amount=150,CG不限制                                                                          n=50,Amount=150,CG<3

     针对不同的图像,要想获得理想的效果,可能还是需要人工的去参与或调解这些参数,我目前尚未发现能自动确定这些参数的文章。

     在代码的编写上,主要的难点就是如何实现快速的计算,因为随着取样半径n的增加,每个像素点涉及到的区域范围就越广,计算量会直线上升,因此需要采用优化的策略。

     均值模糊的优化已经有了很多成熟的方案,普遍的就是先行均值,在列均值,当然还有更好的优化方式。

     而计算局部均方差,则需要一定的水平了,但是同样的原理还是利用每变换一次都有很多的重叠的部分的数据,这个等有空了我会专门写篇文章的。 

     优化后的算法执行时间和n的没有关系的。

     同样,提供个编译好的文件给有兴趣研究该算法的朋友看看效果:http://files.cnblogs.com/Imageshop/LocalAreaContrastEnhance.rar

     当然,这个算法不是对任何图像都能产生正能量,那些本来就很好的图,你用它只会适得其反。

      

2015-06-06 13:31:49 u010839382 阅读数 5438

局部标准差在图像处理邻域具有广泛的应用,但是直接计算非常耗时,本文利用积分图像对局部标准差的计算进行加速。

局部标准差:

标准差定义如下(采用统计学中的定义,分母为):


其中

为了计算图像的局部标准差,首先设定局部区域的大小为 ,则局部区域的像素点个数

对标准差的公式进行化简:


,故:


可以看出,局部标准差计算中需要对图像进行局部求和操作,即

我们可以先分别计算出图像的积分图像,这样就能在常量时间计算出上述的局部和。

时间复杂度:

图像中共个像素点,每个局部区域共有个像素点,直接计算局部标准差的时间复杂度 ,利用积分图像计算局部标准差的时间复杂度为。(计算积分图像的时间为 、计算一个像素点的局部标准差需要常量时间,共个像素点,故最后的计算复杂度仍为 


Matlab程序:

Localstd.m:直接计算图像的局部标准差

function Std=Localstd(Img,d)
[m,n]=size(Img);
Var=zeros(m,n);
Img = padarray(Img,[d d],'symmetric');%边界填充
N=(2*d+1)^2;
for i=1:m
    for j=1:n
        i1=i+d;
        j1=j+d;
        B=Img(i1-d:i1+d,j1-d:j1+d);%取局部区域
        Var(i,j)=sum(sum((B-mean(B(:))).^2))/(N-1);%标准差公式
    end
end
Std=sqrt(Var);

fastLocalstd.m:利用积分图像进行快速计算

function Std=fastLocalstd(Img,d)
[m,n]=size(Img);
Img = padarray(Img,[d+1 d+1],'symmetric');%边界填充
Img2=Img.^2;
Int=integral(Img);%计算Img的积分图像
Int2=integral(Img2);%计算Img^2的积分图像
Var=zeros(m,n);
N=(2*d+1)^2;%局部像素点个数
for i=1:m
    for j=1:n
        i1=i+d+1;
        j1=j+d+1;
        %利用积分图像求局部和
        sumi2=Int2(i1+d,j1+d)+Int2(i1-d-1,j1-d-1)-Int2(i1+d,j1-d-1)-Int2(i1-d-1,j1+d);
        sumi=Int(i1+d,j1+d)+Int(i1-d-1,j1-d-1)-Int(i1+d,j1-d-1)-Int(i1-d-1,j1+d);
        Var(i,j)=(sumi2-sumi^2/N)/(N-1);
    end
end
Std=sqrt(Var);


function I=Integral(Img) 
Img=double(Img);
[m,n]=size(Img);
I=zeros(m,n);
for i=1:m
    for j=1:n
        if i==1 && j==1             %积分图像左上角
            I(i,j)=Img(i,j);
        elseif i==1 && j~=1         %积分图像第一行
            I(i,j)=I(i,j-1)+Img(i,j);
        elseif i~=1 && j==1         %积分图像第一列
            I(i,j)=I(i-1,j)+Img(i,j);
        else                        %积分图像其它像素
            I(i,j)=Img(i,j)+I(i-1,j)+I(i,j-1)-I(i-1,j-1);  
        end
    end
end

main.m:

clear all;
close all;
clc;
Img=double(imread('lena.tif'));
d=1;
tic
J1=Localstd(Img,d);
toc
tic
J2=fastLocalstd(Img,d);
toc
figure;
imshow([Img/max(Img(:)),J1/max(J1(:)),J2/max(J2(:))]);





lena.tif大小为256*256,Localstd运行时间为 0.984451s,而fastLocalstd运行时为0.025528s。

当然,对于局部标准差,Matlab本身提供了函数stdfilt,下面是它的核心代码:

function Std=Stdfilt(I,d)
I=double(I);
h=ones(2*d+1);
n=(2*d+1)^2;
conv1 = imfilter(I.^2,h,'symmetric') / (n-1); 
conv2 = imfilter(I,h,'symmetric').^2 / (n*(n-1));
Std = sqrt(conv1-conv2);

可以看出,它按照化简后的公式直接对图像进行卷积操作,因为imfilter函数是使用c实现的,且内部应该进行了优化,故速度较快。对上面的lena.tif运行时间为:0.064522s。


2015-10-06 20:50:00 weixin_30505751 阅读数 8

基于“局部标准差”的图像增强(原理、算法、代码)

一、理论
         图像增强算法的基本原则是“降低低频区域,突出高频区域”,以此强化边缘,达到增强的目的。最简单的例子就是通过原始图像减去高斯模糊处理后的图像,就能够将边缘强化出来。
         直方图均衡化也是一种非常常见的增强方法。但是为了避免背景的干扰,更倾向于采用“局部”方法进行处理。我们这里着重研究自适应对比度增强(ACE)的相关内容。
        ACE的定义和原理看上去还是比较简单的。这里的都可以根据图像本身计算出来。而则需要单独计算。

          可以为单独的常量,或者通过来代替。这里的D是一个全局的值,比如平均值。

二、实现
        涉及到局部的运算,自然而然会想到使用卷积的方法。更好的是Opencv提供了专门的函数用来做这个工作—BLUR
文档中写到:
那么正是我们想要的结果。
//ace 自适应对比度均衡研究
//by  jsxyhelu
//感谢 imageshop
#include "stdafx.h"
#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace std;
using namespace cv;
//点乘法 elementWiseMultiplication
cv::Mat EWM(cv::Mat m1,cv::Mat m2){
    Mat dst=m1.mul(m2);
    return dst;
}
void main()
{
    Mat src = imread("hand.jpg",0);
    Mat meanMask;
    Mat varMask;
    Mat meanGlobal;
    Mat varGlobal;
    Mat dst;
    Mat tmp;
    Mat tmp2;
    int C = 30;
    int D = 133;
    //全局均值和均方差
    blur(src.clone(),meanGlobal,src.size());
    varGlobal = src - meanGlobal;
    varGlobal = EWM(varGlobal,varGlobal);
    blur(src.clone(),meanMask,Size(50,50));//meanMask为局部均值
    tmp = src - meanMask;                        
    varMask = EWM(tmp,tmp);              
    blur(varMask,varMask,Size(50,50));    //varMask为局部方差
    
    dst = meanMask + C*tmp;
    imshow("src",src);
    imshow("dst",dst);
     
    waitKey();
}
接下来,为了实现那么需要计算局部标准差和全局均值或方差
前面已经计算出了局部均值,那么
tmp = src - meanMask;  
    varMask = EWM(tmp,tmp);         
    blur(varMask,varMask,Size(50,50));    //varMask为局部方差   
计算出局部方差
//换算成局部标准差
    varMask.convertTo(varMask,CV_32F);
    for (int i=0;i<varMask.rows;i++){
        for (int j=0;j<varMask.cols;j++){
            varMask.at<float>(i,j) =  (float)sqrt(varMask.at<float>(i,j));
        }
    }
换算成局部标准差
meanStdDev(src,meanGlobal,varGlobal); //meanGlobal为全局均值 varGlobal为全局标准差
是opencv提供的全局均值和标准差计算函数。
全部代码进行重构后如下
//ace 自适应对比度均衡研究
//by  jsxyhelu
//感谢 imageshop
#include "stdafx.h"
#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace std;
using namespace cv;
//点乘法 elementWiseMultiplication
cv::Mat EWM(cv::Mat m1,cv::Mat m2){
    Mat dst=m1.mul(m2);
    return dst;
}
//图像局部对比度增强算法
cv::Mat ACE(cv::Mat src,int C = 4,int n=20,int MaxCG = 5){
    Mat meanMask;
    Mat varMask;
    Mat meanGlobal;
    Mat varGlobal;
    Mat dst;
    Mat tmp;
    Mat tmp2;
    blur(src.clone(),meanMask,Size(50,50));//meanMask为局部均值 
    tmp = src - meanMask;  
    varMask = EWM(tmp,tmp);         
    blur(varMask,varMask,Size(50,50));    //varMask为局部方差   
    //换算成局部标准差
    varMask.convertTo(varMask,CV_32F);
    for (int i=0;i<varMask.rows;i++){
        for (int j=0;j<varMask.cols;j++){
            varMask.at<float>(i,j) =  (float)sqrt(varMask.at<float>(i,j));
        }
    }
    meanStdDev(src,meanGlobal,varGlobal); //meanGlobal为全局均值 varGlobal为全局标准差
    tmp2 = varGlobal/varMask;
    for (int i=0;i<tmp2.rows;i++){
        for (int j=0;j<tmp2.cols;j++){
            if (tmp2.at<float>(i,j)>MaxCG){
                tmp2.at<float>(i,j) = MaxCG;
            }
        }
    }
    tmp2.convertTo(tmp2,CV_8U);
    tmp2 = EWM(tmp2,tmp);
    dst = meanMask + tmp2;
    imshow("D方法",dst);
    dst = meanMask + C*tmp;
    imshow("C方法",dst);
    return dst;
}
void main()
{
    Mat src = imread("plant.bmp",0); 
    imshow("src",src);
    ACE(src);
    waitKey();
}
三、小结
      从结果上来看,ACE算法对于特定情况下的图片细节增强是显著的,但是并不是适用于所有的情况,并且其参数需要手工进行调整。了解它的特性,就能够解决一系列的问题,有效地增强现实。
 
 





转载于:https://www.cnblogs.com/jsxyhelu/p/4857721.html

图像对焦区域检测

阅读数 1553

没有更多推荐了,返回首页