setwd("D:/Cocorobo/test/RobustVideoMatting/test/")

library(jpeg)
library(plotrix)
xstart=0; xend=90;
ystart=0; yend=55;
plot(NA,NA,xaxt='n', yaxt='n', xlim=c(xstart,xend),ylim=c(ystart,yend),xlab=NA, ylab=NA)
final <- c()
for (pic_input in 1001:1030){
  top_ratio = 0.2; bottom_ratio = 0.7
  #####第1步，打开指定图片，并进行裁剪
  {
    #####模拟化JPG图片，横坐标150，纵坐标按比例
    pic1 <- readJPEG(paste("02Shape/",pic_input,".jpg",sep=""))
    # pic1 <- pic1[,1:1400,]
    top = floor(top_ratio*nrow(pic1)); bottom = ceiling(bottom_ratio*nrow(pic1))
    pic2 <- pic1[,,1]*0.299 + pic1[,,2]*0.587 + pic1[,,3]*0.114;
    pic2[pic2 > 0.985] = 1
    pic3 <- (1-pic2)
    tmp = pic3[bottom,]; tmp[tmp > 0.01] = 1; tmp[tmp <= 0.01] = 0
    names(tmp) = 1:ncol(pic3)
    left = max(1,min(as.numeric(names(tmp[tmp == 1]))) - 100)
    right = min(max(as.numeric(names(tmp[tmp == 1]))) + 100,ncol(pic1))
    # left = as.numeric(names(tmp[tmp > 1][1]))
    # right = as.numeric(names(tmp[tmp > 1][length(tmp[tmp > 1])]))
    pic1 <- pic1[top:bottom,left:right,]
    wi = ceiling(1:ncol(pic1)/(ncol(pic1)/200))
    hi = ceiling(1:nrow(pic1)/(ncol(pic1)/200))
    npic1 <- array(NA, dim = c(max(hi),200,3))
    for (i in 1:200){
      for (j in 1:max(hi)){
        npic1[j,i,1] = mean(pic1[which(hi == j),which(wi == i),1])
        npic1[j,i,2] = mean(pic1[which(hi == j),which(wi == i),2])
        npic1[j,i,3] = mean(pic1[which(hi == j),which(wi == i),3])
      }
    }
    tmp = array(1,dim = c(nrow(npic1)+6,6+ncol(npic1),3))
    tmp[4:(3+nrow(npic1)),4:(3+ncol(npic1)),] <- npic1
    npic1 <- tmp;
    
    #####模拟化黑白图片
    pic2 <- pic1[,,1]*0.299 + pic1[,,2]*0.587 + pic1[,,3]*0.114
    npic2 <- array(NA, dim = c(max(hi),200))
    for (i in 1:200){
      for (j in 1:max(hi)){
        npic2[j,i] = mean(pic2[which(hi == j),which(wi == i)])
      }
    }
    tmp = array(1, dim = c(nrow(npic2)+6,6+ncol(npic2)))
    tmp[4:(3+nrow(npic2)),4:(3+ncol(npic2))] <- npic2
    npic2 <- tmp
    
    rm(tmp,wi,hi,i,j,left,right,pic2)
    print(paste("第",pic_input-1000,"帧，分析中..",sep=""))
  }
  #####第2步，黑白化并进行边缘化修整
  {
    #####黑白二维
    npic3 <- 1-npic2;
    npic3[npic3 < 0.01] = 0
    npic3[npic3 >= 0.01] = 1
    npic2 <- npic2*npic3
    print(paste("第",pic_input-1000,"帧，分析中...",sep=""))
  }
  ####第3步:确定平滑边缘
  {
    ####第三步，两轮裁剪
    {
      input = npic3
      for (n_round in 1:5){
        ####rounding
        {
          tmp <- as.data.frame(input)
          
          dim(tmp)
          ####left
          {
            picl <- cbind(tmp, rep(0,nrow(input)))
            picl <- picl[,2:ncol(picl)]
            dim(picl)
            picl <- as.matrix(picl)
            picll <- abs(input - picl)
          }
          ####right
          {
            picr <- cbind(tmp,rep(0,nrow(input)))
            picr <- picr[,c(ncol(picr),1:(ncol(picr)-2))]
            dim(picr)
            picr <- as.matrix(picr)
            picrr <- abs(input - picr)
          }
          ###bottom
          {
            picb <- rbind(tmp,rep(1,ncol(input)))
            picb <- picb[2:nrow(picb),]
            dim(picb)
            picb <- as.matrix(picb)
            picbb <- abs(input - picb)
          }
          ###top
          {
            pict <- rbind(tmp,rep(0,ncol(input)))
            pict <- pict[c(nrow(pict),1:(nrow(pict)-2)),]
            dim(pict)
            pict <- as.matrix(pict)
            pictt <- abs(input - pict)
          }
          tmp = picll+picrr+picbb+pictt
          tmp[tmp == 0] = 5
          tmp[tmp == 2] = 0
          tmp[tmp == 3] = 0
          tmp[tmp == 4] = 0
          tmp[tmp == 5] = 1
          tmp2 = input*tmp
          input <- tmp2
        }
        output2 <- tmp2
        npic1[,,1] <- npic1[,,1]*output2;
        npic1[,,2] <- npic1[,,2]*output2;
        npic1[,,3] <- npic1[,,3]*output2;
        npic3 <- npic3*output2
        ####确定surface
        {
          bk <- output2
          tmp = bk[nrow(bk)-3,]; names(tmp) = 1:ncol(bk)
          x_start = 1; tmp2 = as.numeric(names(tmp[tmp == 1]));
          for (i in 1:length(tmp2)){
            if (sum((bk[nrow(bk)-3,(tmp2[i]):(tmp2[i]+29)])) == 30){
              x_start = tmp2[i]; break;
            }
          }
          y_start = nrow(bk)-3
          surface <- data.frame()
          surface <- rbind(surface, c(x_start, y_start))
          colnames(surface) <- c("xpos","ypos")
          x_input = x_start; y_input = y_start;
          
          bk2 <- cbind(bk, rep(0,nrow(bk)))
          bk2 <- bk2[,c(ncol(bk2),1:ncol(bk2))]
          bk2 <- rbind(bk2, rep(1,ncol(bk2)))
          bk2 <- rbind(bk2, rep(0,ncol(bk2)))
          bk2 <- bk2[c(nrow(bk2),1:(nrow(bk2)-1)),]
          
          for (num in 1:10000){
            pos <- data.frame(xpos = x_input + c(-1,-1,-1,0,0,1,1,1),
                              ypos = y_input + c(1,0,-1,1,-1,1,0,-1))
            pos <- pos[(pos$xpos >= 1) & (pos$xpos <= ncol(bk)) & (pos$ypos >= 1) & (pos$ypos <= nrow(bk)),]
            pos <- pos[!(paste(pos$xpos,pos$ypos) %in% paste(surface$xpos,surface$ypos)),]
            target <- data.frame()
            for (i in 1:nrow(pos)){
              if (bk[pos$ypos[i],pos$xpos[i]] == 1){
                tpos <- data.frame(xpos = pos$xpos[i] + c(-1,0,1,0),
                                   ypos = pos$ypos[i] + c(0,-1,0,1))
                # tpos <- data.frame(xpos = pos$xpos[i] + c(-1,-1,-1,0,0,1,1,1),
                #                    ypos = pos$ypos[i] + c(1,0,-1,1,-1,1,0,-1))
                va <- c()
                for (j in 1:nrow(tpos)){
                  va = c(va, bk2[tpos$ypos[j]+1, tpos$xpos[j]+ 1])
                }
                if (length(va[va==0])>0){
                  target <- rbind(target, c(pos$xpos[i],pos$ypos[i]))
                }
              }
              
            }
            if (nrow(target) < 1){ break; }
            colnames(target) = c("xpos","ypos")
            ####最短距离
            distance = sqrt((target$xpos - x_input)**2 + (target$ypos - y_input)**2)
            surface <- rbind(surface, c(target$xpos[order(distance)[1]],
                                        target$ypos[order(distance)[1]]))
            x_input = target$xpos[order(distance)[1]];
            y_input = target$ypos[order(distance)[1]];
          }
          if (length(surface$ypos[1:20][surface$ypos[1:20] == (nrow(npic3) - 3)]) > 10){
            surface <- data.frame(xpos = rev(surface$xpos), ypos = rev(surface$ypos))
          }
          surface2 <- surface[surface$ypos != (nrow(npic3) - 3),]
          judge1 = nrow(surface2[(surface2$xpos < 100) & (surface2$ypos < nrow(npic3)) & (surface2$ypos > (nrow(npic3)-20)),])
          judge2 = nrow(surface2[(surface2$xpos > 100) & (surface2$ypos < nrow(npic3)) & (surface2$ypos > (nrow(npic3)-20)),])
        }
        xmiddle = mean(surface2$xpos)
        if ((judge1 > 10) & (judge2 > 10)){
          break;
        }
        input <- output2
      }
    }
    ######
    {
      #####
      surface_left = data.frame()
      for (j in 1:nrow(surface)){
        if (surface$xpos[j] > xmiddle){break;}
        surface_left <- rbind(surface_left, c(surface$xpos[j], surface$ypos[j]))
      }
      #####
      surface_right <- data.frame()
      surface_right <- surface[(j+1):nrow(surface),]
    }
    
    if ((judge1 > 10) & (judge2 > 10) & (nrow(surface_left) > 100) & (nrow(surface_right) > 100)){
      colnames(surface_left) = c("xpos","ypos")
      surface_left <- data.frame(xpos = rev(surface_left$xpos), ypos = rev(surface_left$ypos))
      colnames(surface_right) = c("xpos","ypos")
      #####第五步，计算差异点
      {
        tmp1 = array(1,dim = c(nrow(npic1),ncol(npic1)))
        tmp1[1:(nrow(tmp1)-1),] <- npic1[2:nrow(npic1),,1]
        tmp2 = array(1,dim = c(nrow(npic1),ncol(npic1)))
        tmp2[1:(nrow(tmp2)-1),] <- npic1[2:nrow(npic1),,2]
        tmp3 = array(1,dim = c(nrow(npic1),ncol(npic1)))
        tmp3[1:(nrow(tmp3)-1),] <- npic1[2:nrow(npic1),,3]
        
        tmp = abs(npic1[,,1] - tmp1) + abs(npic1[,,2] - tmp2) + abs(npic1[,,3] - tmp3)
        tmp[tmp >= 0.15] = 1
        tmp[tmp < 1] = 0
        ###去掉边缘点
        target <- data.frame()
        for (i in 3:(nrow(tmp)-3)){
          for (j in 3:(ncol(tmp)-3)){
            if (tmp[i,j] == 1){
              tmp2 = c(npic3[i-1,j+1],npic3[i-1,j],npic3[i-1,j-1],
                       npic3[i,j+1],npic3[i,j-1],
                       npic3[i+1,j+1],npic3[i+1,j],npic3[i+1,j-1]);
              if (length(tmp2[tmp2 == 0]) == 0){
                target <- rbind(target, c(i,j))
              }
            }
          }
        }
        colnames(target) <- c("ypos","xpos")
        test = table(target$ypos)
        candidate <- as.numeric(names(test[test >= 10]))
        rm(i,j,test,tmp,tmp1,tmp2,tmp3)
        print(paste("第",pic_input-1000,"帧，分析中.....",sep=""))
      }
      
      
      if (length(candidate) > 0){    #####第六步，找到脖子曲线对应的边缘曲线
        {
          #####
          pointl <- c(); pointr <- c()
          for (i in 1:length(candidate)){
            for (j in 1:nrow(surface_left)){
              if (surface_left$ypos[j] >= candidate[i]){ break; }
            }
            partl <- surface_left[j:nrow(surface_left),]
            if (nrow(partl) > 40){
              dis1 = (partl$xpos[5] - partl$xpos[10:40]);
              dis2 = sqrt((partl$xpos[5] - partl$xpos[10:40])**2 + (partl$ypos[5] - partl$ypos[10:40])**2);
              angle1 = mean(180*acos(dis1/dis2)/pi)
              if ((angle1 < 70) & (candidate[i] >= min(surface$ypos))){
                pointl <- c(pointl, candidate[i])
              }
            }
            for (j in 1:nrow(surface_right)){
              if (surface_right$ypos[j] >= candidate[i]){ break; }
            }
            partr <- surface_right[j:nrow(surface_right),]
            if (nrow(partr) > 40){
              dis1 = (partr$xpos[5] - partr$xpos[10:40]);
              dis2 = sqrt((partr$xpos[5] - partr$xpos[10:40])**2 + (partr$ypos[1] - partr$ypos[10:40])**2);
              angle1 = mean(180*acos(dis1/dis2)/pi)
              if ((angle1 > 110) & (candidate[i] >= min(surface$ypos))){
                pointr <- c(pointr, candidate[i])
              }
            }
          }
          points <- intersect(pointl, pointr)
          if (length(points) > 1){
            points <- points[points < (min(surface$ypos)+100)]
            points <- sort(points, decreasing = T)
            point_tmp <- data.frame()
            for (i in 1:(length(points)-1)){
              for (j in 1:length(points)){
                if ((i + j) > length(points)){break;}
                if ((points[i] - j) != points[i + j]){ break;}
              }
              point_tmp <- rbind(point_tmp, c(points[i],j))
            }
            colnames(point_tmp) <- c("v1","v2")
            if (nrow(point_tmp[(point_tmp$v2 > 1),]) > 1){
              bottom = point_tmp[(point_tmp$v2 > 1),]$v1[1];
              top = point_tmp[(point_tmp$v2 > 1),]$v1[1] - point_tmp[(point_tmp$v2 > 1),]$v2[1] + 1;
            } else{
              top = point_tmp$v1[1];
              bottom = point_tmp$v1[1];
            }
            
            #####partl
            {
              for (j in 1:nrow(surface_left)){
                if (surface_left$ypos[j] >= top){ break; }
              }
              partl <- surface_left[j:nrow(surface_left),]
            }
            #####partr
            {
              for (j in 1:nrow(surface_right)){
                if (surface_right$ypos[j] >= top){ break; }
              }
              partr <- surface_right[j:nrow(surface_right),]
            }
            #####
            dis1 = xmiddle - min(partl[(partl$ypos <= floor(top + 30)),]$xpos)
            dis2 = max(partr[(partr$ypos <= floor(top + 30)),]$xpos) - xmiddle;
            leftmin = xmiddle - floor(0.8*min(dis1,dis2));
            rightmax = xmiddle + floor(0.8*min(dis1,dis2));
            
            #####partl2
            {
              for (j in 1:nrow(partl)){
                if (partl$xpos[j] < leftmin){ break; }
              }
              partl2 <- partl[1:j,]
            }
            #####partr2
            {
              for (j in 1:nrow(partr)){
                if (partr$xpos[j] > rightmax){ break; }
              }
              partr2 <- partr[1:j,]
            }
            
            #############
            target2 <- target[(target$xpos > partl$xpos[1]) & (target$xpos < partr$xpos[1]) & (target$ypos >= top) & (target$ypos <= bottom),]
            xpos <- as.numeric(names(table(target2$xpos)))
            ypos <- c()
            for (i in 1:length(xpos)){
              ypos <- c(ypos, mean(target2$ypos[which(target2$xpos == xpos[i])]))
            }
            xposl <- as.numeric(names(table(partl2$xpos)))
            yposl <- c()
            for (i in 1:length(xposl)){
              yposl <- c(yposl, mean(partl2$ypos[which(partl2$xpos == xposl[i])]))
            }
            #####
            xposr <- as.numeric(names(table(partr2$xpos)))
            yposr <- c()
            for (i in 1:length(xposr)){
              yposr <- c(yposr, mean(partr2$ypos[which(partr2$xpos == xposr[i])]))
            }
            
            ###############
            xpos2 <- c(xposl,xpos,xposr); ypos2 <- c(yposl, ypos, yposr)
            input <- data.frame(xpos = xpos2, ypos = ypos2)
            loessMod10 <- loess(ypos ~ xpos, data=input, span=0.3)
            
            smoothedleft <- predict(loessMod10, newdata = (xmiddle - 10:floor(0.8*min(dis1,dis2))))
            smoothedright<- predict(loessMod10, newdata = (xmiddle + 10:floor(0.8*min(dis1,dis2))))
            RES <- data.frame(dis = 2*(10:floor(0.8*min(dis1,dis2))),
                              left = smoothedleft, right = smoothedright)
            
            angle = 0;
            for (i in 1:nrow(RES)){
              if ((!is.na(RES$left[i])) & (!is.na(RES$right[i]))){
                ang = atan(abs(RES$left[i] - RES$right[i])/RES$dis[i])*180/pi
                # print(ang)
                if (ang > angle){ angle = ang}
              }
              
            }
          } else{
            #####
            for (j in 1:nrow(surface_left)){
              if (surface_left$ypos[j] > (min(surface$ypos) +30)){ break; }
            }
            partl <- surface_left[1:j,]
            dis1 = xmiddle - partl$xpos[j]
            #####
            for (j in 1:nrow(surface_right)){
              if (surface_right$ypos[j] > (min(surface$ypos) +30)){ break; }
            }
            partr <- surface_right[1:j,]
            
            dis2 = partr$xpos[j] - xmiddle;
            leftmin = xmiddle - floor(0.8*min(dis1,dis2));
            rightmax = xmiddle + floor(0.8*min(dis1,dis2));
            
            ####
            xposl <- as.numeric(names(table(partl$xpos)))
            yposl <- c()
            for (i in 1:length(xposl)){
              yposl <- c(yposl, mean(partl$ypos[which(partl$xpos == xposl[i])]))
            }
            #####
            xposr <- as.numeric(names(table(partr$xpos)))
            yposr <- c()
            for (i in 1:length(xposr)){
              yposr <- c(yposr, mean(partr$ypos[which(partr$xpos == xposr[i])]))
            }
            ###############
            xpos2 <- c(xposl,xposr); ypos2 <- c(yposl, yposr)
            input <- data.frame(xpos = xpos2, ypos = ypos2)
            loessMod10 <- loess(ypos ~ xpos, data=input, span=0.3)
            
            smoothedleft <- predict(loessMod10, newdata = (xmiddle - 10:floor(0.8*min(dis1,dis2))))
            smoothedright<- predict(loessMod10, newdata = (xmiddle + 10:floor(0.8*min(dis1,dis2))))
            RES <- data.frame(dis = 2*(10:floor(0.8*min(dis1,dis2))),
                              left = smoothedleft, right = smoothedright)
            
            angle = 0;
            for (i in 1:nrow(RES)){
              if ((!is.na(RES$left[i])) & (!is.na(RES$right[i]))){
                ang = atan(abs(RES$left[i] - RES$right[i])/RES$dis[i])*180/pi
                # print(ang)
                if (ang > angle){ angle = ang}
              }
              
            }
          }
          
        }} else{ 
          #####
          for (j in 1:nrow(surface_left)){
            if (surface_left$ypos[j] > (min(surface$ypos) +30)){ break; }
          }
          partl <- surface_left[1:j,]
          dis1 = xmiddle - partl$xpos[j]
          #####
          for (j in 1:nrow(surface_right)){
            if (surface_right$ypos[j] > (min(surface$ypos) +30)){ break; }
          }
          partr <- surface_right[1:j,]
          
          dis2 = partr$xpos[j] - xmiddle;
          leftmin = xmiddle - floor(0.8*min(dis1,dis2));
          rightmax = xmiddle + floor(0.8*min(dis1,dis2));
          
          ####
          xposl <- as.numeric(names(table(partl$xpos)))
          yposl <- c()
          for (i in 1:length(xposl)){
            yposl <- c(yposl, mean(partl$ypos[which(partl$xpos == xposl[i])]))
          }
          #####
          xposr <- as.numeric(names(table(partr$xpos)))
          yposr <- c()
          for (i in 1:length(xposr)){
            yposr <- c(yposr, mean(partr$ypos[which(partr$xpos == xposr[i])]))
          }
          ###############
          xpos2 <- c(xposl,xposr); ypos2 <- c(yposl, yposr)
          input <- data.frame(xpos = xpos2, ypos = ypos2)
          loessMod10 <- loess(ypos ~ xpos, data=input, span=0.3)
          
          smoothedleft <- predict(loessMod10, newdata = (xmiddle - 10:floor(0.8*min(dis1,dis2))))
          smoothedright<- predict(loessMod10, newdata = (xmiddle + 10:floor(0.8*min(dis1,dis2))))
          RES <- data.frame(dis = 2*(10:floor(0.8*min(dis1,dis2))),
                            left = smoothedleft, right = smoothedright)
          
          angle = 0;
          for (i in 1:nrow(RES)){
            if ((!is.na(RES$left[i])) & (!is.na(RES$right[i]))){
              ang = atan(abs(RES$left[i] - RES$right[i])/RES$dis[i])*180/pi
              # print(ang)
              if (ang > angle){ angle = ang}
            }
            
          }
        }
      
    } else{
      angle = 0
    }
    
    
  }
  
  #####第七步，输出角度
  {
    
    print(paste("第",pic_input-1000,"帧，偏离角度为",angle,sep=""))
  }
  #####输出波形图
  {
    dev.off()
    xs = 0; xe = 100; width = xe - xs;
    ys = 0; ye = (xe-xs)*(dim(npic1)[1]/dim(npic1)[2]); height = ye - ys;
    
    xstart=0; xend=width;
    ystart=0; yend=height;
    
    plot(NA,NA,xaxt='n', yaxt='n', xlim=c(xstart,xend),ylim=c(ystart,yend),xlab=NA, ylab=NA)
    
    rasterImage(npic2, xleft = 0, xright = 100, ytop = ye, ybottom = ys)
    for (i in 1:nrow(surface)){
      draw.ellipse(x = width*surface$xpos[i]/ncol(npic3),
                   y = ye - height*surface$ypos[i]/nrow(npic3), a = 0.2, b = 0.2, lwd = 1, border = "cyan", col = "cyan")
      
    }
    
    lines( ye - height*smoothedleft/nrow(npic3), x=width*(xmiddle - 10:floor(0.8*min(dis1,dis2)))/ncol(npic3), col="purple", lwd = 10)
    lines( ye - height*smoothedright/nrow(npic3), x=width*(xmiddle + 10:floor(0.8*min(dis1,dis2)))/ncol(npic3), col="purple", lwd = 10)
    segments(x0 = width*xmiddle/ncol(npic3), y0 = ys, y1 = ye, lwd = 10, col = "white")
    text(x = 2, y = height - 5, labels = paste("第",pic_input-1000,"帧"), font = 2, cex = 2, adj = 0, col = "white")
    text(x = 2, y = height - 15, labels = paste("偏离角度",signif(angle, digits = 3)), font = 2, cex = 2, adj = 0, col = "white")
  }
  
  final <- c(final, angle)
}

#####
{
  dev.off()
  
  xstart=0; xend=100;
  ystart=0; yend=100;
  
  plot(NA,NA,xaxt='n', yaxt='n', xlim=c(xstart,xend),ylim=c(ystart,yend),xlab=NA, ylab=NA)
  text(x = 50, y = 60, labels = paste("输出波形图"), font = 2, cex = 4, col = "lightpink4")
  Sys.sleep(1)
}

All <- data.frame()
for (i in 1:23){
  file = paste("final/final",i,".out", sep = "")
  tmp = read.csv(file, header = T);
  All <- rbind(All, tmp$x)
}
value = colMeans(All)

#####输出波形图
{
  dev.off()
  xstart=0; xend=90;
  ystart=0; yend=55;
  jpeg("波形图.jpg", width=(xend-xstart), height=(yend-ystart), unit="cm", res=100)
  plot(NA,NA,xaxt='n', yaxt='n', xlim=c(xstart,xend),ylim=c(ystart,yend),xlab=NA, ylab=NA)
  
  input = read.csv("final.out", header = T)
  ymax = 15; ymin = -1; bwidth = 80; bheight = 40;
  
  segments(x0 = 10, x1 = 10 + bwidth, y0 = 10 + bheight*(5-ymin)/(ymax-ymin), lwd = 5, col = "grey")
  segments(x0 = 10, x1 = 10 + bwidth, y0 = 10 + bheight*(10-ymin)/(ymax-ymin), lwd = 5, col = "grey")
  segments(x0 = 10, x1 = 10 + bwidth, y0 = 10 + bheight*(15-ymin)/(ymax-ymin), lwd = 5, col = "grey")
  segments(x0 = 10, x1 = 10 + bwidth, y0 = 10 + bheight*(20-ymin)/(ymax-ymin), lwd = 5, col = "grey")
  
  
  xpos = 10 + bwidth*(1:30)/31; 
  ypos = 10 + bheight*(value - ymin)/(ymax - ymin)
  
  lines(xpos, ypos, lwd = 3, col = "cornflowerblue")
  
  xpos = 10 + bwidth*(1:30)/31; 
  ypos = 10 + bheight*(final - ymin)/(ymax - ymin)
  
  lines(xpos, ypos, lwd = 5, col = "red")
  
  count = 0;
  for (i in 1:30){
    if ((final[i] > value[i]) & (final[i] > 5)){
      count = count + 1
    }
  }
  rate = signif(1.4*100*count/30, digits = 2)
  
  segments(x0 = 10 + bwidth*c(1,5,10,15,20,25,30)/31, y0 = 10, y1 = 9, lwd = 5)
  text(x = 10 + bwidth*c(1,5,10,15,20,25,30)/31, y = 7, labels = c(1,5,10,15,20,25,30), font = 2, cex = 2)
  # text(x = 10 + bwidth/2, y = 5, labels = "帧数", font = 2, cex = 2)
  rect(xleft = 10, xright = 10 + bwidth, ytop = 10 + bheight, ybottom = 10, lwd = 5)
  
  
  segments(x0 = 10, x1 = 9, y0 = 10 + bheight*(c(0,5,10,15,20,25,30) - ymin)/(ymax - ymin), lwd = 5)
  text(x = 8, y = 10 + bheight*(c(0,5,10,15,20,25,30) - ymin)/(ymax - ymin), labels = c(0,5,10,15,20,25,30), font = 2, cex = 2, adj = 1)
  text(x = 2, y = 10 + bheight/2, labels = "偏离角度", font = 2, cex = 2, srt = 90)
  
  text(x = 10 + bwidth/2, y = 10 + bheight + 3, labels = paste("侧凸概率     ",rate,"%",sep = ""), font = 2, cex = 3)
  
  segments(x0 = 10, x1 = 20, y0 = 2, col = "red", lwd = 5)
  text(x = 21, y = 2, labels = "受检者侧凸曲线", font = 2, cex = 2, adj = 0)
  segments(x0 = 50, x1 = 60, y0 = 2, col = "cornflowerblue", lwd = 5)
  text(x = 61, y = 2, labels = "标准曲线", font = 2, cex = 2, adj = 0)
  dev.off()
}