R/fitBinImage.R

In BayesBD: Bayesian Inference for Image Boundaries

fitBinImage = function(image, gamma.fun = NULL, center = NULL, inimean = NULL, nrun, nburn, J, ordering, mask = NULL, slice, outputAll){

	is.string <- function(input) {
    		is.character(input) & length(input) == 1
	}
	if(is.string(image)){
		if(any(substr(image,start = nchar(image)-2,stop = nchar(image))=="jpg",substr(image,start = nchar(image)-3,stop = nchar(image))=="jpeg")){
			img = readJPEG(image)
		}else if(substr(image,start = nchar(image)-2,stop = nchar(image))=="png"){
			img = readPNG(image)
		}else {
			return('The image file must have a .png or .jpeg/.jpg file type.')
		}
		
		if(length(dim(img))>2){
			img = matrix(img[,,1],dim(img)[1],dim(img)[2])
		}else {
			img = img
		}

		n1 = nrow(img)
		n2 = ncol(img)

		img_flip = img
	
		for(i in 1:n1){
			for(j in 1:n2){
				img_flip[i,j] = img[n1-i+1,j]
			}
		}

		img = img_flip

		if(any(center[1]>n2,center[2]>n1,center[1]<0,center[2]<0)){
			return(paste('The center should be a pixel (x,y) between 0 and ', ncol(img), 'for x, and 0 and ', nrow(img), ' for y.'))
		}

		r.obs = img
		theta.obs = img
		intensity = img
		for(i in 1:n1){
			for(j in 1:n2){
				r.obs[i,j] = sqrt(((i - center[2])/n1)^2 + ((j - center[1])/n2)^2)
				theta.obs[i,j] = atan2((i - center[2])/n1, (j - center[1])/n2)
				theta.obs[theta.obs < 0] = theta.obs[theta.obs < 0] + 2*pi
			}
		}

		r.obs = as.vector(r.obs)
		theta.obs = as.vector(theta.obs)
		intensity = as.vector(intensity)

		obs = list(r.obs = r.obs, theta.obs = theta.obs, intensity = intensity, center=center)

		if(is.null(mask)){
			mask = rep(1,length(obs$intensity))
		}else {
			center = obs$center
			new.r.obs = obs$r.obs[mask==1]
			new.theta.obs = obs$theta.obs[mask==1]
			new.intensity = obs$intensity[mask==1]
			obs = list(r.obs=new.r.obs, theta.obs = new.theta.obs, intensity=new.intensity, center=center)
		}

		if(is.null(inimean)){
			ini.mean.estimator = function(r){
				obs.in = obs$intensity[which(obs$r.obs<r)]
				obs.out = obs$intensity[which(obs$r.obs>=r)]
				p.in=mean(obs.in)
				p.out=mean(obs.out)
				log_lhood = sum(obs.in)*p.in+sum(1-obs.in)*(1-p.in)+sum(obs.out)*p.out+sum(1-obs.out)*(1-p.out)
				if(any(ordering=="I" & p.in>p.out,ordering=="O" & p.in<p.out,ordering=="N")){
					return(log_lhood)
				}else {
					return(-Inf)
				}
			}
			eval_seq = seq(from = min(obs$r.obs)+0.02, to = min(.5,max(obs$r.obs)-0.02), by=0.01 )
			len=length(eval_seq)
			evals=rep(0,len)
			for(i in 1:len){evals[i]=ini.mean.estimator(eval_seq[i])}
			inimean = eval_seq[which.max(evals)] 
		}

		output = BayesBDbinary(obs, inimean, nrun, nburn, J, ordering, mask = rep(1,length(obs$intensity)), slice, outputAll)

		return(list(image = image, center = center, gamma.fun = gamma.fun, output = output, obs = obs))

	}else if(is.list(image)){
		
		if(is.null(mask)){
			mask = rep(1,length(image$intensity))
		}else {
			center = image$center
			new.r.obs = image$r.obs[mask==1]
			new.theta.obs = image$theta.obs[mask==1]
			new.intensity = image$intensity[mask==1]
			image = list(r.obs=new.r.obs, theta.obs = new.theta.obs, intensity=new.intensity, center=center)
		}

		if(is.null(inimean)){
			ini.mean.estimator = function(r){
				obs.in = image$intensity[which(image$r.obs<r)]
				obs.out = image$intensity[which(image$r.obs>=r)]
				p.in=mean(obs.in)
				p.out=mean(obs.out)
				log_lhood = sum(obs.in)*p.in+sum(1-obs.in)*(1-p.in)+sum(obs.out)*p.out+sum(1-obs.out)*(1-p.out)
				if(any(ordering=="I" & p.in>p.out,ordering=="O" & p.in<p.out,ordering=="N")){
					return(log_lhood)
				}else {
					return(-Inf)
				}
			}
			eval_seq = seq(from = min(image$r.obs)+0.02, to = min(.5,max(image$r.obs)-0.02), by=0.01 )
			len=length(eval_seq)
			evals=rep(0,len)
			for(i in 1:len){evals[i]=ini.mean.estimator(eval_seq[i])}
			inimean = eval_seq[which.max(evals)] 
		}

		output = BayesBDbinary(image, inimean, nrun, nburn, J, ordering, mask = rep(1,length(image$intensity)), slice, outputAll)
		return(list(image = image, center = center, gamma.fun = gamma.fun, output = output, obs = list(r.obs=as.vector(image$r.obs), theta.obs=as.vector(image$theta.obs), intensity=as.vector(image$intensity), center=center)))
	}else {
		return("Input image is not a compatible image file nor a compatible list object.")
	}
}