Creation mask for Spike Detection Index

[yup]

Hi folk!
Please help write script which create mask for Spike Detection Index.
Mathematics very simple:
p value luma in current frame
f value luma in next frame
b value luma in previous frame
need calculate
d1=abs(p-f)
d2=abs(p-b)
SDI=1-abs((d1-d2)/(d1+d2)) for d1>t1 or d2>t1 otherwise SDI=0
last operation if SDI>t2 then SDI=1
We can use two thresholds t1 and t2.
May be somebody write script for SDI, please advice.
With kind regards yup.

[gzarkadas]

If you don't need to do it in a single pass, you can:

use a simple WriteFile(... AverageLuma ...) to pass on to ScriptClip() in the first script,
import the log file to a spreadsheet and calculate SDI for each frame, based on the math-specs you supply,
use the results to make a script that masks the specific frames where SDI falls within the criteria you specify.

Especially if your video has many frames and spikes appear rarely, this approach will be faster than doing the entire math in Avisynth.

[yup]

gzarkadas!
Thank you for reply. Now I read documetation for mask tools. I have small expirience for writing relative complex Avisynth script, but I try.
yup.

[gzarkadas]

You opt for one pass as I see :) . I believe the script below will do the job, without the need for masktools:


Function CalcSDI(float b, float p, float f, float t1, float t2) {
d1 = Abs(p - f)
d2 = Abs(p - b)
ret_sdi = (d1 > t1 || d2 > t1) ? 1.0 - Abs((d1 - d2) / (d1 + d2)) : 0.0
return ret_sdi > t2 ? 1.0 : ret_sdi
}

# source to calculate SDI for

clp = AviSource("clip-to-detect.avi")

# make three clips, prev, current, next (discard pad frames later if needed)
# sdi values for first/last 2 frames of *padded* clips will not be valid of course

global prev = clp.Trim(0, -2) + clp
global curr = clp.Trim(0, -1) + clp + clp.Trim(clp.Framecount - 1, -1)
global next = clp + clp.Trim(clp.Framecount - 2, -2)

# Because AverageLuma works only inside the runtime environment,
# we write the script as a multiline string

# put your specific thresholds here

ScriptClip(curr, """
b = AverageLuma(prev)
p = AverageLuma(curr)
f = AverageLuma(next)
sdi = CalcSDI(b, p, f, 15, 0.5)
pad = " "
WriteFile(curr, "__sdi__.txt", "b", "pad", "p", "pad", "f", "pad", "sdi")
""")



The script currently is just a stub that logs sdi values; you must add lines if you want to make an action based on sdi.

Hope that this helps.

[gzarkadas]

Well, after little research I ended up with a simpler (and a bit faster) version. Here it is:


Function CalcSDI2(float d1, float d2, float t1, float t2) {
sdi = (d1 > t1 || d2 > t1) ? 1.0 - Abs((d1 - d2) / (d1 + d2)) : 0.0
return sdi > t2 ? 1.0 : sdi
}

# source to calculate SDI for
global clp = AviSource("clip-to-detect.avi")

# put your specific thresholds here
ScriptClip(clp, """
d1 = YDifferenceToNext()
d2 = YDifferenceFromPrevious()
sdi = CalcSDI2(d1, d2, 15, 0.5)
pad = " "
WriteFile(clp, "__sdi__.txt", "current_frame", "pad", "d1", "pad", "d2", "pad", "sdi")
""")

[yup]

Hi gzarkadas!
Thank you for help :thanks: . I do not think that I get solution my problem fast. I explain my aim. I have noisy video analog capture, noise mainly black and white line. I try many plugin, now best result give ml3dex with motion compensation based on MVTools, but noise introduce error in motion vector estimation and ml3dex work not at full power. I try use prefilter before MVTools, but as my noise have impulse character I need use median filter with big radius 2-3 pixel, this filter blured image and get not true motion vectors (silent circle). I want use SDI for mask and prefiler only pixels with high SDI, other pixels will be as in source untouched. After i call MVAnalyse for calculating motion vectors, make compensation and call ml3dex with source. And now i will be try write script based on mt_masktools and you made part work, thanks one more.
yup.

[gzarkadas]

Sorry, I didn't figure you were aiming for pixels from the begining. Anyway, I tried to think in reverse :D and came up with the following solution:

LoadPlugin("mt_masktools.dll") # v2.0a30

# this is for Abs(p - f) or Abs(p - b)
Function absdiff(clip c1, clip c2) { return mt_lutxy(c1, c2, "x y - abs") }

# this build the expression for the SDI function
Function SDI_RPN(string expr_d1, string expr_d2, string expr_t1, string expr_t2) {
_d1 = expr_d1 + " "
_d2 = expr_d2 + " "
_t1 = expr_t1 + " "
_t2 = expr_t2 + " "
return _d1 + _t1 + "> " + _d2 + _t1 + "> | 1.0 " \
+ _d1 + _d2 + "- " + _d1 + _d2 + "+ / abs - 0.0 ? " \
+ _t2 + "> 1.0 " + _d1 + _t1 + "> " + _d2 + _t1 + "> | 1.0 " \
+ _d1 + _d2 + "- " + _d1 + _d2 + "+ / abs - 0.0 ? "
}

# this creates the mask
Function SDIPixelMask(clip clp, float t1, float t2) {
prev = clp + clp.Trim(clp.Framecount - 2, -2)
curr = clp.Trim(0, -1) + clp + clp.Trim(clp.Framecount - 1, -1)
next = clp.Trim(0, -2) + clp

c_d1 = absdiff(curr, next)
c_d2 = absdiff(curr, prev)

# since SDI returns 0..1 we multiply with 255

mask = mt_lutxy(c_d1, c_d2, SDI_RPN("x", "y", String(t1), String(t2)) + " 255 *")
return mask.Trim(1, -clp.Framecount)
}

# this tests the implementation

clp = AviSource("clip-to-detect.avi")
SDIPixelMask(clp, 128, 0.75) # define your thresholds here


It seems to generally work, but of course it must be verified on your sources, because I haven't something similar to directly experiment with.

[yup]

Hi gzarkadas!
One more :thanks: today I try Your function and make report.
yup.

[yup]

Hi gzarkadas!
I back to this problem after more than one year period. I try make mask using motion compensated frame and not very glad. Please se my script:
LoadPlugin("mt_masktools.dll")
LoadPlugin("VerticalCleanerS.dll")
LoadPlugin("RemoveGrainS.dll")


# this is for Abs(p - f) or Abs(p - b)
Function absdiff(clip c1, clip c2) { return mt_lutxy(c1, c2, "x y - abs") }

# this build the expression for the SDI function
Function SDI_RPN(string expr_d1, string expr_d2, string expr_t1, string expr_t2) {
_d1 = expr_d1 + " "
_d2 = expr_d2 + " "
_t1 = expr_t1 + " "
_t2 = expr_t2 + " "
return _d1 + _t1 + "> " + _d2 + _t1 + "> | 1.0 " \
+ _d1 + _d2 + "- " + _d1 + _d2 + "+ / abs - 0.0 ? " \
+ _t2 + "> 1.0 " + _d1 + _t1 + "> " + _d2 + _t1 + "> | 1.0 " \
+ _d1 + _d2 + "- " + _d1 + _d2 + "+ / abs - 0.0 ? "
}

# this creates the mask
Function SDIPixelMask(clip clp, float t1, float t2) {
filtered=clp.VerticalCleaner(mode=2).RemoveGrain(11)
bv1 = filtered.MVAnalyse(blksize=8, isb = true, truemotion=true, search=2, searchparam=7, delta = 1, overlap=4, dct=0, idx = 1)
fv1 = filtered.MVAnalyse(blksize=8, isb = false, truemotion=true, search=2, searchparam=7, delta = 1, overlap=4, dct=0, idx = 1)
curr=clp
prev =clp.MVCompensate(bv1, idx=2, thSAD=16000)
next = clp.MVCompensate(fv1, idx=2, thSAD=16000)
c_d1 = absdiff(curr, next)
c_d2 = absdiff(curr, prev)

# since SDI returns 0..1 we multiply with 255

mask = mt_lutxy(c_d1, c_d2, SDI_RPN("x", "y", String(t1), String(t2)) + " 255 *")
return mask.Trim(1, -clp.Framecount)
}

# this tests the implementation

AviSource("selnew.avi")
AssumeTFF()
ConvertToYV12(interlaced=true)
SeparateFields()
clp=SelectEven()
sdi=SDIPixelMask(clp, 64, 0.8) # define your thresholds here
StackVertical(clp,sdi)
And advice. I see time shifting beween fields and mask.
Current
http://img234.imageshack.us/img234/6955/currvw5.th.png (http://img234.imageshack.us/my.php?image=currvw5.png)
previous
http://img234.imageshack.us/img234/4145/prevlg0.th.png (http://img234.imageshack.us/my.php?image=prevlg0.png)
This script work only for even fields. Now I write script which remove black lines from VHS capture video http://forum.doom9.org/showthread.php?p=1151994#post1151994, but this script introduce some artifact especialy on artificial objects (buildings, windows and etc.)
With kind regards yup.

[yup]

One more Hii to all!
Folowing code simpler and give better result
LoadPlugin("mt_masktools-25.dll")
LoadPlugin("RemoveGrainSSE2.dll")
LoadPlugin("VerticalCleanerSSE2.dll")
AVISource("sel97.avi")
AssumeTFF()
ConvertToYV12(interlaced=true)
fields=SeparateFields()
fieldsf=fields.VerticalCleaner(mode=2).RemoveGRain(11)
bv1 = fieldsf.MVAnalyse(blksize=8, isb = true, truemotion=true, search=2, delta = 2, idx = 1, overlap=4, dct=0,chroma=false)
fv1 = fieldsf.MVAnalyse(blksize=8, isb = false, truemotion=true, search=2, delta = 2, idx = 1, overlap=4, dct=0,chroma=false)
bc1=fields.MVCompensate(bv1,idx=2,thSAD=16000)
fc1=fields.MVCompensate(fv1,idx=2,thSAD=16000)
mb1=mt_lutxy(fields,bc1,"x y - abs",u=-128,v=-128)
mf1=mt_lutxy(fields,fc1,"x y - abs",u=-128,v=-128)
masksdi=mt_logic(mf1,mb1,"min",u=-128,v=-128)
StackVertical(fields,masksdi)
need only tune threshold for masks mb1 and mf1.
yup.

[yup]

Hi all!
I back to this problem
source=AVISource("sel97.avi")
source=source.AssumeTFF().ConvertToYV12(interlaced=true)
bobnn=source.nnedi3(field=-2,U=false, V=false)
bobnnmed=bobnn.mt_luts(bobnn,mode="median",pixels="0 -2 0 -1 0 0 0 1 0 2",U=2,V=2)
bobnnf=bobnn.dfttest(tbsize=1,ftype=1,sbsize=16,sosize=12,sigma=1000,U=false, V=false)
super=MSuper(bobnn,chroma=false)
superf=MSuper(bobnnf,chroma=false)
bw1 = MAnalyse(superf, blksize=16, isb = true, delta = 2, overlap=8, dct=5,chroma=false)
fw1 = MAnalyse(superf, blksize=16, isb = false,delta = 2, overlap=8, dct=5,chroma=false)
bc1 = MCompensate(bobnn, super, bw1, thSAD=16000, thSCD1=16000)
fc1 = MCompensate(bobnn, super, fw1, thSAD=16000, thSCD1=16000)
bwabs=mt_lutxy(bobnn,bc1,"x y - abs")
fwabs=mt_lutxy(bobnn,fc1,"x y - abs")
thresh=25
TH=string(thresh)
SDI=mt_lutxy(bwabs,fwabs,"x "+TH+" > y "+TH+" > | 1 x y - x y + / abs - 255 * 0 ? ", u=-128, v=-128)
SDIavg=SDI.mt_luts(SDI,mode="avg",pixels=mt_square(3),U=2,V=2)
SDIad=mt_lutxy(SDI,SDIavg,"x 2 y * > x 128 > & 255 0 ?")
SDIad=SDIad.mt_luts(SDIad,mode="median",pixels="0 -1 0 0 0 1",U=-128,V=-128)
fieldSDI=SDIad.AssumeTFF().SeparateFields().SelectEvery(4,0,3)
StackVertical(Separatefields(source),fieldSDI)


http://img12.imageshack.us/img12/6620/sdi.th.png (http://img12.imageshack.us/i/sdi.png/)

This script based on http://www.mee.tcd.ie/~ack/papers/a4ackphd.ps.gz Chapters 6.3-6.3.2.
Temporal threshold for spike 25, for spatial averaging 2 (SDI could be greater 2 mean value).
Median postfilter for remove artifacts during motion compensation (size for this postfilter could be less than for prefilter before motion estimation).
yup.

[yup]

Hi all!
I add block calculation for primary SDI (averaging with overlap) . Also inpand expand approach for remove false small spike. Give little less false alarm.
source=AVISource("selnew.avi")
source=source.AssumeTFF().ConvertToYV12(interlaced=true)
bobnn=source.nnedi3(field=-2,U=false, V=false)
bobnnmed=bobnn.mt_luts(bobnn,mode="median",pixels="0 -3 0 -2 0 -1 0 0 0 1 0 2 0 3",U=2,V=2)
bobnnf=bobnnmed.dfttest(tbsize=1,ftype=1,sbsize=16,sosize=12,sigma=1000,U=false, V=false)
super=MSuper(bobnn,chroma=false)
superf=MSuper(bobnnf,chroma=false)
bw1 = MAnalyse(superf, blksize=8, isb = true, delta = 2, overlap=4, dct=5,chroma=false)
fw1 = MAnalyse(superf, blksize=8, isb = false,delta = 2, overlap=4, dct=5,chroma=false)
bc1 = MCompensate(bobnn, super, bw1, thSAD=16000, thSCD1=16000)
fc1 = MCompensate(bobnn, super, fw1, thSAD=16000, thSCD1=16000)
bwabs=mt_lutxy(bobnn,bc1,"x y - abs",U=-128,V=-128)
fwabs=mt_lutxy(bobnn,fc1,"x y - abs",U=-128,V=-128)
threshsp=10
THP=string(threshsp)
threshavg=2
THAVG=string(threshavg)
SDI=mt_lutxy(bwabs,fwabs,"x "+THP+" > y "+THP+" > | 1 x y - x y + / abs - 255 * 0 ? ", u=-128, v=-128)
SDIavg=SDI.mt_luts(SDI,mode="avg",pixels=mt_square(5),U=-128,V=-128)
SDIAvgblk=SDIavg.PointResize(90,72)
SDIavgblkX=SDIAvgblk.mt_luts(SDIAvgblk,mode="median",pixels="1 -1 0 0 -1 1 1 1 -1 -1",U=-128,V=-128)
SDIavgblkCross=SDIAvgblk.mt_luts(SDIAvgblk,mode="median",pixels="1 0 0 0 -1 0 0 1 0 -1",U=-128,V=-128)
SDIAvgblk=clense(SDIAvgblk,SDIavgblkCross,SDIavgblkX,increment=0, grey=true).mt_lut("x",U=-128,V=-128)
SDIavg=SDIAvgblk.PointResize(720,576)
SDIad=mt_lutxy(SDI,SDIavg,"x "+THAVG+" y * > x 128 > & 255 0 ?",U=-128,V=-128)
SDIad=SDIad.mt_luts(SDIad,mode="median",pixels="0 -2 0 -1 0 0 0 1 0 2",U=-128,V=-128)
SDIad=SDIad.mt_inpand(mode="horizontal",U=-128,V=-128).mt_inpand(mode="horizontal",U=-128,V=-128).mt_expand(mode="horizontal",U=-128,V=-128).mt_expand(mode="horizontal",U=-128,V=-128)
fieldSDI=SDIad.AssumeTFF().SeparateFields().SelectEvery(4,0,3)
StackVertical(Separatefields(source),fieldSDI)
yup.