草庐IT

计算方差图像python

codeneng 2023-03-28 原文

Calculating variance image python

有没有一种简单的方法可以使用 Python/NumPy/Scipy 计算图像上的运行方差过滤器?通过运行方差图像,我的意思是计算图像中每个子窗口 I 的 sum((I - mean(I))^2)/nPixels 的结果。

由于图像非常大(12000x12000 像素),我想避免在格式之间转换数组的开销,以便能够使用不同的库然后再转换回来。

我想我可以通过使用类似

的方法找到平均值来手动执行此操作

1
2
3
4
5
kernel = np.ones((winSize, winSize))/winSize**2
image_mean = scipy.ndimage.convolve(image, kernel)
diff = (image - image_mean)**2
# Calculate sum over winSize*winSize sub-images
# Subsample result

但是如果有像 Matlab 中的 stdfilt 函数这样的东西会更好。

任何人都可以指出具有此功能并支持 numpy 数组的库的方向,或者提示/提供一种在 NumPy/SciPy 中执行此操作的方法吗?

  • 也许这个问题是一个起点。


更简单的解决方案也更快:使用 SciPy\\'s ndimage.uniform_filter

1
2
3
4
5
6
7
8
9
import numpy as np
from scipy import ndimage
rows, cols = 500, 500
win_rows, win_cols = 5, 5

img = np.random.rand(rows, cols)
win_mean = ndimage.uniform_filter(img, (win_rows, win_cols))
win_sqr_mean = ndimage.uniform_filter(img**2, (win_rows, win_cols))
win_var = win_sqr_mean - win_mean**2

"跨步技巧"是个漂亮的技巧,但速度慢了 4 且不那么可读。
generic_filter 比步幅慢 20 倍...

  • 不错的解决方案!不过,基于 OpenCV 的解决方案甚至比基于 ndimage 的解决方案还要快。我根据我的基于 OpenCV 的代码对您的代码进行了计时,示例中的 OpenCV 版本快 4.1 倍。


您可以使用 numpy.lib.stride_tricks.as_strided 来获得图像的窗口视图:

1
2
3
4
5
6
7
8
9
10
11
import numpy as np
from numpy.lib.stride_tricks import as_strided

rows, cols = 500, 500
win_rows, win_cols = 5, 5

img = np.random.rand(rows, cols)

win_img = as_strided(img, shape=(rows-win_rows+1, cols-win_cols+1,
                                 win_rows, win_cols),
                     strides=img.strides*2)

现在 win_img[i, j](win_rows, win_cols) 数组,左上角位于 [i, j] 位置:

1
2
3
4
5
6
7
8
9
10
11
12
>>> img[100:105, 100:105]
array([[ 0.34150754,  0.17888323,  0.67222354,  0.9020784 ,  0.48826682],
       [ 0.68451774,  0.14887515,  0.44892615,  0.33352743,  0.22090103],
       [ 0.41114758,  0.82608407,  0.77190533,  0.42830363,  0.57300759],
       [ 0.68435626,  0.94874394,  0.55238567,  0.40367885,  0.42955156],
       [ 0.59359203,  0.62237553,  0.58428725,  0.58608119,  0.29157555]])
>>> win_img[100,100]
array([[ 0.34150754,  0.17888323,  0.67222354,  0.9020784 ,  0.48826682],
       [ 0.68451774,  0.14887515,  0.44892615,  0.33352743,  0.22090103],
       [ 0.41114758,  0.82608407,  0.77190533,  0.42830363,  0.57300759],
       [ 0.68435626,  0.94874394,  0.55238567,  0.40367885,  0.42955156],
       [ 0.59359203,  0.62237553,  0.58428725,  0.58608119,  0.29157555]])

但是,您必须小心,不要将图像的窗口视图转换为它的窗口副本:在我的示例中,这将需要 25 倍的存储空间。我相信 numpy 1.7 可以让您选择多个轴,因此您可以简单地执行以下操作:

1
>>> np.var(win_img, axis=(-1, -2))

我被 numpy 1.6.2 卡住了,所以我无法测试它。如果我没记错我的数学,另一个选项可能会因窗口不太大而失败:

1
2
3
>>> win_mean = np.sum(np.sum(win_img, axis=-1), axis=-1)/win_rows/win_cols
>>> win_sqr_mean = np.sum(np.sum(win_img**2, axis=-1), axis=-1)/win_rows/win_cols
>>> win_var = win_sqr_mean - win_mean**2

现在 win_var 是一个形状数组

1
2
>>> win_var.shape
(496, 496)

win_var[i, j] 保存 (5, 5) 窗口的方差,左上角在 [i, j]

  • img.shape = (1,1) 时,我得到 win_img 的空数组
  • 例如,如何在每个角落的像素周围获得一个窗口?


经过一些优化,我们想出了一个通用 3D 图像的函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def variance_filter( img, VAR_FILTER_SIZE ):
  from numpy.lib.stride_tricks import as_strided

  WIN_SIZE=(2*VAR_FILTER_SIZE)+1
  if ~ VAR_FILTER_SIZE%2==1:
      print 'Warning, VAR_FILTER_SIZE must be ODD Integer number  '
  # hack -- this could probably be an input to the function but Alessandro is lazy
  WIN_DIMS = [ WIN_SIZE, WIN_SIZE, WIN_SIZE ]


  # Check that there is a 3D image input.
  if len( img.shape ) != 3:
      print"\\t variance_filter: Are you sure that you passed me a 3D image?"
      return -1
  else:
      DIMS = img.shape

  # Set up a windowed view on the data... this will have a border removed compared to the img_in
  img_strided = as_strided(img, shape=(DIMS[0]-WIN_DIMS[0]+1, DIMS[1]-WIN_DIMS[1]+1, DIMS[2]-WIN_DIMS[2]+1, WIN_DIMS[0], WIN_DIMS[1], WIN_DIMS[2] ), strides=img.strides*2)

  # Calculate variance, vectorially
  win_mean = numpy.sum(numpy.sum(numpy.sum(img_strided, axis=-1), axis=-1), axis=-1) / (WIN_DIMS[0]*WIN_DIMS[1]*WIN_DIMS[2])

  # As per http://en.wikipedia.org/wiki/Variance, we are removing the mean from every window,
  #   then squaring the result.
  # Casting to 64 bit float inside, because the numbers (at least for our images) get pretty big
  win_var = numpy.sum(numpy.sum(numpy.sum((( img_strided.T.astype('<f8') - win_mean.T.astype('<f8') )**2).T, axis=-1), axis=-1), axis=-1) / (WIN_DIMS[0]*WIN_DIMS[1]*WIN_DIMS[2])

  # Prepare an output image of the right size, in order to replace the border removed with the windowed view call
  out_img = numpy.zeros( DIMS, dtype='<f8' )
  # copy borders out...
  out_img[ WIN_DIMS[0]/2:DIMS[0]-WIN_DIMS[0]+1+WIN_DIMS[0]/2, WIN_DIMS[1]/2:DIMS[1]-WIN_DIMS[1]+1+WIN_DIMS[1]/2, WIN_DIMS[2]/2:DIMS[2]-WIN_DIMS[2]+1+WIN_DIMS[2]/2, ] = win_var

  # output
  return out_img.astype('>f4')

您可以使用 scipy.ndimage.generic_filter。我无法使用 matlab 进行测试,但也许这可以为您提供所需的内容:

1
2
3
4
5
import numpy as np
import scipy.ndimage as ndimage    
subs = 10  # this is the size of the (square) sub-windows
img = np.random.rand(500, 500)
img_std = ndimage.filters.generic_filter(img, np.std, size=subs)

您可以使用 footprint 关键字制作任意大小的子窗口。有关示例,请参阅此问题。

  • 虽然可以轻松访问所有漂亮的边缘处理功能很好,但这是诸如 vectorize 之类的欺骗性 numpy 函数之一,它们基本上在您的图像上运行 python 循环:在带有 (10, 10) 窗口的 (500, 500) 示例中, generic_filter 在我的回答中比窗口视图方法慢 x10。
  • @Jaime:这是真的。我还没有真正理解为什么它这么慢。但是 OP 想要一个命令 ;-)

有关计算方差图像python的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  2. ruby-on-rails - 使用一系列等级计算字母等级 - 2

    这里是Ruby新手。完成一些练习后碰壁了。练习:计算一系列成绩的字母等级创建一个方法get_grade来接受测试分数数组。数组中的每个分数应介于0和100之间,其中100是最大分数。计算平均分并将字母等级作为字符串返回,即“A”、“B”、“C”、“D”、“E”或“F”。我一直返回错误:avg.rb:1:syntaxerror,unexpectedtLBRACK,expecting')'defget_grade([100,90,80])^avg.rb:1:syntaxerror,unexpected')',expecting$end这是我目前所拥有的。我想坚持使用下面的方法或.join,

  3. ruby-on-rails - 添加回形针新样式不影响旧上传的图像 - 2

    我有带有Logo图像的公司模型has_attached_file:logo我用他们的Logo创建了许多公司。现在,我需要添加新样式has_attached_file:logo,:styles=>{:small=>"30x15>",:medium=>"155x85>"}我是否应该重新上传所有旧数据以重新生成新样式?我不这么认为……或者有什么rake任务可以重新生成样式吗? 最佳答案 参见Thumbnail-Generation.如果rake任务不适合你,你应该能够在控制台中使用一个片段来调用重新处理!关于相关公司

  4. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  5. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  6. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  7. python - 如何读取 MIDI 文件、更改其乐器并将其写回? - 2

    我想解析一个已经存在的.mid文件,改变它的乐器,例如从“acousticgrandpiano”到“violin”,然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容,该乐器通过program_change或patch_change指令进行了更改,但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作,但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成,每个音轨是十六个channel中任何一个上的

  8. 「Python|Selenium|场景案例」如何定位iframe中的元素? - 2

    本文主要介绍在使用Selenium进行自动化测试或者任务时,对于使用了iframe的页面,如何定位iframe中的元素文章目录场景描述解决方案具体代码场景描述当我们在使用Selenium进行自动化测试的时候,可能会遇到一些界面或者窗体是使用HTML的iframe标签进行承载的。对于iframe中的标签,如果直接查找是无法找到的,会抛出没有找到元素的异常。比如近在咫尺的例子就是,CSDN的登录窗体就是使用的iframe,大家可以尝试通过F12开发者模式查看到的tag_name,class_name,id或者xpath来定位中的页面元素,会抛出NoSuchElementException异常。解决

  9. 计算机毕业设计ssm+vue基本微信小程序的小学生兴趣延时班预约小程序 - 2

    项目介绍随着我国经济迅速发展,人们对手机的需求越来越大,各种手机软件也都在被广泛应用,但是对于手机进行数据信息管理,对于手机的各种软件也是备受用户的喜爱小学生兴趣延时班预约小程序的设计与开发被用户普遍使用,为方便用户能够可以随时进行小学生兴趣延时班预约小程序的设计与开发的数据信息管理,特开发了小程序的设计与开发的管理系统。小学生兴趣延时班预约小程序的设计与开发的开发利用现有的成熟技术参考,以源代码为模板,分析功能调整与小学生兴趣延时班预约小程序的设计与开发的实际需求相结合,讨论了小学生兴趣延时班预约小程序的设计与开发的使用。开发环境开发说明:前端使用微信微信小程序开发工具:后端使用ssm:VU

  10. ruby-on-rails - 在 Ruby (on Rails) 中使用 imgur API 获取图像 - 2

    我正在尝试使用Ruby2.0.0和Rails4.0.0提供的API从imgur中提取图像。我已尝试按照Ruby2.0.0文档中列出的各种方式构建http请求,但均无济于事。代码如下:require'net/http'require'net/https'defimgurheaders={"Authorization"=>"Client-ID"+my_client_id}path="/3/gallery/image/#{img_id}.json"uri=URI("https://api.imgur.com"+path)request,data=Net::HTTP::Get.new(path

随机推荐