Wrong results with pmulhw

[bill_baroud]

Ok, it this not a direct avisynth problem, but it will probably make it into a filter, if i've time. I'm posting here, because you guys are the only mmx wizards i know ;)

I'm programming a "low-pass" filter, using a two-pass algorithm and with a variable kernel size from 2 to 64. For the division by any integer in this range, i'm using a tricks found in the AMD optimization guide, multiplying by a big integer and shifting the high part.

The problem i've now, it's that pmulhw doesn't give me the correct result.
exemple :
my number is 0x03CD
i want to divide by 5, so i multiply by 0xCCCC and right shit by 2

0x03CD * 0xCCCC = 0x030A 635C, so pmulhw should give me 0x030A but it give me 0xFF3D.
What i'm doing wrong, and what i've missed ?

thanks.


my code for the horizontal pass looks like this :

pxor mm6, mm6
pxor mm7, mm7
movd mm6, [shift]
movq mm5, [mmxdividebyN]
movq mm2, [mmxsignedconv]
[...]
mov edi, [tab_tmp]
mov esi, [lpPixSRC]
movq mm0, [esi]
mov ecx, [nb_loops]
psubb mm0, mm2 // unsigned <> signed
mov eax, [shift]
movq mm3, mm0
punpcklbw mm0, mm7
punpckhbw mm3, mm7
loop_noyau:
movq mm1, [esi+ecx]
psubb mm1, mm2 // unsigned <> signed
movq mm4, mm1
punpcklbw mm1, mm7
punpckhbw mm4, mm7
paddw mm0, mm1
paddw mm3, mm4
sub ecx, 1
cmp ecx, 0
jnz short loop_noyau
pmulhw mm0, mm5 //
pmulhw mm3, mm5 // mult + shift
psrlw mm0, mm6 // division by N
psrlw mm3, mm6 //
packuswb mm0, mm3
paddb mm0, mm2 // signed <> unsigned
movq [edi], mm0


ok, i'm screwed, i just understood something and i'm utterly wrong...
multiplying by 0xCCCC is not like multiplying by 0xCCCC0000 like i was thinking... shame one me.
(infact it should be 0xCCCCCCCD, but i needed a 16-bits variable)
So 0xFF3D is just because of wrap-around...i'm totally dumb.

Anybody with another idea for the division ?

[mg262]

Grr... I dictated out a long reply to this and it got trashed. Here goes again -- apologies if it's terse.


I'm programming a "low-pass" filter,
I want it! I want it! Any chance of high pass as well?

How about pmaddwd? That takes in 16-bit input and gives back 32-bit output, so you don't have to worry about rounding. That also means you don't have to worry about the sign/unsigned mess as much.

One other thing occurred when seeing your code... I think sub sets the flag registers so you don't need the subsequent cmp. (I think this is not true of dec). You can also put (i.e. move downwards) instructions which don't affect the flag registers (which I think includes all the MMX instructions) between subtraction and conditional jump, to avoid the delay between them that would normally happen due to dependency.

Incidentally, pmaddwd is particularly suited for kernels. You would need to round the kernel size up to the nearest multiple of two (multiple of four for sse), but I don't see that that is particularly a problem....

And the most important thing was that for any questions of this kind I would ask here:
http://videoprocessing.11.forumer.com/index.php

In particular, I really wouldn't be happy telling you to rely on anything I've said until you had checked it there (feel free to quote it).

I have assemblified kernel code to calculate a weighted average of various frames, which uses pmaddwd, and you're welcome to it if it would be any use to look at.

[bill_baroud]

I want it! I want it! Any chance of high pass as well?
Well, let me finish it first, and i'm quite screwed right now :p
And for the high-pass, what do you mean ? a laplacian ? a gradient ? any other kernel ? (sobel etc ?)

How about pmaddwd?
I don't see how it can help me there, if you can explain... In fact i need 32-bits input for the multiplication, so i think i'm going to use the SSE2 registers to process the 8 signed word i've after my mmx accumulation.


And i've still have to sort out the borders case, it's giving me headache with the horizontal borders and mmx (depending on kernel size, i'm not anymore mod8 etc ... argl)

[kassandro]

The problem i've now, it's that pmulhw doesn't give me the correct result.
exemple :
my number is 0x03CD
i want to divide by 5, so i multiply by 0xCCCC and right shit by 2

0x03CD * 0xCCCC = 0x030A 635C, so pmulhw should give me 0x030A but it give me 0xFF3D.
What i'm doing wrong, and what i've missed ?

........

Anybody with another idea for the division ?

You are on right track, but it's much easier as you think. pmulhw is for division of signed numbers and pmulhuw is for division of unsigned numbers. It automatically makes a right shift by 16 binary digits. Thus MMX division of 16 bit entries in mm0 by 5 is simply

__asm pmulhw mm0, mm1

where the register mm1 contains the values

(1u<<16)/5, (1u<<16)/5, (1u<<16)/5, (1u<<16)/5

(these are 16 bit entities). Actually it is more precise to put into the mm1 register the values:

((1u<<16) + 2)/5, ((1u<<16) + 2)/5, ((1u<<16) + 2)/5, ((1u<<16) + 2)/5

Alternatively you can use a memory operand (should be aligned), but it's faster to have it in a register. For more information or similar questions go to the forum mentioned by mg262.

[Didée]

Sorry to jump in with primitive scripting, but

I'm programming a "low-pass" filter
I want it! I want it! Any chance of high pass as well?
IIRC I've already seen one (http://forum.doom9.org/showthread.php?p=677115&highlight=lowpass#post677115). ;)

Only by now I use

before = last

# make a temporalsoften /w big thresholds, limit result to +-3 ...
LIM = "x 3 + y < x 3 + x 3 - y > x 3 - x 49 * y 51 * + 100 / ? ?"
ts = before.temporalsoften(2,23,23,23,2)
yv12lutxy(before,ts,LIM,LIM,LIM,U=3,V=3)

after = last

radius = 2.0 # Try 1.25 to 3.0

# ... and apply a lowpass to "the effect" to avoid smearing
diff = yv12lutxy(before,after,"x y - 128 +","x y - 128 +","x y - 128 +",U=3,V=3)
diff2= diff.RemoveGrain(4).bicubicresize( m4(diff.width/radius), m4(diff.height/radius) )
\ .RemoveGrain(2).bicubicresize( diff.width, diff.height, 1,0 )
diff = yv12lutxy(diff,diff2,"x y 128 - -","x y 128 - -","x y 128 - -",U=3,V=3)

lowpassed = yv12lutxy( before, diff, "x y 128 - -","x y 128 - -","x y 128 - -",U=3,V=3)

lowpassed
instead.

Strictly taken, it's not fully correct (would require two different [0-255] operations instead of that single 2*[0-128] operation [difference centered @ 128], which I didn't for speed reasons), but works rather nice.

Of course a speedy plugin would be much better :) - wanted only to mention that trick with RemoveGrain(), which comes handy for the lowpass most times.

[bill_baroud]

@kassandro : i'm not keeping up with you there :)

well, in fact the only thing i doesn't understand is, what's the "1u" in the mm1 formula ?

..ok, looked in the debugger, it mean "1" ?
hmm, your method looks a lot more simpler than mine (a big array with multiplier and number of shift). Can you explain a little more how that work ? I don't really like using something without understanding it (like i'm doing right now :D)

and i'm registered on your forum.... beware, i'll come nagging you every 5 minutes with a question ;)


btw, "pmulhuw" come with SSE, and i would like to stick with MMX atm, keeping SSE ehancement for the next version. (even if the target processor is a pentium 4)

[Manao]

pmulhuw comes with integer SSE, which is widely supported, and almost mandatory for encoding, since motion estimation relies on sad which is computed twice as fast in isse as in mmx.

[mg262]

Either I'm not on the same wavelength as the rest of you or (very possible) there is a some important fact about lowpass filters that I don't know about. I'm talking about this kind of thing:

http://en.wikipedia.org/wiki/Lowpass

Now because a multiplication in the frequency domain happens to correspond to a convolution in the spatial domain, for a given low pass filter you will be able to find a kernel, convolution with which is equivalent to applying the low pass filter. In practice the kernel will be infinite and you would have to truncate it, but it would still be quite large.

Is this the kind of object we are talking about? Sorry for the confusion...