nx<-7
ny<-7
n.pi<-4
grad.im<-expand.grid(x=c(-nx:nx)/nx*n.pi*pi,y=c(-ny:ny)/ny*n.pi*pi)

grad.im<-cbind(grad.im,
               col=1.5*with(grad.im,abs(cos(x*y))/(abs(x*y)+(3*pi/nx)))+
                 0.5*runif(nrow(grad.im)),
               x.vec=with(grad.im,y),
               y.vec=with(grad.im,x))
# normalize vector
grad.im[,c("x.vec","y.vec")]<-with(grad.im,cbind(x.vec/(sqrt(x.vec^2+y.vec^2)),
                                         y.vec/(sqrt(x.vec^2+y.vec^2))))
bn.wid<-c(diff(range(grad.im$x.vec))/10,diff(range(grad.im$y.vec))/10)
ggplot(grad.im,aes(x=x,y=y,fill=col))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.15,"cm")),size=2)+
  scale_fill_gradient(low="black",high="white")+
  guides(fill=F)+labs(title="Orientation Map")+
  theme_bw(20)

ggplot(grad.im,aes(x=x.vec,y=y.vec))+
  stat_bin2d(binwidth=c(0.25,.25),drop=F)+
  labs(x="Pfx",y="Pfy",fill="Frequency",title="Orientation Histogram")+
  xlim(-1,1)+ylim(-1,1)+
  theme_bw(20)

g.with.thresh<-cbind(grad.im,in.thresh=with(grad.im,x.vec>0.25 & y.vec>0.25))
ggplot(g.with.thresh,
       aes(x=x,y=y,fill=in.thresh,color=in.thresh))+
  geom_tile(alpha=0.5,aes(fill=in.thresh))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.2,"cm")),size=3)+
  labs(color="In Threshold")+guides(fill=FALSE)+
  theme_bw(20)

bn.wid<-c(0.25,.25)
keep.bins<-expand.grid(x.vec=seq(-1,1,bn.wid[1]/10),y.vec=seq(-1,1,bn.wid[2]/10))
keep.bins<-cbind(keep.bins,in.thresh=with(keep.bins,x.vec>0.25 & y.vec>0.25))
ggplot(g.with.thresh,aes(x=x.vec,y=y.vec))+
  stat_bin2d(binwidth=bn.wid,drop=F)+
  geom_tile(data=subset(g.with.thresh,in.thresh),fill="red",alpha=0.4)+
  labs(x="Pfx",y="Pfy",fill="Threshold")+xlim(-1,1)+ylim(-1,1)+
  theme_bw(20)

i.vec<-c(1,0)
j.vec<-c(0,1)
g.cmp.thresh<-cbind(grad.im,in.thresh=with(grad.im,
                                            abs(x.vec*i.vec[1]+y.vec*i.vec[2])==1
                                            ))
ggplot(g.cmp.thresh,
       aes(x=x,y=y,fill=in.thresh,color=in.thresh))+
  geom_tile(alpha=0.5,aes(fill=in.thresh))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.2,"cm")),size=3)+
  labs(color="In Threshold")+guides(fill=FALSE)+
  theme_bw(20)

i.vec<-c(1,0)
j.vec<-c(0,1)
ang.accept<-function(c.ang) g.cmp.thresh<-cbind(grad.im,
                                                ang.val=c.ang,
                                                in.thresh=with(grad.im,
                                            acos(x.vec*i.vec[1]+y.vec*i.vec[2])<=c.ang/180*pi
                                            ))
ang.vals<-seq(5,180,length.out=9)
ggplot(ldply(ang.vals,ang.accept),
       aes(x=x,y=y,fill=in.thresh,color=in.thresh))+
  geom_tile(alpha=0.5,aes(fill=in.thresh))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.2,"cm")),size=3)+
  facet_wrap(~ang.val)+
  labs(color="In Threshold")+guides(fill=FALSE)

  theme_bw(20)

lout<-21
b.grid<-expand.grid(x=seq(-10,10,length.out=lout),y=seq(-10,10,length.out=lout))
ring.image<-mutate(
  b.grid,
  Intensity = abs(cos(sqrt(x^2+y^2)/15*6.28))
  )

ring.image.dist<-ddply(ring.image,.(x,y),function(in.pts) {
  in.pts$Distance<-with(in.pts,sqrt(
      min(
        x^2+y^2
      )
    ))
  in.pts
})

ring.image.seg<-mutate(ring.image.dist,
                       Segmented = ifelse(
                         Intensity>0.5 & Distance>5 & Distance<10,
                         "Foreground","Background")
                       )
ring.image.seg<-mutate(ring.image.seg,
                       x=round(x+11),
                       y=round(y+11))
ring.im<-df.to.im(mutate(
  b.grid,
  val = 127*(abs(cos(sqrt(x^2+y^2)/15*6.28))+runif(nrow(b.grid)))
  ))

ring.img<-im.to.df(ring.im,out.col="Intensity")
blur.img<-im.to.df(gblur(ring.im,0.5),out.col="Intensity")
med.img<-im.to.df(medianFilter(ring.im/255,3)*255.0,out.col="Intensity")

ggplot(ring.image.seg,aes(x,y,fill=Segmented))+
  geom_raster()+
  coord_equal()+
  theme_bw(15)

simple.image<-ggplot(ring.img,aes(x,y,fill=Intensity))+
  geom_tile()+
  coord_equal()+
  theme_bw(15)
simple.image

tr.seg.img<-
  mutate(ring.img,Segmented=ifelse(Intensity>170,"Foreground","Background"))[,
   c("x","y","Segmented")]
gr.truth.img<-ring.image.seg[,c("x","y","Segmented")]
ggplot(tr.seg.img,
       aes(x,y,fill=Segmented))+
  geom_tile(color="black")+
  coord_equal()+
  labs(fill="Threshold\nSegmentation")+
  theme_bw(15)

score.img<-mutate(
    merge(tr.seg.img,gr.truth.img,by=c("x","y"),suffixes=c("C","GT")),
  Match=(SegmentedC==SegmentedGT),
  Category=
    ifelse(Match,
           ifelse(SegmentedC=="Foreground","True Positive","True Negative"),
           ifelse(SegmentedC=="Foreground","False Positive","False Negative")
      )
  )

ggplot(score.img,
       aes(x,y,fill=SegmentedC,color=SegmentedGT))+
  geom_tile(size=1)+
  facet_wrap(~Category)+
  coord_equal()+
  scale_colour_manual(values=c("green","black"))+
  scale_fill_manual(values=c("red","blue"))+
  labs(fill="Threshold\nSegmentation",color="Ground\nTruth")+
  theme_bw(20)

mk.th.img<-function(in.img,steps,smin=0.3,smax=0.7) ldply(seq(smin,smax,length.out=steps),function(c.thresh) {
  mutate(in.img,
         Segmented = ifelse(
                         Intensity>c.thresh*255,
                         "Foreground","Background"),
         Thresh = c.thresh
        )
})

ggplot(mk.th.img(ring.img,9),aes(x,y,fill=Segmented))+
  geom_tile(color="black")+
  coord_equal()+
  facet_wrap(~Thresh)+
  theme_bw(15)

mk.stat.table<-function(in.img,steps,smin=0.3,smax=0.7) 
  ddply(mk.th.img(in.img,steps,smin,smax),.(Thresh),function(c.img) {
  in.df<-merge(c.img[,c("x","y","Segmented")],
               ring.image.seg[,c("x","y","Segmented")],
               by=c("x","y"),suffixes=c("C","GT"))
  in.df<-mutate(in.df,Match=(SegmentedC==SegmentedGT))
  data.frame(TP=nrow(subset(in.df,SegmentedGT=="Foreground" & Match)),
             TN=nrow(subset(in.df,SegmentedGT=="Background" & Match)),
             FP=nrow(subset(in.df,SegmentedC=="Foreground" & !Match)),
             FN=nrow(subset(in.df,SegmentedC=="Background" & !Match))
             )
})
tp.table<-mk.stat.table(ring.img,6,smin=0,smax=1)
kable(tp.table)

ggplot(melt(tp.table,id.vars=c("Thresh")),
       aes(x=Thresh*100,y=100*value/200,color=variable))+
  geom_point()+
  geom_path()+
  labs(color="",y="Percent (%)",x="Threshold Value (%)")+
  theme_bw(20)

prec.table<-mutate(mk.stat.table(ring.img,6),
                  Recall=round(100*TP/(TP+FN)),
                  Precision=round(100*TP/(TP+FP))
                   )
kable(prec.table)

prec.table<-mutate(mk.stat.table(ring.img,100),
                   Recall=TP/(TP+FN),
                   Precision=TP/(TP+FP)
                   )

ggplot(prec.table,aes(x=Recall,y=Precision))+
  geom_point(aes(color=Thresh,shape="Image"),size=5)+
  geom_point(data=data.frame(Recall=1,Precision=1),
             aes(shape="Ideal"),size=5,color="red")+
  geom_path()+
  #coord_equal()+
  scale_color_gradientn(colours=rainbow(6))+
  labs(shape="")+
  theme_bw(20)

all.img<-rbind(
  cbind(mk.th.img(ring.img,3), name="Normal"),
  cbind(mk.th.img(blur.img,3), name="Gaussian"),
  cbind(mk.th.img(med.img,3), name="Median")
  )

simple.image<-ggplot(all.img,aes(x,y,fill=Segmented))+
  geom_tile(color="black")+
  coord_equal()+
  facet_grid(name~Thresh)+
  theme_bw(15)

simple.image

stab<-function(img) mk.stat.table(img,100,smin=0.05,smax=0.95)
all.table<-rbind(
  cbind(stab(ring.img), name="Normal"),
  cbind(stab(blur.img), name="Gaussian"),
  cbind(stab(med.img), name="Median")
  )

prec.table<-mutate(all.table,
                   Recall=TP/(TP+FN),
                  Precision=TP/(TP+FP)
                   )
ggplot(prec.table,aes(x=Recall,y=Precision,color=name))+
  geom_point(aes(shape="Image"),size=5)+
  geom_point(data=data.frame(Recall=1,Precision=1),
             aes(shape="Ideal"),size=5,color="red")+
  geom_path()+
  #coord_equal()+
  labs(shape="",color="Filter")+
  xlim(0.5,1)+ylim(0.5,1)+
  theme_bw(20)

tp.table<-mutate(all.table,
                   TPRate=TP/(TP+FN),
                   FPRate=FP/(FP+TN)
                   )
ggplot(tp.table,aes(x=FPRate,y=TPRate,color=name))+
  geom_point(aes(shape="Image"),size=2)+
  geom_point(data=data.frame(TPRate=1,FPRate=0),
             aes(shape="Ideal"),size=5,color="red")+
  geom_path()+
  geom_segment(x=0,xend=1,y=0,yend=1,aes(color="Random Guess"))+
  xlim(0,1)+ylim(0,1)+coord_equal()+
  labs(shape="",x="False Positive Rate",y="True Positive Rate",linetype="")+
  theme_bw(20)

# image data
road.rgb.arr<-readJPEG(qbi.data("road_image.jpg"))
road.arr<-((road.rgb.arr[,,1]+road.rgb.arr[,,2]+road.rgb.arr[,,3])/3) %>% t %>% flip

## Error in asub.default(x, seq.int(dim(x)[2], 1), 2): could not find function "Quote"

road.image<- road.arr %>% im.to.df

## Error in eval(expr, envir, enclos): object 'road.arr' not found

# ground truth
cal.image<-readJPEG(qbi.data("road_image_street_bw.jpg")) %>% t %>% flip %>% im.to.df %>% 
  mutate(Segmented = ifelse(val,"Road","Background"))

## Error in asub.default(x, seq.int(dim(x)[2], 1), 2): could not find function "Quote"

ggplot(road.image,
       aes(x,y,fill=val))+
  geom_raster()+
  coord_equal()+
  scale_fill_gradientn(colours=c("black","white"))+
  guides(fill=F)+
  theme_bw(20)

ggplot(cal.image,aes(x,y,fill=Segmented,alpha=Segmented))+
  geom_raster()+
  coord_equal()+
  guides(alpha=F)+
  theme_bw(20)

## Error in ggplot(cal.image, aes(x, y, fill = Segmented, alpha = Segmented)): object 'cal.image' not found

tr.road<-mutate(road.image,
                 Segmented=ifelse(val>0.75,"Road","Background")
                 )

## Error in mutate_(.data, .dots = lazyeval::lazy_dots(...)): object 'road.image' not found

ggplot(tr.road,aes(x,y,fill=Segmented,alpha=val))+
  geom_raster()+
  coord_equal()+
  guides(alpha=F)+
  theme_bw(20)

## Error in ggplot(tr.road, aes(x, y, fill = Segmented, alpha = val)): object 'tr.road' not found

score.img<-mutate(
    merge(tr.road[,c("x","y","Segmented")],
          cal.image[,c("x","y","Segmented")],
          by=c("x","y"),suffixes=c("C","GT")),
  Match=(SegmentedC==SegmentedGT),
  Category=
    ifelse(Match,
           ifelse(SegmentedC=="Road","True Positive","True Negative"),
           ifelse(SegmentedC=="Road","False Positive","False Negative")
      )
  )

## Error in merge(tr.road[, c("x", "y", "Segmented")], cal.image[, c("x", : object 'tr.road' not found

ggplot(score.img,
       aes(x,y,fill=SegmentedGT))+
  geom_tile()+
  facet_grid(Category~.)+
  coord_equal()+
  scale_colour_manual(values=c("green","black"))+
  scale_fill_manual(values=c("red","blue"))+
  labs(fill="Ground\nTruth")+
  theme_bw(20)

mk.th.img<-function(in.img,stepseq) ldply(stepseq,function(c.thresh) {
  mutate(in.img,
         Segmented = ifelse(
                         val>c.thresh,
                         "Road","Background"),
         Thresh = c.thresh
        )
})

mk.stat.table<-function(in.img,stepseq) 
  ddply(mk.th.img(in.img,stepseq),.(Thresh),function(c.img) {
  in.df<-merge(c.img[,c("x","y","Segmented")],
               cal.image[,c("x","y","Segmented")],
               by=c("x","y"),suffixes=c("C","GT"))
  in.df<-mutate(in.df,Match=(SegmentedC==SegmentedGT))
  data.frame(TP=nrow(subset(in.df,SegmentedGT=="Road" & Match)),
             TN=nrow(subset(in.df,SegmentedGT=="Background" & Match)),
             FP=nrow(subset(in.df,SegmentedC=="Road" & !Match)),
             FN=nrow(subset(in.df,SegmentedC=="Background" & !Match))
             )
})

step.fcn<-function(vals,n) quantile(vals,seq(0,1,length.out=n+2)[-c(1,n+2)])
abs.steps<-step.fcn(road.image$val,6)

## Error in quantile(vals, seq(0, 1, length.out = n + 2)[-c(1, n + 2)]): object 'road.image' not found

tp.table<-mk.stat.table(road.image,abs.steps)

## Error in inherits(.data, "split"): object 'abs.steps' not found

prec.table<-mutate(tp.table,
                   Thresh=round(100*Thresh),
                   Recall=round(100*TP/(TP+FN)),
                   Precision=round(100*TP/(TP+FP))
                   )
kable(prec.table)

abs.steps<-step.fcn(road.image$val,15)

## Error in quantile(vals, seq(0, 1, length.out = n + 2)[-c(1, n + 2)]): object 'road.image' not found

tp.table<-mk.stat.table(road.image,abs.steps)

## Error in inherits(.data, "split"): object 'abs.steps' not found

ggplot(melt(tp.table,id.vars=c("Thresh")),
       aes(x=Thresh*100,y=value,color=variable))+
  geom_point()+
  geom_path()+
  labs(color="",y="Pixel Count",x="Threshold Value (%)")+
  scale_y_sqrt()+
  theme_bw(20)

## Error in self$trans$transform(x): non-numeric argument to mathematical function

prec.table<-mutate(tp.table,
                   Recall=TP/(TP+FN),
                  Precision=TP/(TP+FP)
                   )

ggplot(prec.table,aes(x=Recall,y=Precision))+
  geom_point(aes(color=Thresh,shape="Image"),size=5)+
  geom_point(data=data.frame(Recall=1,Precision=1),
             aes(shape="Ideal"),size=5,color="red")+
  geom_path()+
  scale_color_gradientn(colours=rainbow(6))+
  scale_x_sqrt()+
  labs(shape="")+
  theme_bw(20)

tpr.table<-mutate(tp.table,
                   TPRate=TP/(TP+FN),
                   FPRate=FP/(FP+TN)
                   )
ggplot(tpr.table,aes(x=FPRate,y=TPRate))+
  geom_point(aes(shape="Image"),size=2)+
  geom_point(data=data.frame(TPRate=1,FPRate=0),
             aes(shape="Ideal"),size=5,color="red")+
  geom_path()+
  geom_segment(x=0,xend=1,y=0,yend=1,aes(color="Random Guess"))+
  xlim(0,1)+ylim(0,1)+coord_equal()+
  labs(shape="",x="False Positive Rate",y="True Positive Rate",linetype="")+
  theme_bw(20)

# subtract blur
#road.blur <- (road.arr - gblur(road.arr,2)) %>% im.to.df
roi.road.image<-road.image %>% 
  mutate(val = 
           ifelse(x>200,ifelse(y<200,val,0),0)
         )

## Error in eval(expr, envir, enclos): object 'road.image' not found

dci.steps<-step.fcn(roi.road.image$val,12)

## Error in quantile(vals, seq(0, 1, length.out = n + 2)[-c(1, n + 2)]): object 'roi.road.image' not found

dci.table<-mk.stat.table(roi.road.image,dci.steps)

## Error in inherits(.data, "split"): object 'dci.steps' not found

prec.dci.table<-mutate(dci.table,
                   Recall=TP/(TP+FN),
                  Precision=TP/(TP+FP)
                   )

## Error in mutate_(.data, .dots = lazyeval::lazy_dots(...)): object 'dci.table' not found

ggplot(roi.road.image,
       aes(x,y,fill=val))+
  geom_tile()+
  coord_equal()+
  scale_fill_gradientn(colours=c("black","white"))+
  guides(fill=F)+
  theme_bw(20)

## Error in ggplot(roi.road.image, aes(x, y, fill = val)): object 'roi.road.image' not found

# subtract blur
# road.blur <- (road.arr - gblur(road.arr,2)) %>% im.to.df
roi.blur.image<-gblur(road.arr,2) %>% im.to.df %>% 
  mutate(val = 
           ifelse(x>200,ifelse(y<200,val,0),0)
         )

## Error in validObject(x): object 'road.arr' not found

blur.steps<-step.fcn(roi.blur.image$val,12)

## Error in quantile(vals, seq(0, 1, length.out = n + 2)[-c(1, n + 2)]): object 'roi.blur.image' not found

blur.table<-mk.stat.table(roi.blur.image,blur.steps)

## Error in inherits(.data, "split"): object 'blur.steps' not found

prec.blur.table<-mutate(blur.table,
                   Recall=TP/(TP+FN),
                  Precision=TP/(TP+FP)
                   )

## Error in mutate_(.data, .dots = lazyeval::lazy_dots(...)): object 'blur.table' not found

ggplot(roi.blur.image,
       aes(x,y,fill=val))+
  geom_tile()+
  coord_equal()+
  scale_fill_gradientn(colours=c("black","white"))+
  guides(fill=F)+
  theme_bw(20)

## Error in ggplot(roi.blur.image, aes(x, y, fill = val)): object 'roi.blur.image' not found

ggplot(rbind(cbind(prec.table,Channel="Standard"),
             cbind(prec.dci.table,Channel="Region"),
             cbind(prec.blur.table,Channel="Region Blur")
             ),
             aes(x=Recall,y=Precision))+
  geom_point(aes(shape="Image"),size=5)+
  geom_point(data=data.frame(Recall=1,Precision=1),
             aes(shape="Ideal"),size=5,color="red")+
  geom_path(aes(group=Channel,color=Channel))+
  scale_y_sqrt()+
  labs(shape="")+
  theme_bw(20)

## Error in cbind(prec.dci.table, Channel = "Region"): object 'prec.dci.table' not found

nx<-4
ny<-4
n.pi<-4
grad.im<-expand.grid(x=c(-nx:nx)/nx*n.pi*pi,
                     y=c(-ny:ny)/ny*n.pi*pi)

grad.im<-cbind(grad.im,
               col=1.5*with(grad.im,abs(cos(x*y))/(abs(x*y)+(3*pi/nx)))+
                 0.5*runif(nrow(grad.im)),
              a=with(grad.im,sqrt(2/(abs(x)+0.5))),
               b=with(grad.im,0.5*sqrt(abs(x)+1)),
              th=0.5*runif(nrow(grad.im)),
              aiso=1,count=1)

create.ellipse.points<-function(x.off=0,y.off=0,a=1,b=NULL,th.off=0,th.max=2*pi,pts=36,...) {
  if (is.null(b)) b<-a
  th<-seq(0,th.max,length.out=pts)
  data.frame(x=a*cos(th.off)*cos(th)+b*sin(th.off)*sin(th)+x.off,
             y=-1*a*sin(th.off)*cos(th)+b*cos(th.off)*sin(th)+y.off,
             id=as.factor(paste(x.off,y.off,a,b,th.off,pts,sep=":")),...)
}
deform.ellipse.draw<-function(c.box) {
  create.ellipse.points(x.off=c.box$x[1],
                        y.off=c.box$y[1],
                        a=c.box$a[1],
                        b=c.box$b[1],
                        th.off=c.box$th[1],
                        col=c.box$col[1])                    
}

# normalize vector
tens.im<-ddply(grad.im,.(x,y),deform.ellipse.draw)

ggplot(tens.im,aes(x=x,y=y,group=as.factor(id),fill=col))+
  geom_polygon(color="black")+coord_fixed(ratio=1)+scale_fill_gradient(low="black",high="white")+guides(fill=F)+
  theme_bw(20)

ggplot(grad.im)+
  geom_histogram(aes(x=a,fill="Width"),alpha=0.7)+
  geom_histogram(aes(x=b,fill="Height"),alpha=0.7)+
  geom_histogram(aes(x=th,fill="Orientation"),alpha=0.7)+
  geom_histogram(aes(x=col,fill="Color"),alpha=0.7)+
  guides(color=F)+labs(fill="Variable")+
  theme_bw(15)

plot(grad.im[,c("a","b","th","col")])

ggplot(shale.vals,aes(x=x,y=y,fill=th.label,alpha=1-val))+coord_fixed(ratio=1)+
  geom_raster()+guides(fill=F,alpha=F)+theme_bw(20)+facet_wrap(~th.label)

make.circle<-function(nx,r,nscale=1,
                      circ.fun=function(x,y,r) ((x^2+y^2)