Find Temporal Statistics

[jmac698]

Hi,
I'm a bit stuck on an analysis I want to write.
I have a reference frame, and want to find the min,max, median, and average differences of n frames to it - temporally.

Let's focus on one pixel to understand this. Let's say ref at x=0 and y=0 is r,g,b=(100,100,100). I want to compare this to a clip, with the corresponding pixel (0,0) at frame n. I can compare just two frames with subtract(ref,clip.trim(n,-1)). But I want to compare all the clip's frames to the ref, so now I have n frames with signed numbers (signed offset from subtract). Now I want to find the abs() min, max, avg, median per pixel of all these n frames temporally. In other words, the temporal noise per pixel. The output should be 4 frames representing that statistic for each type.

I could use Stainless clipblend(0) for the avg, I'm not sure how to do a large temporal median, and I don't know how to find temporal min and max. However, a quick mod to clipblend could return min, max,avg easily enough.

Thx

ps this wouldn't be hard in vapoursynth except you can't just read pixels, the raw frame is exposed and you have to know the pixel layouts for each YV12 etc format, that's my biggest barrier now to pixel level processing.

pps
I thought of a quick way to do approximate median for any sized clip, keep a short histogram of the pixels, for each pixel set 2^n (limited to -16,+15) in an int, so your histogram is 32 bits per pixel. To find the median, do a population count of bits set, then find the pop/2th bit set, and that's your median. Add in scaling to bins up to size 16. Example:

at (0,0), red channel only:
ref=128
frame0=120 (-8)
frame1=131 (+3)
frame2=135 (+7)
step1: subtract 112 from each value (which represents -16 from ref):
frame0=8
frame1=19
frame2=23
step2: set those bits in a 32bit int (by ORing)
2^8+2^19+2^23=00000000 10001000 00000001 00000000
population count=3
find floor(3/2)th bit set, or 1
bitscan... bit 19 is set
find 19-16=3
the median is 3
you need population count because the same bit might be set multiple times, could be 1 to n bits actually set.
The median of even samples is the average of the middle two.
You can use scaling to set a bit if it's in any bucket of a range of 16 values.
I think popcnt is asm, I don't know about finding the nth bit.
If the value is out of range of the histogram, set the nearest value.

This is to analyze hot pixels in cameras. They aren't always just stuck on full, but biased relative to normal pixels.

I could turn this into a spatial problem with weaverows, but then how to find min,max etc along columns?

[johnmeyer]

You should have started with "This is to analyze hot pixels in cameras." I say that because after having skimmed through the main part of your post, I'm not sure your approach is going to yield anything.

I just finished two 4-camera wedding shoots, and the three cameras that were locked down have video that hardly varies from frame to frame. A hot pixel wouldn't look any different than a sparkle on a non-moving dress. Even if you averaged five minutes of footage from my videos, 98% of the pixels would not change, and therefore temporal averages wouldn't mean much.

I think that you inevitably have to introduce some sort way to only do averaging on pixels that are changing from frame to frame, i.e., where there is significant motion.

[Reel.Deel]

@jmac698

I don't know if this will be of any help but there's the HealDeadPixels (https://github.com/ladipro/avisynth_filters/wiki/HealDeadPixels) plugin that might be worth a look.

[jmac698]

Thanks for the input, however you misunderstand me. This is not for normal pictures, but astrophotography and night photography (with minutes long exposures) where noise becomes very important. It is normal to use 3 kinds of calibration frames; bias frames (those with a very short exposure, to simply read the patterns in the row amplifiers), dark frames (which shows the slight variations in individual pixels, due to manufacturing variations), and flat frames (normal pics of a white wall to correct for vignetting). My research involves going even further, because noise depends on temperature as well, and for this reason, people usually have to take the dark frames at the same time as the normal frames (during long exposures), obviously this more than doubles the time, and you don't have that much time to waste, for example a picture of Jupiter which rotates visibly in under 1 minute, or the moon which can move out of sight as you watch (without a fancy tracking mount).

What I think I can do is pre-make calibration frames, then apply them based on temperature. There's some other things I can investigate too. There's poisson and gaussian noise, but poisson varies based on image intensity AFAIK and has to be denoised differently. I wanted to find the cross over point.

So you see, this is all very much more advanced than you might typically think. Also my camera has no hotpixels at all. There is some bias patterns however.

I have a test here in low light, iso1600, it looks like complete coloured noise, yet after averaging 12 frames, I can clearly see an image.

Another application is even for the day time, with a deep IR filter an unmodified camera, you can still take IR pictures, but they are very dim. In fact, it can replace an ND filter to look at the sun.

[StainlessS]

Jmac, what is maximum required temporal radius, (presumably both forward and backwards) ?

EDIT: ClipBlend() averages frames prior to current (with Delay=0 being all frames).
RGBAmplifier() samples with temporal radius (up to 99 frames, I think).
Either could probably be modified to produce required stats (but you would need to decide what constitutes min or median etc).

[StainlessS]

Jmac,
Gonna attempt to create file based count array, 256*sizeof(unsigned int) for each pixel/channel (RGB only), then you can
apply whatever method you want to extract data (defo need this for median, I think).

EDIT: Will likely not be terribly speedy.

[StainlessS]

Jmac, here ClipStats:- http://www.mediafire.com/download/1yeq7q02mr0ar2l/ClipStats-v0.0-20151004.zip

EDIT: Above link edited (got filename month and day back to front).

Slow as hell (about 20 secs per frame [after initial file created, on 2.4Ghz Core Duo]).


/*
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

*/

#include "Compiler.h"

#include <windows.h>
#include <stdio.h>
#include <math.h>

#include "avisynth.h"

#define BUG

#ifdef BUG
int dprintf(char* fmt, ...) {
char printString[512]=""; // MUST NUL TERM THIS
char *p=printString;
va_list argp;
va_start(argp, fmt);
vsprintf(p, fmt, argp);
va_end(argp);
for(;*p++;);
--p; // @ null term
if(printString == p || p[-1] != '\n') {
p[0]='\n'; // append n/l if not there already
p[1]='\0';
}
OutputDebugString(printString);
return p-printString; // strlen printString
}

// OLD C89 Compiler: (eg VS6)
// Call as eg DPRINTF(("Forty Two = %d",42)) // Enclosed in double parentethis, No Semi colon
#define DPRINTF(x) dprintf x;

// C99 Compiler (eg VS2008)
// #define DPRINTF(fmt, ...) dprintf(fmt, ##__VA_ARGS__); // No Semi colon needed
#else
// OLD C89 Compiler: (eg VS6)
#define DPRINTF(x); /* x */

// #define DPRINTF(fmt,...); /* fmt */
#endif


class ClipStats : public GenericVideoFilter {
const int delay;
const int mode;
const double thresh;
const int lo;
const int hi;
const char *fn;
//
FILE *fp;
int num_frames;
int Els_0;
int h_y;
int w_y;
int FramesDone;
int last;

enum {MODE_MIN=0,MODE_MAX,MODE_MINMAX,MODE_MEDIAN,MODE_AVG,MODE_STDEV,MODE_IRNG,MODE_LIMIT=MODE_IRNG};

void Reset(void);
void Add(int n, IScriptEnvironment* env);
void Sub(int n, IScriptEnvironment* env);
PVideoFrame Render(int n, IScriptEnvironment* env);

public:
ClipStats(PClip _child,int _delay,int _mode,double _thresh,int _lo,int _hi,const char *_fn,IScriptEnvironment* env);
~ClipStats();
PVideoFrame __stdcall GetFrame(int n, IScriptEnvironment* env);
};

ClipStats::~ClipStats() {
if(fp!=NULL) {fclose(fp); fp=NULL;};
}

ClipStats::ClipStats(PClip _child,int _delay,int _mode,double _thresh,int _lo,int _hi,const char *_fn,IScriptEnvironment* env) :
GenericVideoFilter(_child),delay(_delay),mode(_mode),thresh(_thresh),lo(_lo),hi(_hi),fn(_fn),fp(NULL) {


if(vi.width==0 || vi.num_frames<=0) env->ThrowError("ClipStats: Clip has no video\n");
if(!vi.IsRGB24()) env->ThrowError("ClipStats: RGB24 Only\n");
if(delay<0) env->ThrowError("ClipStats: +ve delay only please.\n");
if(mode<0||mode>MODE_LIMIT) env->ThrowError("ClipStats: 0 <= Mode <= %d.\n",MODE_LIMIT);
if(thresh<0.0||thresh>100.0) env->ThrowError("ClipStats: 0.0 <= Thresh <= 100.0\n");
if(lo<0||lo>255) env->ThrowError("ClipStats: 0 <= Lo <= 255.\n");
if(hi<lo||hi>255) env->ThrowError("ClipStats: Lo <= Hi <= 255.\n");

num_frames = vi.num_frames;
PVideoFrame src = child->GetFrame(0, env); // get frame 0
h_y = src->GetHeight();
w_y = src->GetRowSize();
Els_0 = h_y * w_y;
if((fp=fopen(fn, "wb+" ))==NULL) env->ThrowError("ClipStats: Cannot open data file %s",fn);
FramesDone=0;
last = -1;
}


void ClipStats::Reset(void) {
unsigned int bf[256];
int i;
for(i=256;--i>=0;bf[i]=0);
rewind(fp);
for(i=Els_0;--i>=0;) {
int wr=fwrite(bf,sizeof(bf),1,fp); // should be return 1
}
FramesDone=0;
}

void ClipStats::Add(int n,IScriptEnvironment* env) {
PVideoFrame src = child->GetFrame(n, env);
const int pitch = src->GetPitch();
const BYTE * srcp = src->GetReadPtr();
++FramesDone; // counting this frame
int y,x;
__int64 pos;
for(pos=0,y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
__int64 off = (pos*256 + srcp[x]) * sizeof(unsigned int);
_fseeki64(fp, off,SEEK_SET);
unsigned int val;
fread(&val,sizeof(val),1,fp);
_fseeki64(fp,off,SEEK_SET);
++val;
fwrite(&val,sizeof(val),1,fp);
++pos;
}
srcp += pitch;
}
}

void ClipStats::Sub(int n,IScriptEnvironment* env) {
PVideoFrame src = child->GetFrame(n, env);
const int pitch = src->GetPitch();
const BYTE * srcp = src->GetReadPtr();
--FramesDone; // counting this frame
int y,x;
__int64 pos;
for(pos=0,y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
__int64 off = (pos*256 + srcp[x]) * sizeof(unsigned int);
_fseeki64(fp, off,SEEK_SET);
unsigned int val;
fread(&val,sizeof(val),1,fp);
_fseeki64(fp, off,SEEK_SET);
--val;
fwrite(&val,sizeof(val),1,fp);
++pos;
}
srcp += pitch;
}
}


PVideoFrame ClipStats::Render(int n,IScriptEnvironment* env) {
PVideoFrame dst = env->NewVideoFrame(vi);
const int dpitch = dst->GetPitch();
BYTE * dstp= dst->GetWritePtr();
int y,x,i;
unsigned int cnt[256];
rewind(fp);
switch(mode) {
case MODE_MIN:
{
const unsigned int lim = (unsigned int)(FramesDone * thresh);
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmin=255;
unsigned int sm=0;
for(i=0;i<256;++i) {
sm += cnt[i];
if(sm> lim) {
rmin=i;
break;
}
}
dstp[x]=rmin;
}
dstp += dpitch;
}
}
break;
case MODE_MAX:
{
const unsigned int lim = (unsigned int)(FramesDone * thresh);
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmax=0;
unsigned int sm=0;
for(i=256;--i>=0;) {
sm += cnt[i];
if(sm > lim) {
rmax=i;
break;
}
}
dstp[x]=rmax;
}
dstp += dpitch;
}
}
break;
case MODE_MINMAX:
{
const unsigned int lim = (unsigned int)(FramesDone * thresh);
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmin=255;
unsigned int sm=0;
for(i=0;i<256;++i) {
sm += cnt[i];
if(sm> lim) {
rmin=i;
break;
}
}
int rmax=0;
sm=0;
for(i=256;--i>=0;) {
sm += cnt[i];
if(sm > lim) {
rmax=i;
break;
}
}
dstp[x]=rmax-rmin;
}
dstp += dpitch;
}
}
break;
case MODE_MEDIAN:
{
const unsigned int lim = FramesDone / 2;
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmed=0;
unsigned int sm=0;
for(i=256;--i>=0;) {
sm += cnt[i];
if(sm > lim) {
rmed=i;
break;
}
}
dstp[x]=rmed;
}
dstp += dpitch;
}
}
break;
case MODE_AVG:
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
double mean=0.0;
__int64 acc = 0;
for(i=256;--i>=0;) {
acc += cnt[i] * i;
}
mean = (double)acc / FramesDone;
dstp[x]=int(mean+0.5);
}
dstp += dpitch;
}
break;
case MODE_STDEV:
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
double mean=0.0;
double std=0.0;
if(FramesDone>1) {
__int64 acc = 0;
for(i=256;--i>=0;) {
acc += cnt[i] * i;
}
mean = (double)acc / FramesDone;
double Sum_DiffSquared = 0.0;
for (i=256;--i>=0;) {
const unsigned int count=cnt[i];
if(count) {
const double diff = (double)i - mean; // x - xbar
Sum_DiffSquared += (diff * diff * double(count)); // Sigma((x-xbar)^2)
}
}
std=(Sum_DiffSquared<=0.0)?0.0:sqrt(Sum_DiffSquared / FramesDone);
}
dstp[x]=int(std+0.5);
}
dstp += dpitch;
}
break;
case MODE_IRNG:
for(y=0;y<h_y;++y) {
for(x=0;x<w_y;++x) {
fread(cnt,sizeof(cnt),1,fp);
unsigned int sm=0;
for(i=lo;i<=hi;++i) {
sm += cnt[i];
}
dstp[x]= int((double(sm) / FramesDone)+0.5);
}
dstp += dpitch;
}
break;
}
return dst;
}




PVideoFrame __stdcall ClipStats::GetFrame(int n, IScriptEnvironment* env) {
n = (n<0) ? 0 : (n>= num_frames) ? num_frames - 1 : n;

PVideoFrame dst;

DPRINTF(("ClipStats: %d] Last=%d FramesDone=%d",n,last,FramesDone))

if(n != 0 && n == last) { // Assume cache failure

DPRINTF(("ClipStats: %d] Cache Fail (repeat request for this frame)",n))
dst = Render(n,env);
} else {
int WantS = (delay==0) ? 0 : max(n - delay,0);
int WantE = n;
int GotS = (FramesDone>0) ? max(last - FramesDone + 1,0) : - 1;
int GotE = (FramesDone>0) ? last : - 1;
int ResCnt = WantE-WantS+1; // frames to read if accumulator reset
int AccCnt= (n<last) ? GotE-WantE + GotS-WantS : WantE-GotE + WantS-GotS; // frames to read if subtract/add accumlator

DPRINTF(("ClipStats: GotS=%d GotE=%d WantS=%d WantE=%d",GotS,GotE,WantS,WantE))
// Frame 0 or jump about, restart:clear Accumulators
if(n==0 || FramesDone<=0 || ResCnt<AccCnt) {
DPRINTF(("ClipStats: %d] RESET (last=%d ResCnt=%d AccCnt=%d)",n,last,ResCnt,AccCnt))
Reset();
GotS = - 1;
GotE = - 1;
} else if(n > last) {
while(GotS<WantS) { // Deduct oldest frames dat from Accumulator
DPRINTF(("ClipStats: DEDUCT_S: %d",GotS))
Sub(GotS,env);
++GotS;
}
} else {
int m=WantE+1;
while(m<=GotE) {
DPRINTF(("ClipStats: DEDUCT_E: %d",m))
Sub(n,env);
++m;
}
GotE=WantE;
}

if(WantS<GotS) {
int m=WantS;
while(m<GotS) {
DPRINTF(("ClipStats: ADD_S: %d (GotS=%d WantS=%d",m,GotS,WantS))
Add(m,env);
++m;
}
GotS=WantS;
}

if(GotE<0) {
GotS = WantS ;
GotE = WantS - 1;
}
if(GotE < WantE) {
while(GotE<WantE) {
++GotE;
if(GotE!=WantE) {
DPRINTF(("ClipStats: ADD_E: %d",GotE))
Add(GotE,env);
} else {
DPRINTF(("ClipStats: ADD_RENDER: %d",GotE))
Add(GotE,env);
dst=Render(GotE,env);
}
}
} else {

DPRINTF(("ClipStats: RENDER: %d",n))
dst=Render(n,env);
}
last = n;
}
DPRINTF(("ClipStats:"))
return dst;
}


AVSValue __cdecl Create_ClipStats(AVSValue args, void* user_data, IScriptEnvironment* env) {
int lo = args[4].AsInt(128); // lo
AVSValue ret = new ClipStats(
args[0].AsClip(),
args[1].AsInt(0), // delay
args[2].AsInt(0), // Mode
args[3].AsFloat(0.0), // Thresh
lo, // lo
args[5].AsInt(lo), // hi
args[6].AsString("ClipStats.Data"), // Filename
env);
return ret;
}

extern "C" __declspec(dllexport) const char* __stdcall AvisynthPluginInit2(IScriptEnvironment* env) {
env->AddFunction("ClipStats", "c[delay]i[Mode]i[Thresh]f[Lo]i[Hi]i[FN]s", Create_ClipStats, 0);

return "`ClipStats' ClipStats plugin";
// A freeform name of the plugin.
}


Test avs


avisource("E:\V\JurassicPark.avi")
ConvertToRGB24
Trim(3000,0)
/*
RGB24 Only.

Clipstats(clip c,int "delay"=0,int "Mode"=0,Float "Thresh"=0.0,Int "lo"=128,Int "Hi"=lo,String "FN"="ClipStats.Data")
Delay:- as in ClipBlend (0=all frames).
Mode:- 0=Min, 1=Max, 2=MinMaxDifference, 3=Median, 4=Avg, 5=Stdev, 6=Inrange
Thresh:- As In YPlaneMin(Threshold) [but of course for R,G,B]
Lo:- As in RT_YInRange.
Hi:- As in RT_YInRange.

*/


A=ClipStats(delay=0,mode=4)
B=ClipBlend()

StackHorizontal(A,B)



Little testing done, [only comparison with Clipblend(delay=0)]

EDIT: Supplied with VS 2008 build files.

EDIT:
I have a reference frame, and want to find the min,max, median, and average differences of n frames to it - temporally.

Oh dear, I should not skim read, not exactly as you require, but could be modded to suit now that its working.

[StainlessS]

Take 4, http://www.mediafire.com/download/cj06ciacuhfq6rc/ClipDiffs-v0.0-TAKE4-20151005.zip

Minus GPL, headers, debug stuff due to size.

class ClipDiffs : public GenericVideoFilter {
PClip child2;
const int delay;
const int mode;
const double thresh;
const int lo;
const int hi;
const char *fn;
//
FILE *fp;
int num_frames;
int Els;
int hit;
int wid;
int FramesDone;
int last;
BYTE Diff[513];

enum {MODE_MIN=0,MODE_MAX,MODE_MINMAX,MODE_MEDIAN,MODE_AVG,MODE_STDEV,MODE_IRNG,MODE_LIMIT=MODE_IRNG};

enum{CNT_BFSZ=512}; // Last one not used [-255 -> 255 == 511]
enum{CNT_REL =255};

void Reset(void);
void Add(int n, IScriptEnvironment* env);
void Sub(int n, IScriptEnvironment* env);
PVideoFrame Render(int n, IScriptEnvironment* env);

public:
ClipDiffs(PClip _child,PClip _child2,int _delay,int _mode,double _thresh,int _lo,int _hi,const char *_fn,IScriptEnvironment* env);
~ClipDiffs();
PVideoFrame __stdcall GetFrame(int n, IScriptEnvironment* env);
};

ClipDiffs::~ClipDiffs() {
if(fp!=NULL) {fclose(fp); fp=NULL;};
}

ClipDiffs::ClipDiffs(PClip _child,PClip _child2,int _delay,int _mode,double _thresh,int _lo,int _hi,const char *_fn,IScriptEnvironment* env) :
GenericVideoFilter(_child),child2(_child2),delay(_delay),mode(_mode),thresh(_thresh),lo(_lo),hi(_hi),fn(_fn),fp(NULL) {
const char *myName="ClipDiffs: ";

if(vi.width==0 || vi.num_frames<=0) env->ThrowError("%sClip has no video",myName);
if(!vi.IsRGB24()) env->ThrowError("%sRGB24 Only",myName);
const VideoInfo &vi2=child2->GetVideoInfo();
if(!vi2.IsRGB24()) env->ThrowError("%sclip 2 RGB24 Only",myName);
if(vi2.width==0 || vi2.num_frames<=0) env->ThrowError("%sclip 2 has no video",myName);

if(vi.height != vi2.height) env->ThrowError("%sclip 2 height mismatch",myName);
if(vi.width != vi2.width) env->ThrowError("%sclip 2 width mismatch",myName);

if(delay<0) env->ThrowError("%s+ve delay only please.",myName);
if(mode<0||mode>MODE_LIMIT) env->ThrowError("%s0 <= Mode <= %d.",myName,MODE_LIMIT);
if(thresh<0.0||thresh>100.0) env->ThrowError("%s0.0 <= Thresh <= 100.0",myName);
if(lo< -255||lo>255) env->ThrowError("%s-255 <= Lo <= 255.",myName);
if(hi<lo||hi>255) env->ThrowError("%sLo <= Hi <= 255.",myName);

num_frames = vi.num_frames;
PVideoFrame src = child->GetFrame(0, env); // get frame 0
hit = src->GetHeight();
wid = src->GetRowSize();
Els = hit * wid;
if((fp=fopen(fn, "wb+" ))==NULL) env->ThrowError("%sCannot open data file %s",myName,fn);
int i;
for (i=0; i<=512; i++) Diff[i] = max(0,min(255,i-129));
FramesDone=0;
last = -1;
}

void ClipDiffs::Reset(void) {
int bf[CNT_BFSZ];
int i;
rewind(fp);
const int bfsz = 512*sizeof(bf); // 1MB
BYTE *bfp = new BYTE[bfsz];
if(bfp==NULL) {
memset(bf,0,sizeof(bf));
for(i=Els;--i>=0;) {
int wr=fwrite(bf,sizeof(bf),1,fp); // should be return 1
}
} else {
memset(bfp,0,bfsz);
__int64 clrbytes = __int64(sizeof(bf)) * Els;
int sods = int(clrbytes / bfsz);
int odds = int(clrbytes % bfsz);
int wr;
for(wr=1,i = sods; wr == 1 && --i>=0; ) {
wr=fwrite(bfp,bfsz,1,fp);
}
if(wr==1 && odds>0) {
wr=fwrite(bfp,odds,1,fp);
}
delete [] bfp;
}
FramesDone=0;
}

void ClipDiffs::Add(int n,IScriptEnvironment* env) {
PVideoFrame src = child->GetFrame(n, env);
const int pitch = src->GetPitch();
const BYTE * srcp = src->GetReadPtr();
PVideoFrame ref = child2->GetFrame(0,env);
const int pitchr = ref->GetPitch();
const BYTE * refp= ref->GetReadPtr();
++FramesDone; // counting this frame
int y,x;
__int64 pos;
for(pos=0,y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
int dif = srcp[x] - refp[x]; // -255 -> 255
__int64 off = (pos*CNT_BFSZ + dif + CNT_REL) * sizeof(int);
_fseeki64(fp, off,SEEK_SET);
int val;
fread(&val,sizeof(val),1,fp);
_fseeki64(fp,off,SEEK_SET);
++val;
fwrite(&val,sizeof(val),1,fp);
++pos;
}
srcp += pitch;
refp += pitchr;
}
}

void ClipDiffs::Sub(int n,IScriptEnvironment* env) {
PVideoFrame src = child->GetFrame(n, env);
const int pitch = src->GetPitch();
const BYTE * srcp = src->GetReadPtr();
PVideoFrame ref = child2->GetFrame(0,env);
const int pitchr = ref->GetPitch();
const BYTE * refp= ref->GetReadPtr();
--FramesDone; // counting this frame
int y,x;
__int64 pos;
for(pos=0,y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
int dif = srcp[x] - refp[x];
__int64 off = (pos*CNT_BFSZ + dif + CNT_REL) * sizeof(int);
_fseeki64(fp, off,SEEK_SET);
int val;
fread(&val,sizeof(val),1,fp);
_fseeki64(fp, off,SEEK_SET);
--val;
fwrite(&val,sizeof(val),1,fp);
++pos;
}
srcp += pitch;
refp += pitchr;
}
}


PVideoFrame ClipDiffs::Render(int n,IScriptEnvironment* env) {
PVideoFrame dst = env->NewVideoFrame(vi);
const int dpitch = dst->GetPitch();
BYTE * dstp= dst->GetWritePtr();
int y,x,i;
int cnt[CNT_BFSZ];
rewind(fp);
switch(mode) {
case MODE_MIN:
{
const int lim = (int)(FramesDone * thresh);
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmin=255;
int sm=0;
for(i=0;i<511;++i) {
sm += cnt[i];
if(sm> lim) {
rmin=i-CNT_REL;
break;
}
}
dstp[x]=Diff[rmin + 128 + 129]; // Min diff can be +ve or -ve, center @ 128
}
dstp += dpitch;
}
}
break;
case MODE_MAX:
{
const int lim = (int)(FramesDone * thresh);
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmax=-255;
int sm=0;
for(i=511; --i>= 0;) {
sm += cnt[i];
if(sm > lim) {
rmax=i-CNT_REL;
break;
}
}
dstp[x]=Diff[rmax + 128 + 129]; // Max diff can be +ve or -ve, center @ 128
}
dstp += dpitch;
}
}
break;
case MODE_MINMAX:
{
const int lim = (int)(FramesDone * thresh);
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmin=255;
int sm=0;
for(i=0;i<511;++i) {
sm += cnt[i];
if(sm> lim) {
rmin=i-CNT_REL;
break;
}
}
int rmax= -255;
sm=0;
for(i=511; --i>= 0;) {
sm += cnt[i];
if(sm > lim) {
rmax=i-CNT_REL;
break;
}
}
int minmxdif = rmax-rmin; // 0 -> 510 (max, 255 - -255)
minmxdif=min(minmxdif,512-129); // Dont stray outside Diff[] due to no add 128
dstp[x]=Diff[minmxdif + 129]; // Dont add 128, center a 0 not 128
}
dstp += dpitch;
}
}
break;
case MODE_MEDIAN:
{
const int lim = FramesDone / 2;
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
int rmed=-255;
int sm=0;
for(i=511; --i>= 0;) {
sm += cnt[i];
if(sm > lim) {
rmed=i-CNT_REL;
break;
}
}
dstp[x]=Diff[rmed + 128 + 129]; // Median diff can be +ve or -ve, center @ 128
}
dstp += dpitch;
}
}
break;
case MODE_AVG:
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
double mean=0.0;
__int64 acc= 0;
for(i=511; --i>= 0;) {
acc += cnt[i] * (i-CNT_REL);
}
mean = double(acc) / FramesDone;
dstp[x]=Diff[int(mean+0.5) + 128 + 129]; // Avg diff can be +ve or -ve, center @ 128
}
dstp += dpitch;
}
break;
case MODE_STDEV:
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
double mean=0.0;
double std=0.0;
if(FramesDone>1) {
__int64 acc = 0;
for(i=511; --i>= 0;) {
acc += cnt[i] * (i-CNT_REL);
}
mean = (double)acc / FramesDone;
double Sum_DiffSquared = 0.0;
for (i=511; --i>= 0;) {
const int count=cnt[i];
if(count) {
const double diff = ((double)i-CNT_REL) - mean; // x - xbar
Sum_DiffSquared += (diff * diff * double(count)); // Sigma((x-xbar)^2)
}
}
std=(Sum_DiffSquared<=0.0)?0.0:sqrt(Sum_DiffSquared / FramesDone);
}
int istd = int(std+0.5);
istd=min(istd,512-129); // Dont stray outside Diff[] due to no add 128
dstp[x]=Diff[istd + 129]; // Std +ve, Dont add 128, center a 0 not 128
}
dstp += dpitch;
}
break;
case MODE_IRNG:
for(y=0;y<hit;++y) {
for(x=0;x<wid;++x) {
fread(cnt,sizeof(cnt),1,fp);
int sm=0;
for(i=lo;i<=hi;++i) {
sm += cnt[i+CNT_REL];
}
int ir = int((sm * 255.0 / FramesDone)+0.5);
dstp[x]= ir; // InRange 0 -> 255
}
dstp += dpitch;
}
break;
}
return dst;
}




PVideoFrame __stdcall ClipDiffs::GetFrame(int n, IScriptEnvironment* env) {
n = (n<0) ? 0 : (n>= num_frames) ? num_frames - 1 : n;
PVideoFrame dst;
DPRINTF("%d] Last=%d FramesDone=%d",n,last,FramesDone)
if(n != 0 && n == last) { // Assume cache failure
DPRINTF("%d] Cache Fail (repeat request for this frame)",n)
} else {
int WantS = (delay==0) ? 0 : max(n - delay,0);
int WantE = n;
int GotS = (FramesDone>0) ? max(last - FramesDone + 1,0) : - 1;
int GotE = (FramesDone>0) ? last : - 1;
int ResCnt = WantE-WantS+1; // frames to read if accumulator reset
int AccCnt= (n<last) ? GotE-WantE + GotS-WantS : WantE-GotE + WantS-GotS; // frames to read if subtract/add accumlator

DPRINTF("%d] GotS=%d GotE=%d WantS=%d WantE=%d",n,GotS,GotE,WantS,WantE)
// Frame 0 or jump about, restart:clear Accumulators
if(n==0 || FramesDone<=0 || ResCnt<AccCnt) {
DPRINTF("%d] RESET (last=%d ResCnt=%d AccCnt=%d)",n,last,ResCnt,AccCnt)
Reset();
GotS = - 1;
GotE = - 1;
} else if(n > last) {
while(GotS<WantS) { // Deduct oldest frames dat from Accumulator
DPRINTF("%d] DEDUCT_S: %d",n,GotS)
Sub(GotS,env);
++GotS;
}
} else {
int m=WantE+1;
while(m<=GotE) {
DPRINTF("%d] DEDUCT_E: %d",n,m)
Sub(m,env);
++m;
}
GotE=WantE;
}

if(WantS<GotS) {
int m=WantS;
while(m<GotS) {
DPRINTF("%d] ADD_S: %d (GotS=%d WantS=%d",n,m,GotS,WantS)
Add(m,env);
++m;
}
GotS=WantS;
}

if(GotE<0) {
GotE = WantS - 1;
}
while(GotE<WantE) {
++GotE;
DPRINTF("%d] ADD_E: %d",n,GotE)
Add(GotE,env);
}
}
DPRINTF("%d] RENDER: %d",n,n)
dst=Render(n,env);
last = n;
DPRINTF("")
return dst;
}


AVSValue __cdecl Create_ClipDiffs(AVSValue args, void* user_data, IScriptEnvironment* env) {
int lo = args[5].AsInt(0); // lo
AVSValue ret = new ClipDiffs(
args[0].AsClip(),
args[1].AsClip(),
args[2].AsInt(0), // delay
args[3].AsInt(0), // Mode
args[4].AsFloat(0.0), // Thresh
lo, // lo
args[6].AsInt(lo), // hi
args[7].AsString("ClipDiffs.Data"), // Filename
env);
return ret;
}

extern "C" __declspec(dllexport) const char* __stdcall AvisynthPluginInit2(IScriptEnvironment* env) {
env->AddFunction("ClipDiffs", "cc[delay]i[Mode]i[Thresh]f[Lo]i[Hi]i[FN]s", Create_ClipDiffs, 0);

return "`ClipDiffs' ClipDiffs plugin";
// A freeform name of the plugin.
}

[StainlessS]

Above post updated.

test client

avisource("E:\V\JurassicPark.avi").ConvertToRGB24
BilinearResize(320,240) # reduced size (quicker)

Trim(3000,-10) # 10 frames
ORG=Last
/*
RGB24 Only.

ClipDiffs(clip c,clip Ref,int "delay"=0,int "Mode"=0,Float "Thresh"=0.0,Int "lo"=0,Int "Hi"=lo,String "FN"="ClipDiffs.Data")
Ref:- Reference clip, frame 0 ONLY used.
Delay:- as in ClipBlend (0=all frames up to current_frame, other = current_frame + previous 'Delay' frames).
Mode:- 0=Min, 1=Max, 2=MinMaxDifference, 3=Median, 4=Avg, 5=Stdev, 6=InRange
Thresh:- (0->100%) As In YPlaneMin(Threshold) [but of course for R,G,B] Modes 0, 1, and 2 only.
A value of Thresh=1.0% for eg mode=MIN would ignore up to 1% of the most -ve differences (or least +ve), ie noise.
Lo:- Similar-ish to Lo in RT_YInRange ( -255 <= Lo <= 255). [NOTE, min is -255, as eg 0 - 255 = -255].
Hi:- Similar-ish to Hi in RT_YInRange. ( Lo <= Hi <= 255). Lo and Hi Mode 6 only.
Lo and Hi set the difference range that will be detected.

ClipDiff modes except for modes 2 (MinMaxDifference), 5(Stdev), and 6 (InRange) return frames where 128 represents no difference
[as in Subtract()].

Mode 2 MinMaxDifference, returns frames relative to 0 (no difference), minmaxdifference differences greater than 255 will be
limited to 255.

Mode 5 Stdev, returns frames relative to 0 (no stdev in difference), stdev differences greater than 255 will be limited to 255.
EDIT: Above in blue, Dont think can be greater than 255

Mode 6 (InRange) returns frames where 255 represents 100% population of sampled frames have difference in range lo -> hi.
So eg, ClipDiffs(Last,Ref,delay=0,mode=6,lo= -255,hi= -1) set pixel/channels to 255 where all pixels/channels (from frames sampled so far)
are less than the same pixel/channel reference frame. Where Lo=0 and Hi=0, would show 255 all pixels (from frames sampled so far) that are
exactly the same as refernce frame pixels.
There is no separate Lo/Hi for individual R,G,B channels, perhaps it would be better to move this over to Planar Y only.

The plugin is compiled with debugging enabled, get DebugView to see what it is doing that is taking so much time :(

*/

Ref = ClipBlend(delay=0).Trim(Framecount-1,-1) # Ref is Average of all frames (Requires ClipBlend plugin)
#return Ref
ClipDiffs(Last,Ref,delay=0,mode=4,lo=0,hi=0,FN="ClipDiffs.Data")

T=StackHorizontal(Ref,Last)

Return T

# Show Amplified difference compared with Clipblend(0) ONLY INTENDED for Mode=4,delay=0 : should be same at last frame.
D=Levels(128-32,1.0,128+32,0,255)
B=StackHorizontal(ORG,D)
StackVertical(T,B)



Plugin supplied with VS 2008 build files.'

[jmac698]

Thanks so much! I don't have time to try it right now but I look forward to it. You're one of the most helpful people in this community :)

[StainlessS]

You're one of the most helpful people in ...

What you mean, one of !

A few words of note.
The ClipDiff thing is a little weird to get your head around.
Thresh applies to frames (whereas in eg YPlanemin it applies to Y pixels), so eg Thresh=1.0% for mode=MIN, means ignore up to 1% of frames (as the differences are counted separately for each and every pixel/channel in the frame). I've implemented almost exactly as in RT_Stats but on differences rather than Y levels.
Also, Mode=STDEV, on first frame will always produce gray $808080 frame as there is only a single difference to the corresponding pixel in reference frame (stdev of differences rather than eg Y levels).

Its real slow, perhaps a good idea to have Data file on Sata-3 Solid state if you have one.

Anyways, have a play and point out what mistakes I've made :)

EDIT: Oh yeah, MinMaxDifference is the MinMaxDifference of Differences.

[ajk]

@jmac698

How far do you need the processing to reach? My median filter can extract the temporal minimum, maximum, average or median up to 25 frames at the moment, but that is just an arbitrary value and can be raised if needed.

[StainlessS]

Jmac, can you clarify, is the reference frame a separate frame not existing in the temporal clip ? (that is what I implemented in ClipDiffs).

[StainlessS]

Posts 8 and 9 updated.

Fixed a couple of problems, wish life were not so entertaining.

Made Mode = 3, MinMaxDifference relative 0, not 128. Same for mode=5, Stdev.

Check if I got em right JMac.

EDIT: This bit

BYTE Diff[513];
for (i=0; i<=512; i++) Diff[i] = max(0,min(255,i-129));


From Layer.cpp, Subtract Filter.
Use as in


dstp[x] = Diff[z + 128 + 129];


To limit range, offset by 128, where z could be -255 -> 255.