Minus values of distortion coefficients cause broken images with OpenGL Vertex Shader

Question

I'm working with python moderngl to generate a simulated scene views from point clouds made from Digital Elevation Model and airborne photos. I implemented a camera model based on OpenCV in GLSL Vertex Shader and it works well usually.

However, some minus values of radial distortion coefficients (k1 or k2) cause broken images.

Here below is the code of Vertex Shader and Fragment Shader.

    prog = ctx.program(
     vertex_shader='''
         #version 330
         in vec3 in_vert;
         in vec3 in_color;
         out vec4 v_color;

         // decrare some values used inside GPU by "uniform"
         // the real values will be later set by CPU
         uniform mat4 proj; // projection matrix
         uniform mat4 view; // model view matrix
         uniform mat4 rotate; // rotation matrix for normals
         uniform float dist_coeffs[5]; // distortion coefficients 
         // distortion coefficients. k1x, k1y, k2x, k2y, p1, p2
         // radial distortion parameters are divided as x-axis and y-axis
         vec4 distort(vec4 view_pos){
          // normalize
          float z =  view_pos[2];
          float z_inv = 1 / z;
          float x1 = view_pos[0] * z_inv;
          float y1 = view_pos[1] * z_inv;

          // precalculations
          float x1_2 = x1*x1;
          float y1_2 = y1*y1;
          float x1_y1 = x1*y1;
          float r2 = x1_2 + y1_2;
          float r4 = r2*r2;

          // radial distortion factor
          float r_dist_x = (1 + dist_coeffs[0] * r2 + dist_coeffs[2] * r4);
          float r_dist_y = (1 + dist_coeffs[1] * r2 + dist_coeffs[3] * r4);

          // full (radial + tangential) distortion
          float x2 = x1*r_dist_x + 2*dist_coeffs[4]*x1_y1 + dist_coeffs[5]*(r2 + 2*x1_2);
          float y2 = y1*r_dist_y + 2*dist_coeffs[5]*x1_y1 + dist_coeffs[4]*(r2 + 2*y1_2);

          // denormalize for projection
          return vec4(x2*z, y2*z, z, view_pos[3]);
          }

         void main() {
             vec4 local_pos = vec4(in_vert, 1.0);
             vec4 view_pos = vec4(view * local_pos);
             vec4 dist_pos = distort(view_pos);
             vec4 v_norm = rotate * vec4(in_color, 1.0); // normals are transformed to camera CRS, -1~1
             v_color = (v_norm + 1) / 2; // normals as rgb, 0~1
             gl_Position = vec4(proj * dist_pos);
         }
     ''',
     fragment_shader='''
         #version 330
         in vec4 v_color;
         layout(location=0)out vec4 f_color;
         void main() {
             f_color = vec4(v_color); // 1,0 added is alpha
         }
     '''
    )

And here below I show the R program which makes modelview, and perspective projection matrix.

# camera intrinsic matrix
projection_mat <- function(fov_x_degree, w, h, near = -1, far = 1
                           , cx = w/2, cy = h/2){
  fov_x <- fov_x_degree * pi /180
  fov_y <- fov_x * h / w
  fx <- 1/(tan(fov_x/2))
  fy <- 1/(tan(fov_y/2))
  matrix(c(fx , 0, (w-2*cx)/w, 0, 
           0, fy, -(h-2*cy)/h, 0,
           0, 0, -(far+near)/(far-near), -2*far*near/(far-near),
           0, 0, -1, 0), 
         ncol = 4, 
         byrow = F) %>% 
    as.vector() %>% 
    reticulate::np_array(dtype = "float32")
}

# camera extrinsic parameter
modelview_mat <- function(pan_degree, tilt_degree, roll_degree, t_x, t_y, t_z){
  y <-  (360 - pan_degree) * pi / 180
  x <-  (tilt_degree) * pi / 180 
  z <-  (roll_degree) * pi / 180

  rmat_z <- c(cos(z), -sin(z), 0, 0,
              sin(z), cos(z), 0, 0,
              0, 0, 1, 0,
              0, 0, 0, 1) %>% matrix(ncol = 4, byrow = T)
  rmat_x <- c(1, 0, 0, 0,
              0, cos(x), -sin(x), 0,
              0, sin(x), cos(x), 0,
              0, 0, 0, 1) %>% matrix(ncol = 4, byrow = T)
  rmat_y <- c(cos(y), 0, sin(y), 0,
              0, 1, 0, 0,
              -sin(y), 0, cos(y), 0,
              0, 0, 0, 1) %>%matrix(ncol = 4, byrow = T)
  rmat <- rmat_x %*% rmat_z %*% rmat_y
  tmat <- matrix(c(1, 0, 0, -t_x,
                   0, 1, 0, -t_z,
                   0, 0, 1, -t_y,
                   0, 0, 0, 1), byrow = T, ncol = 4)
  rmat %*% tmat %>%
    as.vector() %>% 
    reticulate::np_array(dtype = "float32")
}

this is a generated image without error

and this is the broken one