Adding Random valued impulse noise to an image

function [ image_noisy ] = addRVIN( image_input, ND, noise_width )
%ADDRVIN adds Random valued impulse noise to an image
% image_input Input noise free image
% ND is noise density that ranges from 0 to 1
% noise_width is width of random valued impulse noise
[~,x] = imhist(image_input);
init_width = x(noise_width);
image_noisy = image_input;
Narr = rand(size(image_noisy));
p = ND/2;
% image_noisy = (image_input + p*rand(size(image_input)))/(1+p);
image_rand = rand(size(image_input));
range_init = init_width;
image_randInit = image_rand*range_init;
% image_noisy(Narr<ND/2) = (init_width + image_rand(Narr<ND/2));
image_noisy(Narr<ND/2) = (image_randInit(Narr<ND/2));
% image_noisy(Narr<ND/2) = (noise_width + p*rand(size(image_input)))/(1+p);
image_randFinal = image_rand*range_init+ x(end) - x(noise_width+1);
image_noisy((Narr>=ND/2)&(Narr<ND)) = image_randFinal((Narr>=ND/2)&(Narr<ND));
% figure, imhist(image_input)
% figure, imhist(image_noisy)

end