Tracking: vídeo-análise pelo Mathematica

Prof. Douglas M G Leite - ITA

1- Define opções e diretório de trabalho (mesmo deste notebook), mostra lista de arquivos de vídeo no diretório

Off[General::"spell1"]; Off[General::"spell"]; Off[General::sysffmpeg];Off[Remove::rmnsm];
SetSystemOptions["VideoProcessingOptions"->{"TryUsingSystemFFmpeg"->False}];
scrsize=AbsoluteOptions[EvaluationNotebook[],WindowSize][[1,2,1]];
tam=1.5;
tam=1.5;
dtam=0.001;
frame=1;
grand={"t (s)","x (m)","y (m)","r (m)","θ (rad)"};
PointInPoly[{x_,y_},poly_List]:=Module[{i,j,c=False,npol=Length[poly]},For[i=1;j=npol,i<=npol,j=i++,If[((((poly[[i,2]]<=y)&&(y<poly[[j,2]]))||((poly[[j,2]]<=y)&&(y<poly[[i,2]])))&&(x<(poly[[j,1]]-poly[[i,1]])*(y-poly[[i,2]])/(poly[[j,2]]-poly[[i,2]])+poly[[i,1]])),c=¬c];];c]
stilo={FontFamily->"Arial",16,Plain};
opt=SetOptions[{ListPlot,Plot},PlotStyle->{Green,PointSize[0.01]},PlotRange->All,ImageSize->scrsize/4,Frame->True,GridLines->Automatic,BaseStyle->stilo,FrameStyle->Directive[stilo],LabelStyle->stilo];
ParallelEvaluate[SetOptions[ListPlot,opt[[1]]]];
ParallelEvaluate[SetOptions[Plot,opt[[2]]]];
ParallelEvaluate[opt];
diretorio=DirectoryName[ToFileName["FileName"/.NotebookInformation[EvaluationNotebook[]]]];
$Path=Append[$Path,diretorio];
Print["Diretório de Trabalho = ",SetDirectory[diretorio]]
file=Flatten[Join[FileNames["*.mp4"],FileNames["*.avi"]]];
Print[Column[{"Arquivos de vídeo na Pasta de Trabalho:",Transpose[{Range[Length[file]],file}]//TableForm}]]

2 - Escolha do vídeo original

narquivo = 7; (* informar o número correspondente ao nome do arquivo do vídeo a ser analisado conforme ista acima *)
nomevideo=file[[narquivo]];
video0=Video[nomevideo];
info=Flatten[Information[video0,{"Duration","FrameCount","FrameRate","OriginalRasterSize"}]];
{dur,num,frate,size0}={QuantityMagnitude[UnitConvert[info[[1]]]],(info[[2]]//Values)[[1]],QuantityMagnitude[(info[[3]]//Values)[[1]]]*1.0,Flatten[info[[4]]//Values]};
Print[Column[{"Dados do Arquivo:",
Row[{"nome          = ",nomevideo}],
Row[{"resolução     = ",size0," pixels"}],
Row[{"duração       = ",dur," segundos"}],
Row[{"taxa          = ",frate," quadros/seg"}],
Row[{"total quadros = ",num}],
Framed[Thumbnail[VideoExtractFrames[video0,1]]]}]]
Print[">>>>>>> memória em uso / memória disponível = ",IntegerPart[MemoryInUse[]/1000000],"MB / ",IntegerPart[MemoryAvailable[]/1000000], "MB"]

3 - Carrega vídeo em forma de uma sequência de imagens (quadros)

size1=1280;  (* resolução máxima dos quadros - pixels na horizontal *)
usom=1/2; 
mem=MemoryAvailable[];
mrot=IntegerPart[(mem*Min[1,usom])/(size0[[1]]*size0[[2]]*4)];
rot=IntegerPart[num/IntegerPart[(num/mrot)+1]+1];
Clear[fun];
If[size1<size0[[1]],fun[x_]:=ExportString[ImageResize[x,size1],"JPEG"],fun[x_]:=ExportString[x,"JPEG"]];
PrintTemporary["Aguarde..."];
Clear[imgs];
imgs={};
Monitor[tempo= Round[AbsoluteTiming[Do[
imgimp=VideoExtractFrames[video0,Quantity[Interval[{i,Min[i+rot-1,Round[num]]}],"Frames"]];
imgst=Flatten[ParallelMap[fun,imgimp]];
Clear[imgimp];
imgs=Flatten[Join[imgs,imgst]];
Clear[imgst],{i,1,num,rot}];][[1]]],Dynamic[Column[{ProgressIndicator[i,{0,num}],Row[{">>>>>>> memória em uso / memória disponível = ",IntegerPart[MemoryInUse[]/1000000],"MB / ",IntegerPart[MemoryAvailable[]/1000000], "MB"}]}]]];
Print[Length[imgs]," quadros importados em ", tempo," segundos"]
size=ImageDimensions[img0=ImportString[imgs[[1]],"JPEG"]];
imgq=imgi=img0;
framef=num=Length[imgs];
imgf=ImportString[imgs[[num]],"JPEG"];
Print[">>>>>>> memória em uso / memória disponível = ",IntegerPart[MemoryInUse[]/1000000],"MB / ",IntegerPart[MemoryAvailable[]/1000000], "MB"]

4 - Painel de controle para seleção de quadros inicial e final, e definir sistema coordenadas, escala e imagem de referência

in=10;
sizescr=scrsize/2.5;
rt0=Round[sizescr/50.];
Manipulate[If[ValueQ[size1],
If[in≠i,in=i;imgq=ImportString[imgs[[in]],"JPEG"]];
If[vbot1,frame=in;imgi=imgq;vbot1=False];
If[vbot2,framef=in;imgf=imgq;vbot2=False];
th0=vth0*π/180;
pf=vpf;
reg=EuclideanDistance[vrg1,vrg2];
p0=vp0;
rt=vrt;
limitr=Round[{p0-vrt,p0+vrt}];
limits0={vl1,vl2};
imgu={imgq,imgi,imgf};
ref=ImageTrim[imgi,limitr];];
Show[{imgu[[vim]],Graphics[{Blue, Dashed,Text[Style["x",18],vpf+{Cos[th0],Sin[th0]}*size1/10-{12,12}],Arrowheads[0.025],Arrow[{vpf,vpf+{Cos[th0],Sin[th0]}*size1/10}],Text[Style["y",18],vpf+{-Sin[th0],Cos[th0]}*size1/10-{12,12}],Arrow[{vpf,vpf+{-Sin[th0],Cos[th0]}*size1/10}],Line[{vpf,vpf+{Cos[th0],Sin[th0]}*size1}],Line[{vpf,vpf-{Cos[th0],Sin[th0]}*size1}],Line[{vpf,vpf+{Sin[th0],-Cos[th0]}*size1}],Line[{vpf,vpf-{Sin[th0],-Cos[th0]}*size1}]}],
Graphics[{Purple, Dashed,Line[{vrg1,vrg2}],Text[Style[Row[{Around[tam,dtam]," m"}],18,Gray,Background->Black],(vrg1+vrg2)/2-{0,15}]}],Graphics[{Red, Dashed,Line[{vp0,pf}]}],Graphics[{FaceForm[],EdgeForm[{Dashed,Thick,Red}],Rectangle[vp0-vrt,vp0+vrt]}],Graphics[{FaceForm[],EdgeForm[{Dashed,Thick,Darker[Green]}],Rectangle[vl1,vl2]}]},ImageSize->{UpTo[Round[sizescr*1.5]],UpTo[Round[sizescr*2]]}],
{{vbot1,False},None},{{vbot2,False},None},{{vbot3,False},None},
{{vrg1,{size[[1]]/4,size[[2]]/4}},{0,0},{size[[1]],size[[2]]},Locator,Appearance->Style["+",Purple,FontSize->48]},
{{vrg2,{size[[1]]/2.5,size[[2]]/4}},{0,0},{size[[1]],size[[2]]},Locator,Appearance->Style["+",Purple,FontSize->48]},
{{vpf,{size[[1]]/2,size[[2]]/2}},{0,0},{size[[1]],size[[2]]},Locator,Appearance->Style["+",Blue,FontSize->48]},
{{vp0,{size[[1]]/1.5,size[[2]]/1.5}},{0,0},{size[[1]],size[[2]]},Locator},
{{vl1,{rt0,rt0}},{0,0},{size[[1]],size[[2]]},Locator,Appearance->Graphics[{Darker[Green],Disk[{0,0},2]},ImageSize->rt0]},
{{vl2,{size[[1]]-rt0,size[[2]]-rt0}},{0,0},{size[[1]],size[[2]]},Locator,Appearance->Graphics[{Darker[Green],Disk[{0,0},2]},ImageSize->rt0]},
Column[{Control[{{i,1,"Cursor quadro"},1,num,vi,Appearance->"Open",AnimationRate->frate}],
Control[{{vi,10,"Passo entre quadros"},{1,2,5,10,20,50,100,200},SetterBar}],
Button["Fixar como quadro inicial",vbot1=True],
Button["Fixar como quadro final",vbot2=True]," ",
Row[{"Mostra quadro: ",Control[{{vim,1," "},{1->"atual" ,2->"inicial",3->"final"},SetterBar }]}]," ",
Control[{{vth0,0,"Ângulo eixos"},0,360,0.1,Appearance->"Open"}]," ",
Row[{"Régua = ",InputField[Dynamic[tam],Number,FieldSize->{6,1}]," ± ",InputField[Dynamic[dtam],Number,FieldSize->{6,1}]," m"}]," ",
"Imagem a ser rastreada",Framed[Dynamic[Show[{ref,Graphics[{Red,Circle[{rt+1,rt+1},(rt+1)/1.2],Line[{{rt+1,0},{rt+1,2(rt+1)}}],Line[{{0,rt+1},{2(rt+1),rt+1}}]}]},ImageSize->10*rt0]],FrameStyle->{Dashed,Thick,Red},FrameMargins->0],
Control[{{vrt,rt0,"Área objeto"},4,Round[size0[[1]]/5],1,Appearance->"Open"}]},Alignment->Right,Spacings->0.5],ControlPlacement->Left]//Quiet

5 - Realizando Tracking

passo=1; (* passo entre os quadros para rastreio *)
mm=pontos={};
ident=data=rpontos={};
gxp={};
stop=False;
cont=frame;
lastpont=pont={1,p0};
senmax=sen=0.99;
senmin=0.5;
trm=1.5;
ir=1;
trest=100;
vel=tempof=1;
skip = Max[1,IntegerPart[passo]];
limits1={{0,0},limits0[[2]]};
sizescr2=scrsize/3;
Quiet[
tempo=Round[AbsoluteTiming[
imgt=ImageTrim[ImportString[imgs[[cont]],"JPEG"],limits0];
imgx = ColorNegate[ImageCorrelate[imgt, ref, NormalizedSquaredEuclideanDistance]];
mm=ComponentMeasurements[Binarize[imgx,sen],{"Centroid","Count"}] ;
While[Length[mm]<2|| mm[[1,2,2]]<2*rt ,
sen=sen/1.01;
mm=ComponentMeasurements[Binarize[imgx,sen],{"Centroid","Count"}] ];
senmin=sen*1.02;
sen=senmax;
mm={};
Monitor[While[stop==False &&cont≤framef,
tempof=Round[(AbsoluteTiming[
imgw=ImageTrim[ImportString[imgs[[cont]],"JPEG"],limits0];
If[mm=={}||FractionalPart[cont/(Round[frate]*skip)]==0,limits1={{0,0},limits0[[2]]};imgt=imgp=imgw;sen=senmax,
vel=trm*EuclideanDistance[pont[[2,1]],lastpont[[2,1]]];
limits1={pont[[2,1]]-vel-2rt-limits0[[1]],pont[[2,1]]+vel+2rt-limits0[[1]]};
imgt=ImageTrim[imgp=imgw,limits1]];
imgx = ColorNegate[ImageCorrelate[imgt, ref, NormalizedSquaredEuclideanDistance]];
mm=ComponentMeasurements[DeleteBorderComponents[Dilation[Binarize[imgx,sen],DiskMatrix[rt/2]]],{"Centroid","Count"}] ;
While[mm=={}&&sen >senmin*1.02,
sen=sen/1.02;
mm=ComponentMeasurements[DeleteBorderComponents[Dilation[Binarize[imgx,sen],DiskMatrix[rt/2]]],{"Centroid","Count"}] ];
While[(Length[mm]>1|| mm[[1,2,2]]>6(rt/2)^2)&&sen<senmax,
sen=sen*1.02;
If[sen≥1,sen=0.99];
mm=ComponentMeasurements[DeleteBorderComponents[Dilation[Binarize[imgx,sen],DiskMatrix[rt/2]]],{"Centroid","Count"}] ];
If[mm=={},anim=Show[imgp,ImageSize->sizescr2];If[limits1=={{0,0},limits0[[2]]},iden={};ident=AppendTo[ident,iden],cont=cont-skip],
lastpont=pont;
pont={cont-frame,Sort[mm,#1[[All,2,2]]>#2[[All,2,2]]&][[1,2]]+{limits0[[1]]+{Max[limits1[[1,1]]-1,0],Max[limits1[[1,2]]-1,0]},0}};
rpont={act=(cont-frame)*1./frate,RotationTransform[-th0][pont[[2,1]]-pf]*tam/(reg)};
limitr1={pont[[2,1]]-rt-limits0[[1]],pont[[2,1]]+rt-limits0[[1]]};
dat={rpont[[1]],rpont[[2,1]],rpont[[2,2]],Norm[rpont[[2]]],ArcTan[rpont[[2,1]],rpont[[2,2]]],cont,quali=Round[100*(sen-senmin)/(senmax-senmin),10]};
pontos=AppendTo[pontos,pont];
rpontos=AppendTo[rpontos,rpont];
data=AppendTo[data,dat];
anim=Show[imgp,iden={Graphics[{Hue[quali/300],Thick,Circle[pont[[2,1]]-limits0[[1]],(rt+1)/1.2]}],Graphics[{FaceForm[],EdgeForm[{Dashed,Thick,Darker[Green]}],Rectangle[limits1[[1]],limits1[[2]]]}]},ImageSize->sizescr2];
ident=AppendTo[ident,iden];];][[1]]*1000+tempof*10)/11];
If[cont==frame||stop||FractionalPart[cont/(10*skip)]==0,
trest=Round[((framef-cont)*tempof/(skip*1000))];
rdat=data;
cdata=SortBy[GatherBy[data,Last],#[[1,7]]&,Greater];
escala=Map[Framed[Style[NumberForm[#,{3,0},NumberPadding->{" ","0"}],Black],Background->Hue[#/300]]&,cdata[[All,1,7]]];
grafs=Table[ParallelMap[ListPlot[#[[All,{1,i}]],PlotRange->All,PlotStyle->{PointSize[0.01],Hue[#[[1,7]]/300]},FrameLabel->{grand[[1]],grand[[i]]}]&,cdata],{i,2,5}]];
cont=cont+skip],
If[cont==frame,"Aguarde, ajustando sensibilidade...",
Framed[Row[{Column[{"Rastreando, aguarde...",
Row[{"Quadro # ",cont}],
Row[{"t= ",act," s"}],
Row[{"Confiança = ",quali," % "}],
Row[{"Tempo por quadro = ",tempof," ms"}],
Row[{"Tempo restante = ",trest," s"}],
ProgressIndicator[cont,{frame,framef}],
Button["Parar",stop=True],anim,
Pane[Row[{{Map[Show[#,PlotRange->All]&,grafs[[1;;2]]],Map[Show[#,PlotRange->All]&,grafs[[3;;4]]]}//TableForm,"   ",Column[{"Confiança (%)",Column[escala]},Alignment->Center]}],{scrsize/1.5,scrsize/2.5}]}]}]]]];][[1]]]];
Print[Length[pontos]," pontos coletados de um total de ",IntegerPart[(cont-frame)/skip]," quadros analisados em ",tempo," segundos"];
rdat=data;
cdata=SortBy[GatherBy[data,Last],#[[1,7]]&,Greater];
escala=Map[Framed[Style[NumberForm[#,{3,0},NumberPadding->{" ","0"}],Black],Background->Hue[#/300]]&,cdata[[All,1,7]]];
grafs=Table[ParallelMap[ListPlot[#[[All,{1,i}]],PlotRange->All,PlotStyle->{PointSize[0.01],Hue[#[[1,7]]/300]},FrameLabel->{grand[[1]],grand[[i]]}]&,cdata],{i,2,5}];
Row[{{Map[Show[#,PlotRange->All]&,grafs[[1;;2]]],Map[Show[#,PlotRange->All]&,grafs[[3;;4]]]}//TableForm,"   ",Column[{"Confiança (%)",Column[escala]},Alignment->Center]}]

6 - Removendo pontos ruins pela qualidade (semelhança com a referência)

Manipulate[limq=lq;If[ValueQ[rdat],
Row[{{Map[Show[#,PlotRange->All]&,grafs[[1;;2,1;;lq]]],Map[Show[#,PlotRange->All]&,grafs[[3;;4,1;;lq]]]}//TableForm,"   ",Column[{"Confiança (%)",Column[escala]},Alignment->Center]}]],
{{lq,Length[escala],"Limite Confiança"},Table[i->escala[[i]],{i,Length[escala],1,-1}],Setter},Button["Finalizar Seleção",rdat=SortBy[Flatten[cdata[[1;;limq]],1],First]]]

7 - Removendo pontos ruins direto do gráfico

vbot1=vbot2=False;
ndat={};
are=0;
Manipulate[If[ValueQ[are],
If[vbot1,rdat=Delete[rdat,Flatten[Map[Position[art,#]&,ndat],1]];vbot1=False];
If[vbot2,rdat=Map[#-{rdat[[1,1]],0,0,0,0,0,0}&,rdat];are=0;vbot2=False];
art=rdat[[All,{1,var}]];
If[are≠var,are=var;ndat={};range={MinMax[art[[All,1]]],MinMax[art[[All,2]]]};midp=Midpoint[Transpose[range]]];
Manipulate[
ndat=Select[art,PointInPoly[#,{p1,p2,p3,p4}]==True&];
Show[ListPlot[art,PlotRange->range,PlotRangePadding->Scaled[0.1],FrameLabel->{grand[[1]],grand[[var]]},PlotStyle->Green],ListPlot[ndat,PlotStyle->{PointSize[0.01],Darker[Red]}],Graphics[{EdgeForm[{Dashed,Red}],Blue,Opacity[0.1],Polygon[{p1,p2,p3,p4}]}],ImageSize->scrsize/3],
{{p1,midp},Locator},{{p2,midp},Locator},{{p3,midp},Locator},{{p4,midp},Locator},
Button["Remover",vbot1=True],Button["Zerar Abcissa",vbot2=True]]],
{{var,2,"Escolha a grandeza"},{2->" x ",3->" y ",4->" r ",5->" θ "}},AutoAction->False]

8 - Removendo pontos ruins pelo vídeo

torem={};
Manipulate[If[ValueQ[torem],
imgr=ImageTrim[ImportString[imgs[[rdat[[fi,6]]]],"JPEG"],limits0];
Column[{Row[{"pontos a serem removidos = ",torem}],Show[imgr,If[MemberQ[torem,fi],Graphics[{Red,Text[Style["X",rt*2,Red],pontos[[rdat[[fi,6]]-frame+1,2,1]]-limits0[[1]]]}],ident[[rdat[[fi,6]]-frame+1]]],Graphics[{Text[Style[fi,Round[rt/1.5],Green,Background->Black],pontos[[rdat[[fi,6]]-frame+1,2,1]]-{0,rt*1.2}-limits0[[1]]]}],ImageSize->scrsize/3]}]],
{{fi,1,"ponto"},1,Length[rdat],1,Appearance->"Open",AnimationRate->frate*skip},
Button["Remover",torem=AppendTo[torem,fi]],
Button["Limpar Lista",torem={}],
Button["Finalizar",rdat=Delete[rdat,Map[{#}&,torem]];torem={}]]

9 - Exportando em arquivo texto

inc=Norm[{tam/reg,dtam,(rt/8)*0.5*tam/reg}];
medr=Mean[rdat[[All,4]]]-Norm[{StandardDeviation[rdat[[All,4]]],inc}];
incertezas={0.5/frate,inc,inc,inc,ArcTan[medr,inc],0,0};
decdig=-Map[MantissaExponent[#][[2]]-If[MantissaExponent[#][[1]]>0.35,1,2]&,incertezas];
precdig=decdig+Max[3,MantissaExponent[size1*tam/reg][[2]]+1];
edat=Transpose[Table[Map[NumberForm[#,{precdig[[i]],decdig[[i]]},NumberPadding->{" ","0"}]&,Transpose[Join[{incertezas},rdat]][[i]]],{i,1,Length[rdat[[1]]]}]];
(nedat=Join[{{"t(s) ","x(m) ","y(m) ","r(m) ","θ(rad)","Frame ","Conf(%)"}},{edat[[1]]},{"-----"},edat[[2;;All]]])[[1;;Min[10,Length[edat]]]]//TableForm
Export["track_"<>FileBaseName[nomevideo]<>".txt",nedat,"Table"]