【发布时间】:2020-09-20 14:35:27
【问题描述】:
我在推导 NDVI(归一化差异植被指数),它是 (NIR-R)/(NIR+R) 的比率,其中 NIR 是近红外波段,R 是红波段。这个指数范围从-1到1。所以我写了一个pyopencl代码,这就是我所做和观察到的。
Python 代码:
import pyopencl as cl
import cv2
from PIL import Image
import numpy as np
from time import time
import matplotlib.pyplot as plt
#get kernel file
def getKernel():
kernel = open('kernel.c').read()
return kernel
#return images as numpy int32 arrays
def convToArray(im_r,im_nir):
a = np.asarray(im_r).astype(np.int32)
b = np.asarray(im_nir).astype(np.int32)
return a,b
#processing part
def getDerivation(platform,device,im_r,im_nir):
#setting device
pltfrm = cl.get_platforms()[platform]
dev = pltfrm.get_devices()[device]
cntx = cl.Context([dev])
queue = cl.CommandQueue(cntx)
#get 2Darrays
r,nir = convToArray(im_r,im_nir)
#shape of array
x = r.shape[1]
mf = cl.mem_flags
bs = time()
#input images buffer
inR = cl.Buffer(cntx,mf.READ_ONLY | mf.COPY_HOST_PTR,hostbuf=r)
inIR = cl.Buffer(cntx,mf.READ_ONLY | mf.COPY_HOST_PTR,hostbuf=nir)
#output image buffers
ndvi = cl.Buffer(cntx,mf.WRITE_ONLY,r.nbytes)
be = time()
print("Buffering time: " + str(be-bs) + " sec")
ts = time()
#load kernel
task = cl.Program(cntx,getKernel()%(x)).build()
#execute the process
task.derive(queue,r.shape,None,inR,inIR,ndvi)
#create empty buffer to store result
Vout = np.empty_like(r)
#dump output buffers to empty arrays
cl.enqueue_copy(queue,Vout,ndvi)
te = time()
#convert arrays to gray - image compatible formate
NDVI = Vout.astype(np.uint8)
print("Processing time: " + str(te - ts) + " On: " + pltfrm.name + " --> " + dev.name)
return NDVI
def process(platform,device,im_r,im_nir):
NDVI,NDBI,NDWI = getDerivation(platform,device,im_g,im_r,im_nir,im_swir)
print(NDVI)
cv2.imshow("NDVI",NDVI)
cv2.waitKey(0)
if __name__ == '__main__':
R = cv2.imread("BAND3.jpg",0)
NIR = cv2.imread("BAND4.jpg",0)
print(R.dtype) #returns uint8
process(0,0,R,NIR) #(0,0) is my intel gpu
内核代码(C):
__kernel void derive(__global int* inR,__global int* inIR,__global int* ndvi){
int x = get_global_id(0);
int y = get_global_id(1);
int width = %d;
int index = x + y*width;
//ndvi ratio (-1 to 1)
int a = ((inIR[index] - inR[index])/(inIR[index] + inR[index])) * (256);
a = (a < (0) ? (-1*a) : (a));
a = (a > (255) ? (255) : (a));
ndvi[index] = (a);
}
输入图片R:
输入图像近红外:
两张图片的位深均为 8
但我得到的只是一张空白图片。我最初出于调试原因将结果写在命令行上, 命令行输出:
(1151, 1151)
Buffering time: 0.015959739685058594 sec
Processing time: 0.22115755081176758 On: Intel(R) OpenCL --> Intel(R) HD Graphics 520
[[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]
...
[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]]
现在我认为我可能没有为图像使用正确的数据类型?此外,在内核中,((inIR[index] - inR[index])/(inIR[index] + inR[index])) 行将给出一个浮点值,我将其乘以 256 以获得相应浮点值的像素值。那么有问题吗?有谁知道我哪里出错了?
非常感谢您的帮助!
【问题讨论】:
-
你为什么认为你会得到一个浮动?你会得到一个整数。在分割之前尝试将顶部或底部转换为浮动。
-
@MarkSetchell
float num = (inIR[index] - inR[index]); float deno = (inIR[index] + inR[index]); int a = (int)((num/deno) * (256));你在说这个吗? -
另外,您可能不想将范围为 [-1,+1] 的 NDVI 乘以 256。您可能会加 1 以使其范围为 [0,2 ] 并乘以 127.5 使其达到 [0,255] 范围。
-
float num = (float)(inIR[index] - inR[index]); float deno = (float)(inIR[index] + inR[index]); int a = (int)(((num/deno)+1)*127.5);我这样做了,但结果是一样的......空白数组。顺便说一句,这就是你所说的吗? -
我暂时无法尝试,但尝试直接使用 Numpy 进行比较。这只是一行
NDVI = (im_nir.astype(np.float)-im_r.astype(np.float))/...然后看看NDVI.max()和NDVI.min()。
标签: c python-3.x image-processing pyopencl