# Reaping wall width from LOG files from direct output of Metamorph
#   without an intermediate 
# Enter variable parameters here:

Perc<- 0.85
dx<-25

# End of parameter entry


# Function Deffinitions -----------------------------

ask<- function(text) {
cat(text,": ", sep="")
indata<-readLines(con=stdin(),n=1,ok=TRUE)
                     }
                     
AddPeak<- function(iA, A, n) {
if (n==1) cat("Subtracting peak from", iA, sep=" ","\n")
if (n==0) cat("Removing force on", iA, sep=" ","\n")
A[iA,7]<- n
A
                             }
                     
AddPeak2<- function(X,iA, A, n, da) {
if (n==2) {cat("Adding peak to frame", iA, sep=" ","\n")
           layout(matrix(c(1,2,3), 1, 3, byrow = TRUE))
           N<-length(A[,1])
           Nl<- round(length(X[,1])/N)

           Y<-as.matrix(X[(4:Nl)+(iA-1)*Nl,])
           Xa<-spline(as.real(Y[,1]), as.real(Y[,3]))
           max0<- imax<-which.max(Xa$y);                 
           minx<- max(c(1, imax-da));      maxx<- min(c(length(Xa$y), imax+da))
           Xax<- Xa$x[minx:maxx];          Xay<- Xa$y[minx:maxx]
           imax<-which.max(Xay);           imin<-which.min(Xay)
           dimax<- 4
           jmax<- imax+dimax
           while (!dimax==0) {
                             jmax<-jmax+dimax   
                             plot(Xax, Xay, ty='l', lty= 1, main=iA)
                             lines(c(Xax[imax], Xax[imax]), c(Xay[imin], Xay[imax]), col='black', lty=3)
                             lines(c(Xax[jmax], Xax[jmax]), c(Xay[imin], Xay[imax]), col='red', lty=2)
                             dimax<- as.real(ask("Enter a change in 2nd line position"))
                             }
            Difx<- abs(sig<-Xax[jmax]-Xax[imax])
            Dify<- 0;  Diff<- (Difx^2 + Dify^2)^0.5 #   if (Diff==0) Diff<-0
            # n<- abs(n*sig/Difx)/(sig/abs(sig))
            A[iA,2:11]<- c(Xax[imax], Xax[jmax], Difx, Dify, Diff, 3, 0, minx-1+imax, 1, minx-1+jmax )
            
           }
A
                     }
                     
Add_A_Peak<- function(X,iA, A, da) {
imax<- as.real(ask(paste("Enter approximate max for frame", iA)))
           layout(matrix(c(1,2,3), 1, 3, byrow = TRUE))
           N<-length(A[,1])
           Nl<- round(length(X[,1])/N)

           Y<-as.matrix(X[(4:Nl)+(iA-1)*Nl,])
           Xa<-spline(as.real(Y[,1]), as.real(Y[,3]))
           minx<- max(c(1, imax-da));      maxx<- min(c(length(Xa$y), imax+da))
           Xax<- Xa$x[minx:maxx];          Xay<- Xa$y[minx:maxx]
           imax<-which.max(Xay);           imin<-which.min(Xay)

# Fix proximal peak
           dimax<- 1
           jmax<-imax+4
           while (!dimax==0) {
                             imax<-imax+dimax   
                             plot(Xax, Xay, ty='l', lty= 1, main=iA)
                             lines(c(Xax[imax], Xax[imax]), c(Xay[imin], Xay[imax]), col='black', lty=3)
                             lines(c(Xax[jmax], Xax[jmax]), c(Xay[imin], Xay[imax]), col='red', lty=2)
                             dimax<- as.real(ask("Enter a change in proximal line position"))
                             }
# Fix distal peak
           dimax<- 4
           jmax<- imax+dimax
           while (!dimax==0) {
                             jmax<-jmax+dimax   
                             plot(Xax, Xay, ty='l', lty= 1, main=iA)
                             lines(c(Xax[imax], Xax[imax]), c(Xay[imin], Xay[imax]), col='black', lty=3)
                             lines(c(Xax[jmax], Xax[jmax]), c(Xay[imin], Xay[imax]), col='red', lty=2)
                             dimax<- as.real(ask("Enter a change in distal line position"))
                             }
            Difx<- abs(sig<-Xax[jmax]-Xax[imax])
            Dify<- 0;  Diff<- (Difx^2 + Dify^2)^0.5 ;   # if (Diff==0) Diff<-0
            # n<- abs(n*sig/Difx)/(sig/abs(sig))
            A[iA,2:11]<- c(Xax[imax], Xax[jmax], Difx, Dify, Diff, 3, 0, minx-1+imax, 1, minx-1+jmax )
layout(matrix(c(1:10), 1, 10, byrow = TRUE))             
A
                     }
                     
Sort_Out<- function(A) {
  k<-length(A[,1])
  for (i in 1:k) { As<- sort(as.matrix(A[i,c(9,11)]))
                   A[i,c(9,11)]<-As[1:2]
                 }
  A                   }

# End of Function Deffinitions #----------------------------------------        
filnam<-"junk"
notuseit<- TRUE
if(fexist<- file.exists("LastFile")) {load("LastFile"); cat("Last file name used =", filnam,".")
	                         notuseit<-('n'== ask("Use it?(y/n)")) }
	                         if(notuseit) {
                                          print(list.files("./","log"))
                                          filnam<- print(ask("Enter log file name"))
	                                      }
	                                      # filnam<-"0614006qall.log"
XX<-readLines(filnam)
save(filnam, file="LastFile")
Nframes<-length(grep("Image Plane", XX))   # Automatic frame counting
XX<-read.csv(filnam, header=FALSE)
cat("This file has", Nframes, "Image Frames.\n")
iOff<- 0
ibrowse<-0
# Three types of startups:
# (1) Startup as a new data set 
# (2) Startup with a previous Out file saved as a textfile -or- 
# (3) Restart after a crash using current Out matrix

if ('y'==ask("Is this a new data set?(y/n)")) {
Nlines<- round(length(XX[,1])/Nframes)
Out<-matrix(0, Nframes,12)
Out[,1]<-1:Nframes   # index of frames
Out[,10]<-1           #  initialise force peak to: '0'= not forced ('1'= force 1 peak, '2'= force 2nd peak)
Out[,12]<-1           #  initialise jpat offset to: '0'= no offset
irep<-1
                                        } else {
cat("Available files: \n")
print(list.files("./",filnam))
iin<- as.real(ask("Enter the input index to load"))
if('n'== ask("Use the Out matrix in memory?")) { Out<-read.table(file=paste(filnam,iin,".txt", sep=''))} 
# Always give opportunity to set a new output suffix number 
}
                                            
irep<- as.real(ask("Enter the output index to append"))
# irep is the numeric suffix used to label the Output of the script.
# irepeat choses whether to analyze the data in another cycle
# Nframes are the number of image frames analyzed
# secondset

irepeat<-1
# Analysis Loops 
 while( irepeat > 0) {
par(mar=c(2.5, 2, 2, 0.1) + 0.1)
layout(matrix(c(1:10), 1, 10, byrow = TRUE)) 

i<-0
while( i < Nframes+10) {
i<-i +1
if( i<Nframes+1) {
secondset<- FALSE; nosecond<- TRUE
YY<-as.matrix(XX[(4:Nlines)+(i-1)*Nlines,])
Xo<-spline(as.real(YY[,1]), as.real(YY[,3]))
Yo<-spline(as.real(YY[,1]), as.real(YY[,2]))
Yos<-Yo$y
if(abs(Out[i,7])< 2) {imax<-which.max(Xo$y); Out[i, c(9,11)]<- imax}                
if(abs(Out[i,7])==3) {secondset<-TRUE; nosecond<- FALSE; imax<- Out[i, 9]; jmax<- Out[i, 11]; dij<- jmax-imax}
if(abs(Out[i,7])==2) {secondset<-TRUE; nosecond<- FALSE; imax<- Out[i, 9]; jmax<- Out[i, 11]; dij<- jmax-imax} 
minx<- max(c(1, sx<-imax-dx));      maxx<- min(c(length(Xo$y), imax+dx))
Xox<- Xo$x[minx:maxx];          Xoy<- Xo$y[minx:maxx]
Yoy<- Yos[minx:maxx];           XoyCrit<-Perc*max(Xoy)

if(abs(Out[i,7]) < 2) {imax<-which.max(Xoy); } else {if(sx< 1) {imax<- dx+1+sx; jmax<- imax + dij} else { imax<- dx+1; jmax<- imax + dij }}

imin<-which.min(Xoy)  
plot(Xox, Xoy, ty='l', lty= 1, main=i)
lines(c(Xox[imax], Xox[imax]), c(Xoy[imin], Xoy[imax]), col='black', lty=3)
lines(Xox[c(1,maxx-minx)],c(XoyCrit, XoyCrit), col="orange", lty=2, lwd=2)
maxflat<-FALSE;   maxi<-length(Xox); lowi<-imax; hii<-imax  
if (Out[i, 7]==1) {nosecond<-TRUE} else  {
       nosecond<- FALSE;
       if((abs(Out[i, 7])) > 1) {maxflat<-TRUE      # changed from ==2                          
                        secondset<- TRUE
                        # sig<- Out[i, 7]/2
                        # dif1<-Xox[imax+1]-Xox[imax]
                        npeak<-2
                               }  else {
       while (maxflat==FALSE) {
	           lowi<-lowi-1
	           hii<-hii+1 
	           if (lowi==0) {lowi<-lowi+1}
	           if (hii > maxi) {hii<-hii-1}
	           if (lowi==1) { if (hii==maxi) {maxflat<-TRUE; nosecond<-TRUE}}
	           Xoyh<- Xoy[hii-1]
	           Xoyl<-Xoy[lowi+1]
 	           if(Xoy[hii]<= XoyCrit) {Xoy[hii]<-XoyCrit}
	           if(Xoy[lowi]<= XoyCrit) {Xoy[lowi]<-XoyCrit}
               if(Xoy[hii]<= Xoyh) {Xoy[imax:hii]<-Xoy[hii]} else {maxflat<-TRUE}
               if(Xoy[lowi]<= Xoyl) {Xoy[lowi:imax]<-Xoy[lowi]} else {maxflat<-TRUE}          
               if(maxflat==TRUE) {Xoy[lowi:hii]<- min(c(Xoy[hii],Xoy[lowi]))}
	                    } # end while maxflat==FALSE
	                    jmax<-which.max(Xoy); npeak<- 2; Out[i,7]<- 2; Out[i,9]<- minx-1+imax; Out[i,11]<- minx-1+jmax
	                                    } # end else
	                    } # end else nosecond == FALSE

if(nosecond) {if(secondset==FALSE) {jmax<-imax; npeak<- 1; Out[i,7]<- 1; }} #  else {jmax<-which.max(Xoy); npeak<- 2}}

if(!nosecond) {lines(c(Xox[jmax], Xox[jmax]), c(Xoy[imin], Xoy[jmax]), col='green', lty=4)}
mtext(npeak, line=-3, adj=1, cex=2, col = 'blue') 
if(npeak==1) Out[i,9] <- Out[i,11]
cat(i, "Max ix,y:",Out[i,9],Out[i,11],"| type:", Out[i,7], Xox[imax], Xox[jmax], "dif(X,Y)=", Difx<- abs(Xox[imax]-Xox[jmax]),
      Dify<- abs(Yoy[imax]-Yoy[jmax]), "diff=", Diff<- (Difx^2 + Dify^2)^0.5, "\n")
      Out[i,2:6]<- c(Xox[imax], Xox[jmax], Difx, Dify, Diff)
        } # end of "if i < Nframes" allowing truncating profiles at end of frames 
         
 if (i==ibrowse) {browser(); ibrowse<-0}
 if (((i+iOff)/10)== round((i+iOff)/10)) {mtext(filnam, line=-2, adj=0, cex=1, col = 'red')
            cat(" back | fwd | 98% dx | 102% dx | tog-5 | +2nd | -2nd | no_force | goto | browse | cntr | sort | wid | nar | abort |\n")
      char<-ask("  <   |  >  |   -    |    +    |   5   |   2  |   1  |    0     |   g  |   b    |  c   |  s   |  w  |  n  |   *   |")
                            if (i>10) { if (char=='<'){i<-i-20}                                      }
                            if (char=='b'){ibrowse<- as.real(ask("Enter frame number to browse at"));i<-i-10}
                            if (char=='c'){Out<- Add_A_Peak(XX, as.real(ask("Enter frame number at which to insert a peak")), Out, dx);i<-i-10}
                            if (char=='g'){i<- as.real(ask("Goto frame number "));i<-i-5; cat("Go to frame",i,"as center frame. \n")}
                            if (char=='-'){ Perc<-0.98*Perc; i<-i-10; cat("Criterion increased to",Perc,"of peak. \n")}
                            if (char=='5'){ if(iOff==0) {iOff<-5} else {iOff<-0}; i<-i-5}
                            if (char=='+'){ Perc<-Perc/0.98; i<-i-10; cat("Criterion decreased to",Perc,"of peak. \n")}
                            if (char=='w'){ dx<-dx+1; i<-i-10; cat("dx widened to",dx,"\n")}
                            if (char=='n'){ dx<-dx-1; i<-i-10; cat("dx narrowed to",dx,"\n")}
                            if (char=='*'){ stop("Abort this program")}
                            if (char=='2'){Out<- AddPeak2(XX, as.real(ask("Enter frame number to add peak")),Out,2, dx);i<-i-10
                                           layout(matrix(c(1:10), 1, 10, byrow = TRUE))
                                          }
                            if (char=='1'){Out<- AddPeak(as.real(ask("Enter frame number to subtract peak")),Out,1);i<-i-10}
                            if (char=='0'){Out<- AddPeak(as.real(ask("Enter frame number to remove force")),Out,0);i<-i-10}
                            if (char=='s'){Out<- Sort_Out(Out);i<-i-10}
                            
                            }  # end i/10
                                     }   # end of 'while i < NFrames + 1'
mtext(filnam, line=-2, adj=0, cex=1, col = 'red')
write.table(Out, file=paste(filnam,irep,".txt", sep=''))
browser()
par(mar=c(4, 4, 4, 4) + 0.1)
layout(matrix(c(1:10), 1, 1, byrow = TRUE))

plot(Out[,1], Out[,6], main = paste(filnam,"DIFF", "rep =",irep), ty='l')
points(Out[,1], Out[,6], col='red')
        irepeat<- as.real(ask("Enter 0 to end | 1 to add another analysis"))
        if (irepeat>0) irep<-irep+1
         }  # end of while irepeat > 0