(*
<pre>
Image Sequence Discrimination Model of
Robert Eriksson and Al Ahumada *)
(* Misc. functions *)

Ds[x_] := Dimensions[x]          
Sqr[x_]:= x*x

Mean2D[image_] :=
        Mean[Map[Mean,image]]

Mean3D[seq_] :=
        Mean[Map[Mean2D,seq]]

grey2lum[sequence_] :=
        intrcpt + slope sequence

(*-----------------------------------------------*) 
	
(* Energy *)

sigE[sig_, displayParms_] :=
        Apply[Plus,
Sqr[Flatten[sig]]]/Apply[Times,displayParms[[{1,2,3}]]]//N


signalEnergy[signal_, dt_:1/76., resolution_:64] :=
        (* Computes energy in an arbitrary  contrast signal: *)

        Module[{sigsum,loge},
        Apply[Plus, Sqr[Flatten[signal]]]*dt/resolution^2
        ];

dBB[e_] :=
        (* dB Barlow *)
        20.(Log[10,e] + 6.0);
(* Computes energy in an arbitrary  contrast signal: *)
(*------------------------------------------------------------*)

(* Filtering *)

tFiltPar[st_,fs_]:= Exp[-Log[2.] / (fs*st)];

gaussFreqList[nseq_, s_] := RotateRight[gaussWindow[nseq, s],
                                        Floor[nseq/2]]

sFilterG[imageDim_, sx_, sy_:sx] := 
  Outer[Times, gaussFreqList[imageDim[[1]], sy], 
   gaussFreqList[imageDim[[2]], sx]]

(*-----------------------------------------------*)
(* Models *)

visLinModelE[sequence0_,sequence1_, visParms_, displayParms_] :=
	Module[{length,rows,cols, sampParms, te, fSeq0, fSeq1,
	coneI, ac,ah,ap, filtc,filth,filtm,filtq, sigm,sigp, gm,gp},
			
	(* convert to sampled frequencies: *)
	{length,rows,cols} = Dimensions[sequence0];
	sampParms = {	visParms[[1]]*displayParms[[1]], 
			visParms[[2]]*rows/displayParms[[2]],
			visParms[[2]]*cols/displayParms[[3]], 
			visParms[[3]]};
	
	(* process the background sequence *)
	filtc = sFilterG[{rows,cols}, sampParms[[3,1]], sampParms[[2,1]] ];
	fSeq0 = (filtc*InverseFourier[#])& /@ sequence0;
	fSeq1 = (filtc*InverseFourier[#])& /@ sequence1;

	coneI = filtc*InverseFourier[Table[grey2lum[displayParms[[4]]],
                {rows},{cols}] ];
	
	ac = tFiltPar[sampParms[[1,1]], displayParms[[1]] ];
	ah = tFiltPar[sampParms[[1,2]], displayParms[[1]] ];
	ap = tFiltPar[sampParms[[1,3]], displayParms[[1]] ];

	te = 1+Floor[Log[1./E//N]/Log[ap] ]; (* ensures drop to 1./E *)

	filth = sFilterG[{rows,cols}, sampParms[[3,2]], sampParms[[2,2]] ];
	filtm = sFilterG[{rows,cols}, sampParms[[3,3]], sampParms[[2,3]] ];
	filtq = sFilterG[{rows,cols}, sampParms[[3,4]], sampParms[[2,4]] ];
	sigm =  sampParms[[4,3]] ;
	sigp =  sampParms[[4,4]] ;
	
	gm = sampParms[[4,2]];
	gp = sampParms[[4,1]];
        visLinModelE1[te, fcontrSeq0, fcontrSeq1,
          coneI, ac,ah,ap, filth,filtm,filtq, sigm,sigp, gm,gp]
  ];
(*-------------------------------------------------*)
	
visLinModelE1[te_, fSeq0_, fSeq1_,
  coneI_, ac_,ah_,ap_, filth_,filtm_,filtq_, sigm_,sigp_, gm_,gp_] :=
	Module[{Em,Ep,
	  in,cone0,cone1,horiz,bipolar0,magno0,parvo0,parvo1},
        (* IN:  te = time for filters to wind down after image sequence *)
	(* fSeq = blurred sequence in spatial Fourier domain *)
        (*   coneI = initial cone state in above domain *)
        (*      ac,ap = temporal filter parameters *)
        (*      filtm,filtq = magno,mask spatial filters *)
        (*      sigm,sigp = magno,parvo masking contrast cutoffs *)
        (*      gm,gp = gain parameters *)
        (* OUT: cumulated magno and parvo energies: {Em, Ep} *)

	cone0 = coneI; cone1 = cone0; horiz = cone0;
	parvo0 = 0.*cone0; parvo1 = parvo0;
	Em = 0.; Ep = 0.;
        For[t=0, Less[t,Dimensions[fSeq0][[1]]], t++,
	  cone0 = ac*cone0+(1.-ac)*fSeq0[[t+1]];
	  cone1 = ac*cone1+(1.-ac)*fSeq1[[t+1]];
	  horiz = ah*horiz+(1.-ah)*Re[Fourier[cone0*filth]];
	  bipolar0 =  Re[Fourier[cone0] ]/horiz - 1.;
	  magno0 = fFilt[bipolar0,filtm];
	  Em+=sigE[(magno0-fFilt[Re[Fourier[cone1]]/horiz-1.,filtm])/
	     Sqrt[sigm*sigm+fFilt[Sqr[magno0],filtq]]];
	  parvo0 = ap*parvo0+(1.-ap)*bipolar0;
	  parvo1 = ap*parvo1+(1.-ap)*bipolar1;
	  Ep += sigE[(parvo0-parvo1)/
	     Sqrt[sigp*sigp+fFilt[Sqr[parvo0],filtq] ] ];
	]; (* end sequence For loop *)
        For[t=0, Less[t,te], t++,
	  cone0 = ac*cone0;
	  cone1 = ac*cone1;
	  horizontal = ah*horizontal+(1.-ah)*Re[Fourier[cone0*filth]];
	  bipolar0 =  Re[Fourier[cone0] ]/horizontal - 1.;
	  magno0 = fFilt[bipolar0,filtm];
	  Em += sigE[(magno0-fFilt[Re[Fourier[cone1]]/horizontal-1.])/
	     Sqrt[sigm*sigm+fFilt[Sqr[magno0],filtq] ] ] ;
	  parvo0 = ap*parvo0+(1.-ap)*bipolar0;
	  parvo1 = ap*parvo1+(1.-ap)*bipolar1;
	  Ep += sigE[(parvo0-parvo1)/
	     Sqrt[sigp*sigp+fFilt[Sqr[parvo0],filtq] ] ];
	]; (* end wind down For loop *)
	{ gm*gm*sigm*sigm*Em, gp*gp*sigp*sigp*Ep }
	];
(*-------------------------------------------------*)
	

visThreshModelE[contrSeq_, visParms_, displayParms_] :=
	Module[{length,rows,cols, sampParms, te, fSeq,
	 ac,ap, filtc,filtm, gm,gp},
			
	(* convert to sampled frequencies: *)
	{length,rows,cols} = Dimensions[sequence0];
	sampParms = {	visParms[[1]]*displayParms[[1]], 
			visParms[[2]]*rows/displayParms[[2]],
			visParms[[2]]*cols/displayParms[[3]], 
			visParms[[3]]};
	
	(* process the background sequence *)
	filtc = sFilterG[{rows,cols}, sampParms[[3,1]], sampParms[[2,1]] ];
	fSeq = (filtc*InverseFourier[#])& /@ contrSeq;

	ac = tFiltPar[sampParms[[1,1]], displayParms[[1]] ];
	ap = tFiltPar[sampParms[[1,3]], displayParms[[1]] ];

	te = 1+Floor[Log[1./E//N]/Log[ap] ]; (* ensures drop to 1./E *)

	filtm = sFilterG[{rows,cols}, sampParms[[3,3]], sampParms[[2,3]] ];
	
	gm = sampParms[[4,2]];
	gp = sampParms[[4,1]];
	visThreshModelE1[te, fSeq, ac, ap, filtm, gm, gp] 
  ];
	
visThreshModelE1[te_, fcontrSeq_,
   ac_, ap_, filtm_, gm_, gp_] :=
	Module[{bipolar, parvo, Em,Ep},
        (* IN:  te = time for filters to wind down after image sequence *)
    (* fcontrSeq = blurred contrast sequence in spatial Fourier domain *)
        (*      ac,ap = temporal filter parameters *)
        (*      filtm = magno spatial filter *)
        (*      gm,gp = gain parameters *)
        (* OUT: cumulated magno and parvo energies: {Em, Ep} *)

	bipolar = 0.*fcontrSeq[[1]]; parvo = bipolar;
	Em = 0.; Ep = 0.;
        For[t=0, Less[t,Dimensions[fcontrSeq][[1]]], t++,
	  bipolar = ac*bipolar+(1.-ac)* fcontrSeq[[t+1]];
	  Em += sigE[ bipolar*filtm ] ;
	  parvo = ap*parvo+(1.-ap)*bipolar;
	  Ep += sigE[parvo];
	]; (* end sequence For loop *)
        For[t=0, Less[t,te], t++,
	  bipolar = ab*bipolar;
	  Em += sigE[ bipolar*filtm ] ;
	  parvo = ap*parvo+(1.-ap)*bipolar;
	  Ep += sigE[parvo];
	]; (* end wind down For loop *)
	{ gm*gm*Em, gp*gp*Ep }
	];
	
	
	
(*-------------------------------------------------
</pre>
*)