I enjoyed your discussion very much, I just finished my page yesterday, Please check lucky nugget win cash webpage devoted to lucky nugget win cash. lucky nugget great chance webpage devoted to lucky nugget great chance. !

The grades should be posted as of yesterday (april 4th). Sorry for the delay.

-Neil

What happened to our grades?!

When will we get our grades?

Thanks!

He said that the final will be heavily focused on Just in Time Compilation in Matlab and on the Sprite Operating System.

could someone post any notes they took at today's Q & A session?

When we refer to the two windows, (the second of which is shifted by the previous result in MSLK)

The gradient of the first window is used in producing A,
The difference of the two windows is used in producing B...

Is this correct?

OK, seems every time I make a post I end up confusing things for people. The optical flow results I posted for the corridor are from images 1 and 2 instead of 0 and 1. So if you're wondering why the optical flow looks funny that's why.

-Neil

Brian,

Yeah, you probably don't need to explicitly apply a Gaussian filter, but I was just making sure...

----

Yuzhe,

In my implmentation, if a window goes outside the image, I truncate it so it only contains valid portions of the image (an equivalent technique would be to zero pad the image). I found this helped a lot because otherwise the edges of the image don't get an optical flow estimate at low-res. I recommend doing this if you haven't already, but it's not really stated anywhere that you should do this so consider it optional.

Hi Neil,

For the sphere image, you used lambda2 > .01 & lambda1/lambda2 < 100, windowsize = 15x15, 6 levels. Well, the sphere image is 200*200, then the pyramid should be 200, 100, 50, 25, 12, 6. Then how could we use the 15*15 window on a 6*6 image on the 6th level? I think the window has fixed size at each level, right?

Thank you!

Is explicit gaussian smoothing necessary? If we use imresize(), isn't it supposed to do that automatically?

From the matlab doc:

When the specified output size is smaller than the size of
the input image, and METHOD is 'bilinear' or 'bicubic',
IMRESIZE applies a lowpass filter before interpolation to
reduce aliasing. The default filter size is 11-by-11.

Oh, one more thing I should mention: I'm using Gaussian smoothing at each level of the pyramid..

-Neil

OK, I guess image tags don't work in this discussion board...here's links to the images:

sphere result

corridor result

Ali,

In my implementation I used 5-6 levels (at least for the sphere and corridor images).

-----------------------------

Ramsin,

When the threshold test fails, you should just keep your previous estimate of the optical flow at that pixel. One of the advantages of multi-scale optical flow is that it can work in regions with little texture (without a huge window) because the optical flow is propagated from the low-res estimation.

Also I think you're confused about the thresholds..there should only be two thresholds. You should check that lambda2 > th1 and that lambda1/lambda2 < th2. I found that values of th1 = .01 and th2 = 100 worked fairly well.

----------------------------

Nikhil,

Below are two results from my implementation. For the corridor image I am using lambda2 > .1 & lambda1/lambda2 < 100 with a window size of 9x9 and 7 levels. For the sphere image, I am using lambda2 > .01 & lambda1/lambda2 < 100, windowsize = 15x15, 6 levels.

<img src="http://www.neilalldrin.com/shared/files/of_sphere.jpg" />

<img src="http://www.neilalldrin.com/shared/files/of_corridor.jpg" />

Neil,
How well do we expect the MSLK to work for the sphere and corridor images?
Thanks

Hey Neil,

What do we do in the cases where the thresholds fail? Do we just set u and v to 0 at those places? Your explanation is still really fuzzy to me. Also, how many thresholds are we supposed to have? Your response to Ali made it sound like we need to check if lambda1 is large, lambda2/lambda1 is less than but close to 1, lambda2 is large, and lambda1/lambda2 is small but not smaller than 1. Looking at my lambda values for the corridor images, those 4 thresholds will never hold together, so I'd pretty much create a u,v field of all zeros. Is that what you want? Sorry, this isn't too clear to me.

neil:
how many layers should the gaussian pyrimad have? 3,4,5?

Ryan,

If you are comparing windows and one of the thresholds fail, it means you shouldn't update your flow field at that position because you will get a bad result (A'*A is either ill-conditioned or the conditioning is due to noise).

--------------

Ali,

It sounds like you understand what's going on. The best optical flow estimates will come from regions with lots of texture; in this case lambda1 will be large and lambda2/lambda1 will be close to 1. Basically, the higher lambda2 is the better and the lower lambda1/lambda2 is the better (where lambda1 > lambda2).

--------------

Yuzhe,

Suppose you are computing the optical flow at position (x,y) and the current optical flow estimate is given by u(x,y) and v(x,y). Then you would use windows of (window size = [5 5]),

W0 = I0(y-2:y+2,x-2:x+2);
W1 = I1(y-2+v(x,y):x+2+u(x,y));

So basically, the optical flow at the center pixel is used to shift the entire window the same amount.

-Neil

Hi everyone,

While implementing the algorithm, shall we shift the window pairs in a "pixel by pixel" manner, or should we divide the whole image into blocks and then figure out the velocity for each block? I think according to our assumption, pixels within a certain window are supposed to have the same velocity. So, if we do it in a pixel by pixel fashion, is it a contradiction to our assumption? I got confused here.

Thanks!

Edited by author 03-11-2006 09:30 PM
Neil:
I still dont get why we need setting thresholds for lambda1 and lamda2. Well I do understand when lamda1 and lamda2 are small it means low textured region(i.e. constant region) so there wont be much motion in those regions. I do also understand that when lamda1 is much larger than lamda2 we have an egde and gradients will point most in one direction and could cause matrix A to become singular. And also when both lamda1 and lamda2 are large it refers to a high textured region. So what do the thresholds that you talked about in the assignment mean? and how could they help us
Thanks

Neil,

So, I understand a bit more about the eigenvalues of A'A. However, I still do not understand what to do once I find that, for example, lambda1/lambda2 is too large. I understand that if lambda1/lambda2 is large then this corresponds to the gradient being the same throughout the whole window, but how does knowing this help modify A'A? Thanks

The free program FineView (scroll down on the following page: http://www.snapfiles.com/Freeware/gmm/fwviewers.html) will do it as well.

Ryan,

I used my ACDSee 7 to convert those files into .bmp format. Let me know if you still need them.:)

Just as a note to all, Matlab 6 (Release 12) DOES NOT support .ppm or .pgm images. Can someone save these images as a different file extension? If not, I guess I will have to find a newer version of Matlab...

Hey Ryan,

Here's one way to think about it. What kind of image windows will yield poor optical flow estimate? One is where the spatial gradients of every pixel point in the same direction. lambda1 / lambda2 indicates these type of regions (if this number is high, then the rank of A'*A is close to 1). Another case where the estimation will fail is if the spatial gradient is close to zero at every pixel (ie, rank(A'*A) is close to 0). Because of noise, such a region may be well conditioned (ie, lambda1 / lambda2 close to 1). That's where the check on lambda2 comes in.

Neil,

I vaguely remember the discussion in class about the eigenvalues of A. Could you explain a little more about why they would help us with teh Lucas-Kanade algorithm and how we can use them as thresholds?

OK, I apologize for the confusion that may be caused by my example. I based it on my implementation and thought the example had an unambiguous solution, but it turns out there are multiple paths with the same lowest cost. As Alex figured out the hard way, if given a choice for the "best" path (ie, more than one of left, up, or upleft have the same path cost) I choose upleft over left and left over up.

-Neil

Just thought I'd point this out for anyone who's not getting Neil's newsgroup example exactly, like if you're getting one too many -2's on the end. Here's a fix a few of us have found.

In the pseudocode in the paper, you'll see something like the following:

if(min1 == cmin) Mij = 1
if(min2 == cmin) Mij = 2
if(min3 == cmin) Mij = 3

If you code this directly, then you'll get situations where ties between min1 and min3 go to min3, because Mij gets overwritten. Essentially, ties go to the occlusion. Change these lines to an if-else clause so that ties go to min1 instead. Then you'll get the right answer. Shouldn't change your picture results much, if at all.

Hope this helps someone,
-Alex

OK, it seems lots of people are having problems with the ARROYO images. Therefore, feel free to optionally substitute the TREES images in place of the ARROYO images.

Note: this is optional, so if you've already completed everything using the ARROYO images, don't sweat it.

-Neil

C gets runtimes around a second if you optimize a bit.

Edited by author 03-05-2006 09:37 PM
Yeah, I think the blotchy occlusions are normal. If you work really hard at optimizing your code, rewriting some of the matlab built-in functions, and pre-processing, etc, you can get runtimes well under a minute within matlab. (~15 seconds)

Regarding Part C, when using an occlusion penalty of 0.1 and a window size of 9, I get really splotchy occlusions, 10s of pixels thick. I understand a larger window is supposed to blur my results, but is this normal? Are the occlusions supposed to thicken or disappear?


BTW, for anyone doing the C++ version of this project, it seems WAY faster. Times are under a minute for 9x9 windows.

.01, but I understand some tuning is necessary.

I get about 30 minutes for each graph in Part 3. What occlusion_penalty are you using to get results? The results that I get are filled with occlusions.

What sort of runtimes are people getting for part C? I'm trying my best to optimize my code but SSD runs are still taking around 30 minutes using reasonable disparity limits and a 5x5 window...

I don't know about these "arroyo" images.. I think the difference in lighting conditions really messes things up here in Part B. Maybe that's the point of the assignment, to see how lighting conditions can foul up disparity matching?

Anyway, regarding the equalization: if you crop the arroyo images to remove the black bars from rectification, then do a histeq() equalization, the lighting conditions match a little better between the images (not perfect, but better). This also seems to reduce excess occlusions (using penalty 0.01).

I tried stereo2d on the "trees" and "corridor" images, which have notably better lighting than "arroyo." The results also look better, suggesting that the difference in lighting is a major problem here.

Yeah, the left image is darker... should have noticed that when choosing the stereo pair. Anyway, it might work somewhat better to scale the image intensity of the left image like you say..

Hi there,

I find a problem in Part C, ARROYO. The left image is obviously darker than the right one. If we directly manipulate the intensities of pixels, I think we might get lots of false "matching" and/or false "occlusion". I guess a reasonable way to fix this problem is that we adjust one of them to be more like the other, in terms of "average" brightness. Then, apply our algorithms to the "normalized" image pair.

Is this idea OK?

Thanks~

Ryan: Enforcing the uniqueness constraint is relatively straightforward... if you're storing the best path correctly, this should just boil down to how you do the backtracking. As for part C, you will need to adjust the occlusion penalty (at least for the unnormalized versions). Try some different values and report the value that works best.

Bhavjit: Yes, it does depend on the occlusion penalty. I'm using an occlusion penalty of .01.

Ryan: It could be that you are calculating the Sum of the difference of the 5*5 matrixes and then squaring it. It is the Sum of the squared(.^2) differences.

Does any one know a good occlusion penalty for the other 2 cost fuctions(SAD and Normalized SSD)? should we, for example, choose a smaller penalty for the SAD becease we are not squaring it?

Bhavjit: It is specific; you can get Neil's answer by using an occlusion penalty of 0.01.

Neil, as for the 3rd point you made, isnt the answer very specific to the value of occlusion penalty chosen?

Is anyone getting their disparity maps to look correct when using the 5x5 or 9x9 windows? My part B looks great, but for the Part C I continue to get large areas of occlusions. I have only modified the stereo_2d and stereo_1d algorithm slightly.

Neil,

Regarding the uniqueness constraint for the pixel-by-pixel case, how are we supposed to enforce it?

We need to fill out the path matrix (the M(i,j) matrix in the paper) in order to determine the optimal path, which means some pixels in the left can map to some pixels in the right. So, how do we get around this or am I understanding this correctly?

And for your Note 2 below, did you give that to us just for our information or do we need to implement it into our program?

Regarding HW3, I just wanted to clarify a few things and standardize a couple of ambiguous definitions to make both our lives easier.

1) Suppose (u_left,v) and (u_right,v) are corresponding pixels on scanline v. The disparity can be defined as either u_right-u_left or u_left-u_right. Following the book's definition, please use d = u_right - u_left and use the left image as the reference frame (ie, your disparity map should contain the disparity for each pixel in the left image).

 1a) Note 1: This definition of disparity means a typical stereo pair will have negative disparity for most pixels. When displaying your disparity maps, you may have to negate the pixel values.

 1b) Note 2: The depth at a point given baseline B and disparity d is z = -B/d.

2) Enforce the uniqueness constraint. This means each pixel in the left can map to at most 1 pixel in the right scanline and vice-versa. Occluded pixels are ones that don't have a corresponding pixel in the other scanline.

3) Try your 1D code on the following input:

v1 = [1 0 1 0 1 1 0 1 1 1 0 1 0 1 1 0];
v2 = [1 0 1 1 1 1 1 1 0 1 1 0 1 0 1 0];

You should get

d = [0 0 0 X -1 -1 X -2 -2 -2 -2 -1 -1 -1 0 0];

where X indicates an occlusion.

Arash, have you found a solution? I couldn't find where mine went wrong.

I hope I can get some volcanic crater showing on my screen by 3:30 today ...

I had the same problem. My issue was that there was a small error in my rotation matrix. Test your rotation matrix by multiplying R*R'. If it's not the identity matrix, double check your work. Otherwise, I don't know. The best way I could describe those images is like a volcano crater.

Patrick, I'm having the same issue. Let me know if you found a solution...

Patrick,

You should get similar results to those shown in the paper... if you want send me an email with your results and I'll take a look.

As for equation 4, it's just a restatement of equation 3 for the diffuse channels U and V. D is the diffuse surface color in rgb space and r_2 and r_3 yield the U and V components of D. dot(n,l) = cos(theta) is just the standard Lambertian foreshortening term. Thus we get I_u and I_v, the pixel intensities in U and V color space.

is it too late to start on the project? !!!!!

When I integrate the notmals to get the surface in part3 (for sphere, using horn integration), the picture comes out as a flat surface in the middle of the sphere, and a lot of spikes in the boundry. I wonder what might be going wrong. The RGB->SUV conversion seems to be able to separate the specular component from the diffuse component as far as I know.

In the paper for equation 4, I am wondering if all the terms after r2' just means the color value of a pixel.

A.M.

You can expect pretty crummy results using a naive integration technique on real data (which the sphere and pear are). Averaging across multiple integration paths might improve results -- if you want, try it out. As for the "mask", this is just a binary image the same size as your input images that specifies whether a given pixel is part of the surface you're reconstructing or not. This is needed because some of the pixels are not part of the surface. One naive way to get the mask is to threshold pixel intensities so that everything below some intensity value will be considered "background".

-Neil

Billy,
In quiver3(Z,U,V,W), Z is the 2-D surface matrix you computed with the integration technique. U, V, and W are the normal vectors in each direction (x, y, and z), which you computed when you solved for the "b" vector at each pixel. The dimensions of Z, U, V, and W are all 100x100, i.e. the size of the image.
So quiver3() essentially plots the normal vector at each pixel point; Z just determines the surface along which to place the arrows.
FYI: your plot may looked cramped if you don't subsample. You might try indexing each matrix with (1:5:end, 1:5:end) to space out the arrows.

I don't understand what you specify in the quiver function

T:

It turned out my problem had nothing to do with rounding, i was simply calculating the wrong angle. In my case, the matrix ended up acting like a projection instead of a rotation.

Edited by author 02-11-2006 06:25 PM
Neil or anyone else who can help me out:
when doing the integration the naive way to recover the surfaces for the real images I get somthing that is not even close to a sphere or a pear. Is this expected?Also, If I do the intergration along different ways and average them will that make it better?
In the Horn integration approach( the code is provided) what is "mask" and how do we obtain it?

Thanks

So, I'm getting the same problem as someone who posted earlier. My rotations for RGB->SUV don't come out exact. They're all off by less than .1 but by a non-trivial amount. The trig for figuring out the angles is right and the rotation matrices around the axies are right. Is there some magic declaration to make to have it keep matrices in symbolic mode or something that I should be doing?

In fact it is :).

-Neil

is the deadline REALLY postponed!!!

It has to do with the direction of your nx and ny normals and how Matlab plots images. Just experiment with flipping around the two normals in quiver3.

When I plot the normal field with 'quiver3' for each result, I found that the arrows were not pointing to the normal direction, but the tangent of the surface.

Does anyone have the same problem? Where did I get wrong?

Thanks.

Who should we email the report to?
The TA or the professor?

Oh, one more thing, I got a followup email from Satya and yes, the reason the images have small negative values is because he subtracted out the ambient lighting (the negative values are from noise). So just use Inoneg = max(I,0) to remove any negative components.

-Neil

The best way is to convert to double AND scale the numbers so they range from 0 to 1. You can do this by issuing the command,

Idouble = double(I)/intmax('uint8');

or equivalently,

Idouble = double(I)/255;

-Neil

Your variables need to be of type double, rather than uint8. So, you need to typecast either the vector E or matrix S to a double by the following:

double (E);
double (S);

You can see in Matlab by looking at the workspace what the type of the variables are.

Hope this helps.

I create a function to computer a matrix. Inside the function, E is a vector and S is a matrix. the line b = E * S cause matlab to complain
"Function '*' not defined for variables of class 'uint8'"
What does it mean? Could anyone help? Thanks

I'm almost positive the reason you're seeing small negative values in the images is that the ambient lighting component was subtracted out (it's practically difficult to generate lighting from a single point). The negative values are probably the result of noise in this process. I've emailed Satya to verify this.

The negative values shouldn't affect your results, but it might be prudent to use im_noneg = max(im,0) to get rid of negative components.

-Neil

The reason the images are dark is that they are high dynamic range images and the images are scaled to a small value. If you want to see the image on your screen you need to scale it. Try imshow(im1*2^6) and you'll see the image.

For those interested, HDR imaging involves merging multiple images taken at different exposures. The prevents saturation and undersaturation (which occur because the image sensor has a limited range of light intensities it can measure). HDR images are typically stored as floating point numbers instead of the standard 8-bit images most cameras take.

Hope that helps!
-Neil

Yeah, the images have negative values, but they display fine with rgb2gray(). I'm wary of scaling the image values because the rotation in the algorithm depends on the color of the light source. Is it ok to scale the images without changing the source color?

Edited by author 02-04-2006 09:12 PM
I got the same problem---black images.

But I just scaled these values by some factor, then you can get some dark red spheres somewhat like the spheres appearing at the end of the paper.

Well, could someone have more reasonable solution?

Is anyone able to view the RGB images in any of the .mat files?

When I look at the images, most of them have very small pixel values (e.g. on the order of 10^-3) and thus, I get a black image when I show the image using imshow().

It should be possible to get it exact (or at least closer than that).

-Neil

In RGB->SUV conversion, the rotation of the light source coordinates (R*S) should map onto some axis (eg red = [1 0 0]'), correct?

I think i've implemented this correctly, but the rotation is slightly off the axis ([0.9898 0 -0.1425]'). Can I blame this on numerical underflow and call it close enough? Or is it more likely a bug?

part 2 and 3 of the photometric stereo assignment have not been covered in class yet and the HW is due on thursday! Is this somekind of joke?

Well, I know it looks like a lot, but it's hopefully not as hard as it looks. Photometric stereo at it's core is a relatively simple idea as is the specularity removal technique. Just start trying to work through each part of the assignment in turn and I think you'll be ok.

-Neil

you've got to be joking... where do we even start on the project?

For HW1 problem 4, there is a little bit of confusion regarding the position of the image plane. Basically, the image plane in this problem is a _virtual_ image plane because it is in front of the camera. Virtual image planes are equivalent to real image planes except the resulting images are flipped.

-Neil

why is this class so hard?

Welcome to the CSE 252A discussion board. Feel free to post questions or comments related to the class.

Enjoy :).