# Enjoy

Member

37

355 Neutral

• Rank
Member

• Interests
Programming

## Social

• Github
https://github.com/suzdalnitski
1. ## A simpler introduction to Non-Local means denoising filter

Non-Local means is a great algorithm to denoise an image. It's quite straightforward, however it may be difficult to understand for people not familiar with image processing. I did my best to explain the algorithm without complex math equations. I would highly recommend first watching the following video for an introduction to the algorithm by Prof. Fred Hamprecht: The algorithm is based on the idea that every image has similar regions. I.e in the following image regions q1 and q2 are similar to region p. If we take similar regions and average them, then we'll get a denoised version of that region. In the image above regions q1 and q2 are very similar to region p, however region q3 is very different, therefore regions q1 and q2 will have a higher weight, and region q3 will have a lower weight. The denoised version of p can be described by the following simple formula: [indent=1]p = V(q1)*W(q1) + V(q2)*W(q2) + V(q3)*W(q3)... V() is the grayscale value of a central pixel of a certain region, and W() is the weight of the whole region (higher weight = more similar to the original region). The weight W() will depend on the similarity between corresponding pixels in the original and searched region as well as the distance between the original and search regions. Example: Let's say that e is the original pixel that we want to denoise, a-i is the region around that pixel: [indent=1]1 2 3[color=#fff0f5]_____[/color]a b c [indent=1]4 5 6[color=#fff0f5]_____[/color]d e f [indent=1]7 8 9[color=#fff0f5]_____[/color]g h i 1-9 is the searched region. The weight of the searched region 1-9 will depend on the distance between pixels 5 and e (central pixels) and on the sum of differences in the grayscale value of 1 and a, 2 and b, and so on. The non-local means algorithm uses a region of 7x7 pixels (around the pixel to be denoised), and looks for similar patterns in a 21x21 pixels neighbourhood. Search region weight W = ( ( F(1,a) + F(2,b) + ...) / Fn ) * G(1,a); G is gaussian function and F is the non-local means function: [indent=1]G = ( 1 / (2*PI*sigma^2) ) * exp(- (x^2 + y^2) / (2*sigma^2) ) [indent=1]F = exp( -( V(p1)-V(p2) )^2/h^2 ) In the above equations sigma is the gaussian strength. In our algorithm it should be equal to 10 (radius of the 21x21 search neighbourhood). x and y is the distance vector of a searched region's center (pixel 5) from the pixel's to be denoised (pixel e) center. h is the denoising strength. Fn is the normalizing factor for the searched region: Fn = F1 + F2 + F3... The following pseudo-code will calculate the weight of a certain window (c and d is the offset in the 21x21 search neighbourhood):for (int a = c; a < c + 7; a++) { for (int b = d; b < d + 7; b++) { int2 srcCoord = (int2)(coord.x + a - c - 3, coord.y + b - d - 3); int2 nhCoord = (int2)(coord.x + a - 10, coord.y + b - 10); //Value of a corresponding pixel surrounding the pixel to be denoised (a-i) float pValue = srcImg.GetPixel(srcCoord); //Value of a corresponding pixel in the neighbourhood (1-9) float nValue = srcImg.GetPixel(nhCoord); float deltaValue = pValue - nValue; float curWeight = exp( - (deltaValue * deltaValue)/ (h*h) ); windowWeight += curWeight; }} Now windowWeight will have the 0..1 weight of the current searched region. Then we simply have to multiply the value of the searched region's (1-9) central pixel by window weight and add it to the totalValue:windowWeight *= gaussianWeight; totalWeight += windowWeight;totalValue += windowCenterPixel * windowWeight; totalWeight is the normalization factor that will be used later. totalValue will be used for calculating a denoised value of the pixel. After going over all possible 7x7 search regions in the 21x21 neighbourhood we normalize the totalValue: totalValue /= totalWeight; totalValue is the denoised pixel. I hope this introduction was helpful. Here's an OpenCL kernel, that implements the ideas described in this article://This is free and unencumbered software released into the public domain.////Anyone is free to copy, modify, publish, use, compile, sell, or//distribute this software, either in source code form or as a compiled//binary, for any purpose, commercial or non-commercial, and by any//means.////In jurisdictions that recognize copyright laws, the author or authors//of this software dedicate any and all copyright interest in the//software to the public domain. We make this dedication for the benefit//of the public at large and to the detriment of our heirs and//successors. We intend this dedication to be an overt act of//relinquishment in perpetuity of all present and future rights to this//software under copyright law.////THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,//EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF//MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.//IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR//OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,//ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR//OTHER DEALINGS IN THE SOFTWARE.////For more information, please refer to int2 GetWrappedCoord(const int2 _input, const int2 _imageDimensions) { int2 c = _input; while (c.x < 0) { c.x += _imageDimensions.x; } while (c.x >= _imageDimensions.x) { c.x -= _imageDimensions.x; } while (c.y < 0) { c.y += _imageDimensions.y; } while (c.y >= _imageDimensions.y) { c.y -= _imageDimensions.y; } return c;}__kernel void denoiseNLM(__read_only image2d_t srcImg, __write_only image2d_t dstImg, const uint xOffset, const uint yOffset, const float filteringStrength){ const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | //Natural coordinates CLK_ADDRESS_CLAMP_TO_EDGE | //Clamp to zeros CLK_FILTER_NEAREST; int2 imageDim = get_image_dim( srcImg ); int2 coord = (int2)(get_global_id(0) + xOffset, get_global_id(1) + yOffset); float totalValue = 0.0f; float totalWeight = 0.0f; int searchArea = 21; int halfSearchArea = (searchArea - 1) / 2; int searchAreaMaxStart = searchArea - 7; float gaussSigmaSQR = convert_float(halfSearchArea); gaussSigmaSQR *= gaussSigmaSQR; for (int c = 0; c < searchAreaMaxStart; c++) { for (int d = 0; d < searchAreaMaxStart; d++) { float windowWeight = 0.0f; float windowValue = 0.0f; for (int a = c; a < c + 7; a++) { for (int b = d; b < d + 7; b++) { int2 pWrapped = GetWrappedCoord( (int2)(coord.x + a - c - 3, coord.y + b - d - 3), imageDim ); int2 nhWrapped = GetWrappedCoord( (int2)(coord.x + a - halfSearchArea, coord.y + b - halfSearchArea), imageDim ); //Point float pValue = native_divide(convert_float( read_imageui(srcImg, smp, pWrapped).z ), 255.0f); //Neighbourhood float nValue = native_divide(convert_float( read_imageui(srcImg, smp, nhWrapped).z ), 255.0f); float deltaValue = pValue - nValue; float curWeight = native_exp( - native_divide( (deltaValue * deltaValue), (filteringStrength * filteringStrength) ) ); windowWeight += curWeight; } } //int2 windowCenterCoord = GetWrappedCoord( (int2)(coord.x + c - 7, coord.y + d - 7), imageDim ); int2 windowCenterCoord = GetWrappedCoord( (int2)(coord.x + c - halfSearchArea + 3, coord.y + d - halfSearchArea + 3), imageDim ); float windowCenterPixel = native_divide(convert_float( read_imageui(srcImg, smp, windowCenterCoord).z ), 255.0f); //Add distance factor to initial patch float2 diffVector = convert_float(windowCenterCoord - coord); float gaussianWeight = native_exp( - native_divide( dot(diffVector, diffVector), (2.0f * gaussSigmaSQR) ) ); gaussianWeight = native_divide(gaussianWeight, 2.0f * M_PI * gaussSigmaSQR); windowWeight *= gaussianWeight; totalWeight += windowWeight; totalValue += windowCenterPixel * windowWeight; } } totalValue = native_divide(totalValue, totalWeight); uint4 bgra; bgra.x = bgra.y = bgra.z = (uint) (totalValue * 255.0f); bgra.w = 255; write_imageui(dstImg, coord, bgra);}
2. ## GPGPU image processing basics using OpenCL.NET

OpenCL is a cross-platform framework used mostly for GPGPU (General-purpose computing on graphics processing units). There are plenty of tutorials available on image processing with OpenCL using C/C++, however there's not much information that would cover OpenCL image processing with .NET. I won't go into details about OpenCL kernels/queues/etc. (there's plenty of information available on the internet), however I'll provide you with a bare minimum code required to load an image from disk, process it with OpenCL on the GPU and save it back to a file. Before we get started, make sure that you download the source code of OpenCL.NET from http://openclnet.codeplex.com/ and add it to your project. We'll use a simple OpenCL kernel that converts an input image into a grayscale image. The kernel should be saved to a separate file. Kernel source code: __kernel void imagingTest(__read_only image2d_t srcImg, __write_only image2d_t dstImg) { const sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | //Natural coordinates CLK_ADDRESS_CLAMP_TO_EDGE | //Clamp to zeros CLK_FILTER_LINEAR; int2 coord = (int2)(get_global_id(0), get_global_id(1)); uint4 bgra = read_imageui(srcImg, smp, coord); //The byte order is BGRA float4 bgrafloat = convert_float4(bgra) / 255.0f; //Convert to normalized [0..1] float //Convert RGB to luminance (make the image grayscale). float luminance = sqrt(0.241f * bgrafloat.z * bgrafloat.z + 0.691f * bgrafloat.y * bgrafloat.y + 0.068f * bgrafloat.x * bgrafloat.x); bgra.x = bgra.y = bgra.z = (uint) (luminance * 255.0f); bgra.w = 255; write_imageui(dstImg, coord, bgra); } Namespaces Used using System; using System.Collections; using System.Collections.Generic; using System.Drawing; using System.Drawing.Imaging; using System.IO; using System.Runtime.InteropServices; using OpenCL.Net; Error handling Since OpenCL.NET is a wrapper for C API, we'll have to do all the error checking on our own. I'm using the following two methods: private void CheckErr(Cl.ErrorCode err, string name) { if (err != Cl.ErrorCode.Success) { Console.WriteLine("ERROR: " + name + " (" + err.ToString() + ")"); } } private void ContextNotify(string errInfo, byte[] data, IntPtr cb, IntPtr userData) { Console.WriteLine("OpenCL Notification: " + errInfo); } Setting Up The following two variables should be declared in the class itself and will be shared across all of the methods: private Cl.Context _context; private Cl.Device _device; And this is the method that sets up OpenCL: private void Setup () { Cl.ErrorCode error; Cl.Platform[] platforms = Cl.GetPlatformIDs (out error); List devicesList = new List (); CheckErr (error, "Cl.GetPlatformIDs"); foreach (Cl.Platform platform in platforms) { string platformName = Cl.GetPlatformInfo (platform, Cl.PlatformInfo.Name, out error).ToString (); Console.WriteLine ("Platform: " + platformName); CheckErr (error, "Cl.GetPlatformInfo"); //We will be looking only for GPU devices foreach (Cl.Device device in Cl.GetDeviceIDs(platform, Cl.DeviceType.Gpu, out error)) { CheckErr (error, "Cl.GetDeviceIDs"); Console.WriteLine ("Device: " + device.ToString ()); devicesList.Add (device); } } if (devicesList.Count
3. ## Flowmaps distortion

I was able to solve the issue. In case somebody else encounters the problem, you just have to set the cycle to a much lower value. The optimal value for cycle is somewhere around: [CODE]cycle = 15.0f;[/CODE]
4. ## Flowmaps distortion

[quote name='Postie' timestamp='1355402409' post='5010175'] My guess is that setting the offset back to 0 when it exceeds the cycle length is the culprit. At faster flow rates the offset will overshoot by more than with a slower rate, so there will be a more obvious discontinuity. Try subtracting cycle from the offset instead of setting it to 0. [/quote] No luck so far... Any other ideas? Using higher resolution flowmap didn't solve the problem.
5. ## Flowmaps distortion

I've implemented flowmaps based on the following article: [url="http://graphicsrunner.blogspot.ru/2010/08/water-using-flow-maps.html"]http://graphicsrunner.blogspot.ru/2010/08/water-using-flow-maps.html[/url] For smaller values of flowSpeed everything is fine, but for larger values the water surface is getting extremely distorted. Is there any solution to this problem? Here's how do the FlowMapOffset calculation: [CODE] flowMapOffset0 += flowSpeed * Time.deltaTime; flowMapOffset1 += flowSpeed * Time.deltaTime; if ( flowMapOffset0 >= cycle ) flowMapOffset0 = .0f; [/CODE] Small flowSpeed: [img]https://dl.dropbox.com/u/16265778/Water_undistorted.jpg[/img] Large flowSpeed: [img]https://dl.dropbox.com/u/16265778/Water_distorted.jpg[/img]
6. ## Tangent space rotation

You're right, reflection vector should be in world space, and not in tangent space.
7. ## Tangent space rotation

I'm trying to rotate a view direction vector in my shader into tangent space. I need this for cubemap reflections to work on an arbitrary normal mapped surface. However for some strange reason the resulting reflection is rotated by 90 degrees in the YZ plane. It used to work correctly until I introduced the tangent space rotation. I suppose that the rotated viewDir vector is rotated incorrectly by 90 degrees. What's the reason and how can I fix it? Here's a snippet from my CG shader: [CODE] v2f vert (appdata_tan v) { v2f o; o.pos = mul (UNITY_MATRIX_MVP, v.vertex); float3 binormal = cross( v.normal, v.tangent.xyz ) * v.tangent.w; float3x3 rotation = float3x3( v.tangent.xyz, binormal, v.normal ); o.viewDir = mul (rotation, ObjSpaceViewDir(v.vertex)); } fixed4 frag (v2f i) : COLOR { half3 refl = reflect( normalize(i.viewDir), half3(0.0, 1.0, 0.0) ); refl.x = -refl.x; texcol.rgb = texCUBE( _Cube , refl ).rgb; return texcol; } [/CODE]
8. ## Velocity of a rolling ball

Hello, I need help with calculating angular and CM (center of mass) velocity of a rolling ball. The ball is rolling on an inclined surface. What I have: - CM velocity (velocity of a ball itself) - Angular velocity ( w ) - Coefficient of friction ( coeff ) - Normal of a plane What I need to get: - New velocity of center of mass (ball itself) - New angular velocity I have already implemented collision detection and collision response (only calculating new velocity of a ball after collision, based on restitution). But I can't figure out how to make the ball roll properly on the surface. I want to take into account sliding of a ball (if tangent velocity + ball velocity is not equal to 0 then a ball slides on the surface) Thank you in advance.
9. ## Collision detection between spheres

Hello, I want to find when a collision between a static and a moving ball occurs, but the algorithm I came up with, sometimes doesn't detect a collision and the moving ball goes through the static one. The moving ball is affected by gravity and the static one is not. Here's my collision detection code: GLfloat whenSpheresCollide(const sphere2d &firstSphere, const sphere2d &secondSphere) { Vector2f relativePosition = subtractVectors(firstSphere.vPosition, secondSphere.vPosition); Vector2f relativeVelocity = subtractVectors(firstSphere.vVelocity, secondSphere.vVelocity); GLfloat radiusSum = firstSphere.radius + secondSphere.radius; //We'll find the time when objects collide if a collision takes place //r(t) = P[0] + t * V[0] // //d^2(t) = P[0]^2 + 2 * t * P[0] * V[0] + t^2 * V[0]^2 // //d^2(t) = V[0]^2 * t^2 + 2t( P[0] . V[0] ) + P[0]^2 // //d(t) = R // //d(t)^2 = R^2 // //V[0]^2 * t^2 + 2t( P[0] . V[0] ) + P[0]^2 - R^2 = 0 // //delta = ( P[0] . V[0] )^2 - V[0]^2 * (P[0]^2 - R^2) // // We are interested in the lowest t: // //t = ( -( P[0] . V[0] ) - sqrt(delta) ) / V[0]^2 // GLfloat equationDelta = squaref( dotProduct(relativePosition, relativeVelocity) ) - squarev( relativeVelocity ) * ( squarev( relativePosition ) - squaref(radiusSum) ); if (equationDelta >= 0.0f) { GLfloat collisionTime = ( - dotProduct(relativePosition, relativeVelocity) - sqrtf(equationDelta) ) / squarev(relativeVelocity); if (collisionTime >= 0.0f && collisionTime <= 1.0f / GAME_FPS) { return collisionTime; } } return -1.0f; } And here is the updating function that calls collision detection: void GamePhysicsManager::updateBallPhysics() { // //Update velocity vVelocity->y -= constG / GAME_FPS; //v = a * t = g * 1 sec / (updates per second) shouldApplyForcesToBall = TRUE; vPosition->x += vVelocity->x / GAME_FPS; vPosition->y += vVelocity->y / GAME_FPS; if ( distanceBetweenVectors( *pBall->getPositionVector(), *pBasket->getPositionVector() ) <= pBasket->getRadius() + vectorLength(*vVelocity) / GAME_FPS ) { //Ball sphere sphere2d ballSphere; ballSphere.radius = pBall->getRadius(); ballSphere.mass = 1.0f; ballSphere.vPosition = *( pBall->getPositionVector() ); ballSphere.vVelocity = *( pBall->getVelocityVector() ); sphere2d ringSphereRight; ringSphereRight.radius = 0.05f; ringSphereRight.mass = -1.0f; ringSphereRight.vPosition = *( pBasket->getPositionVector() ); //ringSphereRight.vPosition.x += pBasket->getRadius(); ringSphereRight.vPosition.x += (pBasket->getRadius() - ringSphereRight.radius); ringSphereRight.vVelocity = zeroVector(); GLfloat collisionTime = whenSpheresCollide(ballSphere, ringSphereRight); if ( collisionTime >= 0.0f ) { DebugLog("collision"); respondToCollision(&ballSphere, &ringSphereRight, collisionTime, pBall->getRestitution() * 0.75f ); } // //Implement selection of the results that are first to collide collision vVelocity->x = ballSphere.vVelocity.x; vVelocity->y = ballSphere.vVelocity.y; vPosition->x = ballSphere.vPosition.x; vPosition->y = ballSphere.vPosition.y; } Problem seems to be solved when I change FPS from 30 to 10. How does FPS affect my collision detection? Why isn't the collision being detected in 100% of cases? It's being detected only in 70% of cases. Thanks.
10. ## OpenGL Index buffer in OpenGL

Hello, How is an index buffer in OpenGL organized? I understand, how vertex arrays are passed to OpenGL, but I don't understand how to explain OpenGL, which vertex corresponds to which triangle. For example, I need to pass a square. How could this be done, using index buffers? Thank you in advance, Ilya.
11. ## Joining points into ITN

Hello, Is there a fast way to join points into Irregular Triangulated Network? I need an algorithm, that will enable me to do it in real time. Many thanks, Ilya.
12. ## bottom-up terrain

Hello, I am interested in bottom-up terrain LOD algorithms. Where can I read something on this topic? I found only up-down algorithms, where patches are being subdivided into triangles, basing on a needed LOD. I need something reverse - points are being joined into triangles. I hope everything is clear.
13. ## Cloudy fog

Hello, I need to create an effect of volumetric fog. And I want it to be "cloudy". Something like on this image: http://www.gebauz.com/uploads/images/screenshots/pentag/pentag05.png How could this be done?
14. ## Orientation in code

BeauMN, thank you, I will consider your advices and review my architecture. Also thank you for a link. I will try it. I am using c++, working under XCode in MAC OS X.
15. ## Orientation in code

Hello, I am working on a project and it is getting big. For now I have around 5000 lines of code and 20 source files. It is getting really hard for me to orientate in my code - to add a new feachure or just to update something. Really hard to find what I need in that huge classes. I want to know if there is some software, that could help me or an article. Something that could help me. Thank you in advance, Ilya.