//$MENU window=fprintimage;popup=Normalization;insertafter=auto assign using standard (curves);name=Auto assign (sequencer profiles)

string st,st1,comments;
integer refsyscount,lanecount,laneresol,lanenr,i,j,p,cert,bandcount,ok,refband1,refband2,band1,band2;
float h,w,refps1,refps2,ps1,ps2,a,b;
float refsyspos[],refsysdiff[],bandh[],diff,score,bestscore;
integer bandpos[],bestband1,bestband2;
DIALOG dlg;


if AnFprGetCurrent="" then AnFprAttach;
if AnFprGetCurrent="" then { message("No fingerprint analysis present"); stop; }

st=AnFprGetRefSystem;
refsyscount=0;
while st<>"" do {
   st1=splitstring(st,"	");
   splitstring(st1,":");
   refsyscount=refsyscount+1;
   refsyspos[refsyscount]=val(st1);
}
if refsyscount<2 then { message("Insufficient number of reference bands"); stop; }

for i=1 to refsyscount do {
   if (i>1) and (i<refsyscount) then refsysdiff[i]=(refsyspos[i+1]-refsyspos[i-1])/2;
   if i=1 then refsysdiff[i]=refsyspos[i+1]-refsyspos[i];
   if i=refsyscount then refsysdiff[i]=refsyspos[i]-refsyspos[i-1];
   refsysdiff[i]=refsysdiff[i]/3;
}

lanecount=AnFprGetLaneCount;
for lanenr=1 to lanecount do {
   for i=1 to refsyscount do AnFprDelRefLink(lanenr,i);
   laneresol=AnFprGetCurveResol(lanenr);
   bandcount=0;
   for i=1 to AnFprGetBandCount(lanenr) do {
      AnFprGetBandInfo(lanenr,i,p,h,w,cert);
      if cert then { bandcount=bandcount+1; bandpos[bandcount]=p; bandh[bandcount]=h; }
   }
   if bandcount=0 then comments=comments+"Lane "+str(lanenr,0,0)+" has no bands assigned!~n";
   i=1;
   while i<=bandcount do {
      ok=1;
      if i>1 then if bandh[i]<bandh[i-1]*0.3 then ok=0;
      if i>2 then if bandh[i]<bandh[i-2]*0.2 then ok=0;
      if i<bandcount then if bandh[i]<bandh[i+1]*0.3 then ok=0;
      if i<bandcount-1 then if bandh[i]<bandh[i+2]*0.2 then ok=0;
      if i>1 then if bandpos[i]-bandpos[i-1]<laneresol/500.0 then if bandh[i]<bandh[i-1] then ok=0;
      if i<bandcount then if bandpos[i+1]-bandpos[i]<laneresol/500.0 then if bandh[i]<bandh[i+1] then ok=0;
      if not(ok) then {
         bandcount=bandcount-1;
         for j=i to bandcount do { 
            bandh[j]=bandh[j+1];
            bandpos[j]=bandpos[j+1];
         }
      }
      else i=i+1;
   }
   
   refband1=3;
   refband2=refsyscount-3;
   refps1=refsyspos[refband1];
   refps2=refsyspos[refband2];
   bestscore=refsyscount*0.5;bestband1=-1;
   for band1=1 to bandcount-1 do
      for band2=band1+1 to bandcount do {
         ps1=bandpos[band1]/laneresol*100;         
         ps2=bandpos[band2]/laneresol*100;
         b=(refps2-refps1)/(ps2-ps1);
         a=refps1-ps1*b;
         if (b>1.0/1.5) and (b<1.5) then {
            j=1;score=0;
            for i=1 to refsyscount do {
               while (j<bandcount) and 
                     (abs(refsyspos[i]-a-b*bandpos[j+1]/laneresol*100) <
                      abs(refsyspos[i]-a-b*bandpos[j]/laneresol*100)) do j=j+1;
               diff=abs(refsyspos[i]-a-b*bandpos[j]/laneresol*100);
               diff=1.0-diff/refsysdiff[i];
               if diff>0 then score=score+diff*diff;
            }
            if score>bestscore then {
               bestscore=score;
               bestband1=band1;
               bestband2=band2;
            }
         }
      }
      
   if bestband1>0 then {
      ps1=bandpos[bestband1]/laneresol*100;         
      ps2=bandpos[bestband2]/laneresol*100;
      b=(refps2-refps1)/(ps2-ps1);
      a=refps1-ps1*b;
      
      j=1;score=0;
      for i=1 to refsyscount do {
         while (j<bandcount) and 
               (abs(refsyspos[i]-a-b*bandpos[j+1]/laneresol*100) <
                abs(refsyspos[i]-a-b*bandpos[j]/laneresol*100)) and
               (abs(refsyspos[i]-a-b*bandpos[j+1]/laneresol*100) <
                abs(refsyspos[i+1]-a-b*bandpos[j+1]/laneresol*100)) do j=j+1;
         diff=abs(refsyspos[i]-a-b*bandpos[j]/laneresol*100);
         diff=1.0-diff/refsysdiff[i];
         if diff>0 then AnFprSetRefLink(lanenr,i,bandpos[j]);
      }
      
//      AnFprSetRefLink(lanenr,refband1,bandpos[bestband1]);
//      AnFprSetRefLink(lanenr,refband2,bandpos[bestband2]);
   }
    
}

if comments<>"" then {
   DlgAddEdit(dlg,comments,15,15,400,300);
   DlgAddButton(dlg,"Close",11,350,325,65,25);
   DlgShow(dlg,"Error report",450,400);
}