(*
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 } ]; (*-------------------------------------------------*)