R аварийно завершает работу / прерывает работу, используя Rcpp с вводом NA

Я хочу обработать два растровых изображения (Ra и Rb), где Ra — это значение пикселя, а Rb — значения его соседей. взяв сумму в качестве примера, предполагая 3 * 3 соседей, для каждого пикселя в Ra я добавлю свое значение к значениям соседних пикселей в Rb, и, наконец, я получу другое изображение.

Растровый пакет R предоставляет функцию фокусировки, которая работает только на одном входе изображения, я попытался изменить код C ++ (введите описание ссылки здесь) принять два ввода изображения, используя Rcpp. Модифицированный код работает хорошо, если во входном изображении Rb нет пропущенных значений. Однако R всегда прерывается, если в Rb есть NA. В частности, прервать на втором или третьем тесте. это может быть похоже на эта почта. однако, это не вылетало, если нет NA на входе Rb. Кажется, я не справился с АН правильно. У меня нет глубоких знаний о C ++. Может кто-нибудь помочь мне проверить это?

вот мой файл cpp:

#include <Rcpp.h>
#include <R.h>
#include <Rinternals.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <Rmath.h>
#include "Rdefines.h"#include "R_ext/Rdynload.h"
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector focal_quantile(NumericVector xd, int ngbb, NumericVector sf) {
//the imges are transfered to vector, ngbb is the size of the window
R_len_t i, j, k, q;
int wrows = ngbb;
int wcols = ngbb;
int wn = wrows * wcols;

int nrow = 6;//the input raste has 6 rows
int ncol = 7;//the input raste has 7 cols

int n = nrow * ncol;
NumericVector xans(n);
NumericVector xx(wn);

int wr = floor(wrows / 2);
int wc = floor(wcols / 2);

int nwc = ncol - wc - 1;
int col = 0;

// first rows
for (i = 0; i < ncol*wr; i++) {// the first row, the resutl is set as NA as the neighbor does not have nine values
xans[i] = R_NaReal;
}

for (i = ncol*wr; i < (ncol * (nrow-wr)); i++) {//start from the second row
col = i % ncol;
if ((col < wc) | (col > nwc)) {//the first pixel of the second is also set as NA
xans[i] = R_NaReal;
} else {// to get the nine values in the 3*3 windows
q = 0;
for (j = -wr; j <= wr; j++) {
for (k = -wc; k <= wc; k++) {
xx[q] = xd[j * ncol + k + i];
q++;
}
}
xx = na_omit(xx);
int n_qt = xx.size();
if (n_qt > 0){//
xans[i]=sum(xx)+100*sf[i];// here is the calculation, my goal is more complicated than this example
} else {
xans[i] = R_NaReal;//R_NaReal
}

}
}
// last rows
for (i = ncol * (nrow-wr); i < n; i++) {
xans[i] = R_NaReal;
}
return(xans);
}

Затем скомпилируйте его с помощью sourceCpp

сгенерировать пример данных для проверки

  rr=raster(nrow=6,ncol=7)## example for Ra
projection(rr)="+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"rr[]=(2:43)*10
rrqt=rr/43 ## example for Rb
##it works fine, if there is no NA in Ra
#rr[1:10]=NA #window of global enviornment is refleshing and then aborts with such NAs
focal_quantile(rr[],3,rrqt[])

Пример результатов

 [1]       NA       NA       NA       NA       NA       NA       NA       NA 118918.6 130810.5 142702.3 154594.2 166486.0       NA       NA
[16] 202161.6 214053.5 225945.3 237837.2 249729.1       NA       NA 285404.7 297296.5 309188.4 321080.2 332972.1       NA       NA 368647.7
[31] 380539.5 392431.4 404323.3 416215.1       NA       NA       NA       NA       NA       NA       NA       NA

полученный NA приемлем, так как в окнах нет девяти значений.
Для такого примера я изменяю значения растра rr (без NA). это работает гладко. когда я введу NA в rr, например шестую строку кодов выше. окно глобальной среды обновляется, а Rstudio прерывается.

информация о сеансе

R version 3.3.0 (2016-05-03)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] Rcpp_0.12.11 raster_2.5-8 sp_1.2-3

loaded via a namespace (and not attached):
[1] rgdal_1.2-5     tools_3.3.0     grid_3.3.0      lattice_0.20-35

Большое спасибо!