Spherical harmonics, meshes and some WebGL and JavaScript : toxiclibs.js, Processing.js, Three.js, ThreeNodes.js

1--------Spherical Harmonics Mesh Builder

http://www.chromeexperiments.com/detail/spherical-harmonics-mesh-builder/?f=

Spherical Harmonics Mesh Builder is an example of Paul Bourke's "spherical harmonics" function for creating organic meshes. This example shows the use of Toxiclibs.js for creating meshes inside Three.js

(Technology: JavaScript, WebGL, toxiclibs.js, three.js).

Spherical harmonics are the angular portion of a set of solutions to Laplace's equation. Spherical harmonics are important in many theoretical and practical applications, particularly in the computation of atomic orbital electron configurations, representation of gravitational fields, geoids, and the magnetic fields of planetary bodies and stars, and characterization of the cosmic microwave background radiation. In 3D computer graphics, spherical harmonics play a special role in a wide variety of topics including indirect lighting (ambient occlusion, global illumination, precomputed radiance transfer, etc.) and recognition of 3D shapes.

The 19th century development of Fourier series made possible the solution of a wide variety of physical problems in rectangular domains, such as the solution of the heat equation and wave equation. This could be achieved by expansion of functions in series of trigonometric functions. Whereas the trigonometric functions in a Fourier series represent the fundamental modes of vibration in a string, the spherical harmonics represent the fundamental modes of vibration of a sphere in much the same way.

is called a spherical harmonic function of degree ℓ and order m.

For a given value of ℓ, there are 2ℓ+1 independent solutions of this form, one for each integer m with −ℓ ≤ m ≤ ℓ. 

A quite good article http://en.wikipedia.org/wiki/Spherical_harmonics

and in french: http://fr.wikipedia.org/wiki/Harmonique_sph%C3%A9rique

http://mathworld.wolfram.com/SphericalHarmonic.html

http://demonstrations.wolfram.com/NodalDomainsOfSphericalHarmonics/

2------------http://haptic-data.com/

http://haptic-data.com/sketches/simplex_trails.html

http://haptic-data.com/toxiclibsjs/

Toxiclibs.js is an open-source computational design library ported to javascript by Kyle Phillips originally written by Karsten Schmidt for Java and Processing. Examples of the original library can be found at http://toxiclibs.org

Toxiclibs.js works great with Canvas, with SVG or any ordinary DOM element. with Processing.js, Three.js, or Raphael.js for SVG...

great example with additive waves

http://haptic-data.com/toxiclibsjs/examples/AdditiveWaves_pjs-webgl.html

3--------three.js

https://github.com/mrdoob/three.js

JavaScript 3D library

The aim of the project is to create a lightweight 3D library with a very low level of complexity — in other words, for dummies. The library provides ,

--ThreeNodes.js

http://idflood.github.com/ThreeNodes.js/public/index.html

Threenodes.js is an attempt to create a dataflow environment in javascript and html5. It's a way to see what is possible to do in the browser. I had a little list of interesting javascript libraries to try but never had the chance to play with them.

The project is still in an experimental state but you can already see a live demo there: http://idflood.github.com/ThreeNodes.js/

coffeescript, compass and haml
These are languages that compile to javascript, css and html. Their syntax are much shorter than their equivalent and they offer some great features. If you work a lot with js, css or html you should definitely check these.
http://jashkenas.github.com/coffee-script/
http://compass-style.org/
http://haml-lang.com/

Three.js
The project is using webgl to render the scene. Three.js was obviously the 3d engine of choice. It provide many features like shaders, geometries, postprocessing filters and so on. A great project :)
https://github.com/mrdoob/three.js/

node.js
This is a javascript server. It's used to automatically compile the source files and serve them. I could have used some simple bash script to compile the files but found it would be interesting to use this for this "node" based project.
http://nodejs.org/

Require.js
The application was quickly growing and I needed a way to make it more modular. It's still not used at his full potential in threenode.js but already help make the code cleaner. One nice feature is that it's possible to require text files, or in this case little html template files.
http://requirejs.org/

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

http://processingjs.org/

Processing.js is the sister project of the popular Processing visual programming language, designed for the web. Processing.js makes your data visualizations, digital art, interactive animations, educational graphs, video games, etc. work using web standards and without any plug-ins. You write code using the Processing language, include it in your web page, and Processing.js does the rest. It's not magic, but almost.

Is it possible to use a canvas to do 3d?; 3d web development with Google's O3D ; a canvas application which can read the contents of a 3D .OBJ file and display the results in real-time.

the new O3D Project Hosting site: http://code.google.com/p/o3d/

3d web development with Google's O3D: WebGL implementation of O3D

The old O3D plugin API is Deprecated:  
http://code.google.com/apis/o3d/docs/samplesdirectory.html

Originally built as a browser plug-in, this new implementation of O3D is a JavaScript library implemented on top of WebGL.


Introduction
The WebGL implementation of O3D is a JavaScript library built on top of WebGL. It implements a subset of the original O3D plug-in API, which includes all O3D interfaces that have an analogous interface in WebGL. For example, the render graph, transform graph, and texture classes are all implemented in the WebGL implementation of O3D. Animation, parameter operations, 2D canvas, and file I/O classes are not included. See Functional Groupings of O3D for details.

One major difference between the two implementations is that the WebGL implementation of O3D uses the GLSL shader required by WebGL. This release includes a Cg-to-GLSL converter script to aid in this conversion.



Here (http://code.google.com/p/o3d/) you can:

• Download the WebGL implementation of O3D.
• Read about how to convert your O3D plug-in application to the WebGL implementation of O3D.
• Check out the samples.
• Read about how to load COLLADA files into the WebGL implementation of O3D.
• Browse the source code.
• File bugs and submit patches.

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

O3D is a good development environment for 3d.
You will be surprised by how easy it was compared to OpenGL.
For example, entire 3d models could be loaded using a single command, while as with OpenGL you would have had to write your own low-level model loader.
Not to mention, Google made integration with Sketchup easy.

ex-ample:

http://src.chromium.org/viewvc/chrome/trunk/src/o3d/samples/o3d-webgl-samples

• Hello Cube, v. 2
• Spinning Cube, v. 2 (minimal changes)
• Primitives, v. 2

http://code.google.com/p/o3d/wiki/SampleCodeWalkthrough

Converting and Loading COLLADA Models:
http://code.google.com/p/o3d/wiki/ColladaConverter

-------
"Is it possible to use a canvas to do 3d?"
3d and HTML5's Canvas element together...
The canvas element is becoming more widely adopted than O3D.  It could draw polygons.

The future of 3d web applications lies with WebGL : http://www.khronos.org/webgl/.
Using HTML5's canvas unlike WebGL: canvas will be able to run on any HTML5 enabled browser (software rendered) while as WebGL runs on computers with hardware acceleration.

Ex-ample:

3D OBJ Viewer – 3D, Applications - Canvas Demos:
http://www.canvasdemos.com/2010/05/05/3d-obj-viewer/
http://www.canvasdemos.com/userdemos/toxicgonzo/3dobjviewer.html
a canvas application which can read the contents of a 3D .OBJ file and display the results in real-time.
This particular demo takes Blender’s monkey model and rotates it. It also displays how many FPS in the upper left corner.

This program:

1.  A 3d .OBJ object is stored as a string in the javascript file
2. The .OBJ string is decoded to hold information about vertex position and which vertices form the polygons
3.  Every polygon is assigned a random color
4. Matrix multiplication happens: Convert object from model space -> camera space -> clip space -> screen space
5.  Z-sort the polygons based on the polygon’s centroid
6. Draw the polygons from back to front

webGL skull (1162 faces) rendering.

http://www.chromeexperiments.com/detail/webgl-skull/?f=


Low poly normal mapped model exported from Blender, rendered in Three.js using Blinn-Phong shader with baked ambient occlusion texture, one directional and one point light.

Skull model by Daniel FR Gordillo:
http://www.blendswap.com/3D-models/misc-objects/craneo/
1162 faces, 924 KB total (including textures)
Technology: JavaScript, WebGL, Three.js, Blender

7th Symposium on Trends in Unstructured Mesh Generation; 2009

MeshTrends VII

7th Symposium on Trends in Unstructured Mesh Generation

July 16-19, 2009 Columbus, Ohio

Symposium Schedule

The list of speakers and titles is now availble here (pdf).

The Symposium on Trends in Unstructured Mesh Generation brings together a wide variety of disciplines for the exchange of technical information related to unstructured mesh generation. It is a symposium traditionally held in conjunction with the national and international computational mechanics congresses. The following is a list of previous MeshTrends symposia:

MeshTrends I	1997	Joint ASME/ASCE/SES Summer Meeting	Northwestern University
MeshTrends II	1999	5th US National Congress on Computational Mechanics	University of Colorado, Boulder
MeshTrends III	2001	6th US National Congress on Computational Mechanics	Dearborn, Michigan
MeshTrends IV	2003	7th US National Congress on Computational Mechanics	Albuquerque, New Mexico
MeshTrends V	2006	7th World Congress on Computational Mechanics	Los Angeles, California
MeshTrends VI	2007	9th US National Congress on Computational Mechanics	San Francisco, California

Scope of Symposium

Automatic unstructured mesh generation continues to be a vital technology in computational field simulations. As computing technology continues to advance and modeling requirements become more precise, automatic mesh generation techniques must rise to fulfill ever-increasing and diverse expectations. This symposium is a forum for exploring and synthesizing many of technologies needed to develop a computational grid suitable for simulation.

All abstracts related to geometry and mesh generation for computational simulation are welcome. In this symposium we are soliciting, in particular, advancements and trends from academics and industry in the following areas:

• Hexahedral mesh generation: including theoretical foundations and new algorithms for automatic all-hex methods.
• Meshing tools and applications: including commercial meshing tools and their application to current problems in industry.
• Multiphysics meshing issues: including tools and methods for managing meshing and geometry for mutiscale, mutliphysics applications.
• Infrastructure and tools for meshing: including APIs and tools for managing and interfacing meshing tools.
• Adaptive meshing tools and applications: including tools and methods for adaptively modifying mesh and geometry based on run-time results or optimization parameters.
• Meshing and geometry for geophysics applications: including geometry and meshing technologies for subsurface modeling and simulation.
• Meshing and geometry for biomedical applications: including geometry and meshing technologies for biomedical applications.
• Geometry repair and improvement for mesh generation: including tools and methods for characterizing dirty geometry and improvement techniques for mesh generation.
  

Abstract Submission

Abstracts are required for the conference and will be included in the conference proceedings. A one-page abstract must be submitted electronically through the USNCCM10 website. The deadline for abstracts was Feb 28, 2009.

Once you enter the abstract submission webpage at http://usnccm-10.eng.ohio-state.edu/abstractsub.html you will be asked to create a login and password. Note the number of the Trends in Unstructured Mesh Generation Symposium is 2.18.2. You will need to use that number in your abstract file name to ensure your abstract is submitted to the correct symposium (see below for more information).

The following information regarding abstract submission also appears on the abstract submission webpage.

Abstract Submission Information and Related Policies (Please read carefully) Each paid registrant at USNCCM-10 will be limited to one presentation. Please identify the presenting author in your abstract by underlining the last name. Participants may be the author of multiple abstracts. However, he/she will be limited to being the presenting author on only one abstract. Please submit a one-page abstract using the format shown here. Your submission must be in a pdf format to be submitted. Your filename should contain the last names of the authors as they appear in the title, followed by the number of your selected minisymposium. Example: Ghosh Joseph 2.18.2.pdf (Note the Trends in Unstructured Mesh Generation Symposium is number 2.18.2) The following specifications should be followed in the abstract pdf. Each abstract should be limited to 1-page. The file size should be limited to 1 MB. Recommended font: Times Roman; Title font size: 14 Boldface; Abstract font size: 12; Line spacing: single space; Number of References: maximum of 2. Important: Save your logon and password; you may use them again when registering to attend the Congress in Spring 2009. Once submitted, your abstract will be reviewed by the organizers of your selected mini-symposium. You will receive an email notification that the abstract has been received. We anticipate that notification of your abstract acceptance or rejection will be made 30 days after the close of abstract submission as published on the conference web site.

Paper Submission

As part of this symposium, full papers will be solicited from the accepted presentations for inclusion in a peer-reviewed special journal edition of Engineering With Computers. Publication solicitation will be based on the interest of the participating authors and the technical merit of the presentation. Invitations for paper submissions will be made following the USNCCM.

Important Dates

Abstract submission opened on USNCCM web site	October 20, 2008
Deadline for receipt of one-page abstracts	February 28, 2009
Notification of abstract acceptance	March 15, 2009
Deadline for early registration	May 1, 2009
USNCCM X technical program	July 16-19 2009

Symposium Organizers

Steven J. Owen, Ph.D.
Computational Modeling Sciences Department
Sandia National Laboratories
Albuquerque, New Mexico, U.S.A.
Phone: (505) 284-6599
Email: sjowen@sandia.gov

Mark S. Shephard, Ph.D.
Director, Scientific Computation Research Center
Rensselaer Polytechnic Institute
Troy New York, U.S.A.
Phone: (518) 276-6795
Fax: (518) 276-4886
Email: shephard@scorec.rpi.edu

Matthew L. Staten
Carnegie Mellon University and
Sandia National Laboratories
Email: mlstate@sandia.gov

Additional Information

Additional information on the conference can be found at:
http://usnccm-10.eng.ohio-state.edu/index.html

http://www-users.informatik.rwth-aachen.de/~roberts/meshgeneration.html

Mesh Generation & Grid Generation on the Web

The aim of this document is to provide information on mesh and grid generation: people working in the field, research groups, books and conferences. It is maintained by Robert Schneiders.

Mesh generation is an interdisciplinary area, and people from different departments are working on it: Mathematicians, computer scientists, engineers from many disciplines. Despite the fact that surprisingly many people are active in the field, often there are few contacts between researchers. The aim of this page is to improve communication between research groups and to help people to get an overview of the field.

The page is organized as follows:

<!-- -->

People and research groups: Info on meshing research at universities, companies, government labs etc. List of people: A directory of people working on mesh generation. Latest news: What's up in mesh generation. Software: A list of programs, both public domain and commercial. Conferences: Information on conferences, summerschools, short courses etc. Literature: Books, reviews, online sources and course materials. Open positions: Career opportunities for people with background in mesh generation. Information on related topics: Pages with information on CFD, scientific computing, computational geometry and other fields related to mesh generation.

<!-- --> Service for frequent readers: You can find all entries, sorted by time of insertion, here (there is also an archive page). <!-- Click here to see the latest updates (there is also an archive page). -->

Research on mesh generation is abundant, and I don't claim to give a complete overview. In order make this page a useful service for the mesh generation community, I need help from other people. So if you are interested in getting put on the list, or if you have any comments or hints on other sources of information on mesh generation in the net, please send me an email (robert.schneiders@arcor.de).

A valuable source of information is the Meshing Research Corner, a comprehensive database with literature on mesh generation. It is maintained by Steve Owen.

ParaView is an application framework as well as a turn-key application; Modeling software OpenFOAM and ParaView ; mesh processing: generation, manipulation, conversion

ParaView is an open source, multi-platform data analysis and visualization application. It has a client-server architecture to facilitate remote visualization of datasets

It is an application built on top of the Visualization Tool Kit (VTK) libraries.

The ParaView code base is designed in such a way that all of its components can be reused to quickly develop vertical applications. This flexibility allows ParaView developers to quickly develop applications.

Input/Output and File Format

• Supports a variety of file formats including: VTK (new and legacy, all types including parallel, ascii and binary, can read and written).
• Various polygonal file formats including STL and BYU (by default, read only, other VTK writers can be added by writing XML description).
• Many other file formats are supported. See ParaView Readers and ParaView Writers for a full list.

CMake is a family of tools designed to build, test and package software. CMake is used to control the software compilation process using simple platform and compiler independent configuration files. CMake generates native makefiles and workspaces that can be used in the compiler environment of your choice. ParaView utilizes CMake for the software compilation process.

--------

ParaView is used as the visualization platform for the Modeling software OpenFOAM (Open Field Operation and Manipulation).It is primarily a C++ toolbox for the customisation and extension of numerical solvers for continuum mechanics problems, including computational fluid dynamics (CFD). It comes with a growing collection of pre-written solvers applicable to a wide range of problems.

First major general-purpose CFD package to use polyhedral cells. This functionality is a natural consequence of the hierarchical description of simulation objects.

OpenFOAM compares favourably with the capabilities of most leading general-purpose commercial closed-source CFD packages. It relies on the user's choice of third party pre- and post-processing utilities, and ships with:

• a plugin (paraFoam) for visualisation of solution data and meshes in ParaView.
• a wide range of mesh converters allowing import from a number of leading commercial packages
• an automatic hexahedral mesher to mesh engineering configurations

OpenFOAM was conceived as a continuum mechanics platform but is ideal for building multi-physics simulations.

OpenCFD develop OpenFOAM in the Linux/UNIX operating system because: we believe it is the best platform for this kind of high end simulation code development and operation; Linux is efficient, robust, reliable and flexible and undergoes rapid development and improvement; Linux is open source, like OpenFOAM; Linux is very effective for parallel operation on Beowulf clusters.

OpenFOAM is open source software so people can freely compile it on any operating system they choose. Most OpenFOAM users are running Linux, so this site offers the download of binaries for selected Linux systems.

As the present time we are unaware of any binary distributions for Windows or MacOSX. However, ports to these operating systems have been the subject of debate on the OpenFOAM discussion site, which may provide the best source of information on the matter.

http://www.opencfd.co.uk/openfoam/

OpenFOAM uses finite volume numerics to solve systems of partial differential equations ascribed on any 3D unstructured mesh of polyhedral cells.

Mesh generation

OpenFOAM applications handle unstructured meshes of mixed polyhedra with any number of faces: hexahedra, tetrahedra, degenerate cells, basically anything.

Mesh generation is made simple by the fact that a cell is simply represented as a list of faces and a face as a list of vertices: this makes mesh handling very easy even for complex meshes with, say, embedded refinement or complex shapes near the boundary.

OpenFOAM is supplied with the following mesh generator tools that run in parallel.

Mesh generation tools
blockMesh	A multi-block mesh generator
extrude2DMesh	Takes 2D mesh (all faces 2 points only, no front and back faces) and creates a 3D mesh by extruding with specified thickness
extrudeMesh	Extrude mesh from existing patch (by default outwards facing normals; optional flips faces) or from patch read from file
snappyHexMesh	Automatic split hex mesher. Refines and snaps to surface

The main mesh generators cover two extremes: snappyHexMesh, that can mesh to complex CAD surfaces; blockMesh a simple file-driven block mesh generator.

Mesh manipulation

OpenFOAM is supplied with several utilties that perform mesh checking and manipulation. The full list of utilties is given below

Mesh manipulation
attachMesh	Attach topologically detached mesh using prescribed mesh modifiers
autoPatch	Divides external faces into patches based on (user supplied) feature angle
cellSet	Selects a cell set through a dictionary
checkMesh	Checks validity of a mesh
createBaffles	Makes internal faces into boundary faces. Does not duplicate points, unlike mergeOrSplitBaffles
createPatch	Utility to create patches out of selected boundary faces. Faces come either from existing patches or from a faceSet
deformedGeom	Deforms a polyMesh using a displacement field U and a scaling factor supplied as an argument
faceSet	Selects a face set through a dictionary
flattenMesh	Flattens the front and back planes of a 2D cartesian mesh
insideCells	Picks up cells with cell centre ’inside’ of surface. Requires surface to be closed and singly connected
mergeMeshes	Merge two meshes
mergeOrSplitBaffles	Detects faces that share points (baffles). Either merge them or duplicate the points
mirrorMesh	Mirrors a mesh around a given plane
moveDynamicMesh	Mesh motion and topological mesh changes utility
moveEngineMesh	Solver for moving meshes for engine calculations.
moveMesh	Solver for moving meshes
objToVTK	Read obj line (not surface!) file and convert into vtk
pointSet	Selects a point set through a dictionary
refineMesh	Utility to refine cells in multiple directions
renumberMesh	Renumbers the cell list in order to reduce the bandwidth, reading and renumbering all fields from all the time directories
rotateMesh	Rotates the mesh and fields from the direcion to the direction
setSet	Manipulate a cell/face/point set interactively
setsToZones	Add pointZones/faceZones/cellZones to the mesh from similar named pointSets/faceSets/cellSets
splitMesh	Splits mesh by making internal faces external. Uses attachDetach
splitMeshRegions	Splits mesh into multiple regions
stitchMesh	’Stitches’ a mesh
subsetMesh	Selects a section of mesh based on a cellSet
transformPoints	Transforms the mesh points in the polyMesh directory according to the translate, rotate and scale options
zipUpMesh	Reads in a mesh with hanging vertices and zips up the cells to guarantee that all polyhedral cells of valid shape are closed

Mesh motion

OpenFOAM adopts a novel approach to mesh motion by defining it in terms of the boundary motion which is extremely robust.

The solver need only define the the motion of the boundary and everything else will be done automatically. The open architecture of OpenFOAM solver codes allows quick and efficient implementation: mesh motion can be based on any solution variable, either local or integrated and by dynamically adjusted during the run.

Mesh motion is also transparently integrated with top-level models: the model writer does not see the additional complexity, which is conveniently packaged within the discretisation operators.

For examples of automated mesh motion in OpenFOAM, see Solutions

Mesh conversion

OpenFOAM accepts meshes generated by any of the major mesh generators and CAD systems. Listed below are converter utlities for the major commercial mesh generators. Note that it is also possible to import the meshes from most general purpose mesh generators since they will export in a format read by one of the converters.

Mesh converters
ansysToFoam	Converts an ANSYS input mesh file, exported from I-DEAS, to OPENFOAM®format
cfx4ToFoam	Converts a CFX 4 mesh to OPENFOAM®format
fluent3DMeshToFoam	Converts a Fluent mesh to OPENFOAM®format
fluentMeshToFoam	Converts a Fluent mesh to OPENFOAM®format including multiple region and region boundary handling
foamMeshToFluent	Writes out the OPENFOAM®mesh in Fluent mesh format
foamToStarMesh	Reads an OPENFOAM®mesh and writes a PROSTAR (v4) bnd/cel/vrt format
gambitToFoam	Converts a GAMBIT mesh to OPENFOAM®format
gmshToFoam	Reads .msh file as written by Gmsh
ideasUnvToFoam	I-Deas unv format mesh conversion
kivaToFoam	Converts a KIVA grid to OPENFOAM®format
mshToFoam	Converts .msh file generated by the Adventure system
netgenNeutralToFoam	Converts neutral file format as written by Netgen v4.4
plot3dToFoam	Plot3d mesh (ascii/formatted format) converter
polyDualMesh	Calculate the dual of a polyMesh. Adheres to all the feature and patch edges
sammToFoam	Converts a STAR-CD SAMM mesh to OPENFOAM®format
star4ToFoam	Converts a STAR-CD (v4) PROSTAR mesh into OPENFOAM®format
starToFoam	Converts a STAR-CD PROSTAR mesh into OPENFOAM®format
tetgenToFoam	Converts .ele and .node and .face files, written by tetgen
writeMeshObj	For mesh debugging: writes mesh as three separate OBJ files which can be viewed with e.g. javaview

FEM, finite element method, PDE, partial differential equation;mathematica

http://library.wolfram.com/infocenter/MathSource/4478/
This package allows to solve second order elliptic differential equations in two variables:

div(a*grad u) - b*u = f in the domain domain u = gD Dirichlet boundary conditions on first part of boundary a*du/dn = gN Neumann condition on the other part of the boundary

If the functions a, b f, gD and gN are given, then a numerical approximation is computed, using the method of finite elements. To generate meshes the programm EasyMesh can be used.

MESH GENERATION;mathematica

http://library.wolfram.com/infocenter/Conferences/5771/

The upcoming release of Mathematica (2005) establishes the next level of sophistication in scientific function and data visualization. All the plotting functions have been rewritten to incorporate state-of-the-art algorithms in numeric and algebraic analysis, mesh generation, and computational geometry. Arbitrary-precision function evaluation, region, mesh overlays, piecewise functions, and singularities handling are just part of the new functionality. All of this makes Mathematica a unique system for function and data visualization with features not found in any other algebraic computational system. In this presentation we will give a general description and present several examples of the new capabilities.

mesh generation, literate programming, matlab; mathematica; persson

http://library.wolfram.com/infocenter/MathSource/5475/

This Mathematica notebook is an effort to transcribe the MATLAB code of a 2-D mesh generation algorithm as described explicitly in Persson and Strang's paper [1]. The goal is to make the algorithm executable in Mathematica so that its users can also experiment with the algorithm.

Since the algorithm was expressed very clearly from their original paper [1] including the MATLAB code, which is a perfect example of literate programming in MATLAB, it is pretty easy to translate the MATLAB code "literally" into Mathematica. Such translation is virtually always possible in either direction even without human interference. And such a Rosetta Stone kind of translation might be useful if one species of people coding in either MATLAB or Mathematica were to disappear, future generations would still be able to rediscover one programming language by reading its interpretation in the other one.

However, it is so tempting to present the literate programming capability of Mathematica by following its general principles; that is, (a) documentation mingles with code and both get pretty-printed; (b) shuffle code pieces for human readability. I decided to transcribe the code manually.

The original MATLAB code was documented as 8 steps (sections) in sequential order, which is easy to follow because the ideas behind the code were explained beforehand in early parts of the paper. So it is recommended that you read part 1 and 2 of the original paper. Instead of following the MATLAB code literally in 8 steps, this notebook breaks the code pieces apart and examines each of them separately.

Mesh Generation using Matlab ; mesh generator, 2D and 3D;

Mesh Generation using Matlab

These days most of the research in the field of fluids, structures, porous media, brain computer interfacing you name it, uses numerical simulations. Reason: It is much cheaper and many times faster compared to experiments. Mesh Generation forms an integral part of numerical analysis/simulation.

Although, there are plenty of commercial softwares based on Finite Element Methods and Finite Volume Methods like COMSOL, FLEUNT, ANSYS, NUMECA and many more with exceptional Mesh/Grid Generation features (amira).

But many times its difficult to use the meshes generated by these softwares which suits to your particular simulation need.

Reason: Many of the exsiting software don't have this feature where you can create a mesh and use it some which have such features requires you to do some complicated modifications in your code to import these meshes. There are although loads of mesh generator available some of which open source and free to download.

But, then again problem comes does these free source code suits your purpose. I encountred this problem over the last couple of months.

I am doing research in the field of Petroleum Reservoir Simulation and I need to test a lot of numerical examples on different sorts of meshes/grids in 2 and 3D. I do most of my simulation work in MATLAB, some people might argue that MATLAB is slow and all sorts of reason about other programming languages are faster like C++ and Fortran. I don't deny that fact but on the other hand the library of existing function which matlab has is amazing and its Array handling feature and sparse code it amazing too. The only and important reason I use MATLAB is its capability to handle array operations.

{In my simulation code I have to solve at times 9 simulatneous equations in 2D and 27 equations in 3D, which maximizes use of array operations. I also frequently use MAPLE to do my algebra and other good thing about MATLAB is that I can directly import the MAPLE algebra in Array Format into MATLAB which suits my purpose}

Now comming to the meshes in MATLAB, try doing a google on 'meshes in MATLAB' or 'grid generation in 2 and 3D in MATLAB', a invested a lot of time to find some unseful source code in matlab searching on google groups etc the only useful package I found was by Per-Olof Persson titled 'DistMesh - A Simple Mesh Generator in MATLAB'.

YES!

see: http://persson.berkeley.edu/software.html

No doubt its an amazing piece of work but again it didnt realy suits my purpose. The reason being I needed unstructured meshes of different element types in 3D like prisms, hex, tetra and pyramids. In 2D also I needed meshes which are boundary aligned to control volume and are matching to the underlying medium.

So, What next ? I started from scratch and now I have come up with stand alone code in MATLAB which has functionality to create different kind of meshes in 2D and 3D. These are structured and Unstructured meshes, perturbed and bondary aligned too. If any one is in need of such meshes in 2 and 3D please have a look at:

www.mayurpal.com

Then you can drop me an email and I will get back to you and will help you and if required will also provide you with the source code if it suits your purpose.

Posted by C2EC- Swansea

http://meshgeneration.blogspot.com/2006/10/mesh-generation.html

Have you seen this MATLAB mesh generator:

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=10307&objectType=FILE

iso2mesh: a Matlab/Octave-based mesh generator: Examples Segmented images

iso2mesh: a Matlab/Octave-based mesh generator: Examples Segmented images

iso2mesh: a Matlab/Octave-based mesh generator: Doc/Examples: "Example Scripts for Using iso2mesh

1. Demonstration for surf2mesh
2. Demonstration 1 for vol2mesh
3. Demonstration 2 for vol2mesh
4. Demonstration 3 for vol2mesh

We provided a couple of example scripts to demonstrate the basic utilities of this toolbox. They can be found at sample/ subdirectory from iso2mesh's installation directory. In this page, we will briefly describe the input/output for each of the example and show you the process of the mesh creation using iso2mesh. Most the examples come from author's realistic research problems. For example, the surf2mesh example came from vessel network modeling; demo1 for vol2mesh came from small animal florescences imaging, and demo3 of vol2mesh originated from a combined image reconstruction of brain functional imaging using MRI, MEG and diffuse optical tomography.

1. Demonstration for surf2mesh
The sample script 'demo_surf2mesh_ex1.m' is to demonstrate how to use surf2mesh for meshing a simple open-ended surface. The surface is saved in 'surfmesh_demo.mat'. Loading this mat file, you will find two variables, f, and v, they record a tube surface mesh. One can view this surface by

trisurf(f,v(:,1),v(:,2),v(:,3)); axis equal;

The surface looks like upload:tube_surf.png

then, we call surf2mesh to produce volumetric mesh from this surface. This includes 3 sub steps

1. in the first step, the fine surface, f and v, are re-sampled to only keep 10% of the original elements
2. because the surface has openings, in the second step, iso2mesh added a cubic frame outside the tube, and punch openings on the bounding box and merge with the tube surface to create a closed surface
3. in the third step, a volumetric mesh is generated within the enclosed volume bounded by the bounding box and the tube.

After step 1, the surface mesh look like upload:tube_surf_sampled.png
After step 2, the surface mesh look like upload:tube_mesh.png
After step 3, the volume mesh look like upload:tube_mesh2.png

I want to mention that at this point, iso2mesh can only mesh surfaces with their edges (borders of the open regions) strictly located within a bounding box surface. It can even handle the complex cases where a opening edge contour located on multiple sides of the bounding box.

2. Demonstration 1 for vol2mesh
The first vol2mesh demonstration is to create an FEM mesh from a CT scan of a rat-head. In this example, the volumetric data is a thresholded CT image of a rat. In a vol2mesh session, we first wrap the volume with a layer of zeros, in order to make the non-zero region close. Then we extract the boundary of the non-zero regions by calling Matlab's isosurface(). The created iso-surface, in the resolution of voxel sizes, will be re-sampled to specified mesh density. Before and after the simplification, we perform a mesh validation and repair process to avoid the appearance of 'non-manifold' nodes. Once the surface mesh is completed, we produce volumetric mesh within the region bounded by the re-sampled surface mesh.

The voxelated rat head binary image looks like upload:mouse_vox.png
The re-sampled rat surface mesh is shown below upload:mouse_mesh.png
The cross-section of the 3D volumetric mesh looks like upload:mouse_head.png
3. Demonstration 2 for vol2mesh
In this example, we try to make a little bit fun using iso2mesh while showing off the extreme flexibility and simplicity of making a complex mesh using this toolbox. We will make a mesh on a slab-shaped object with 'ISO2Mesh' etched on one of the surface, just like the header image of this document showing.

To make this 'object', we first create a black-white image, 'iso2mesh_bar.tif', with an image editor. This image is then read as a 2D array into matlab, we then invert the 0-1 values, and stack it to 30 layers to make a etched patten, then, append another 30 layers of pure 1 image at the end of the stack as the base. Now we have this binary array with etched letter patterns.

The original 2D binary pattern is upload:iso2mesh_bar_2d.png
and the 3D object surface look like upload:iso2mesh_bar_3d_vox_small.png

We pass this array to vol2mesh, and set the surface simplification ratio to 10%, and maximum volume of 100 (in voxel unit), after going through similar processes as in demo1, we will get a volumetric mesh just like what you see from the header bar image.

The generated surface mesh is upload:iso2mesh_bar_1.png

4. Demonstration 3 for vol2mesh
Demonstration 3 is another challenging example showing you how to apply iso2mesh in complex geometries such as head and brain imaging. The 3D head and brain images were extracted from a real head scan from MRI. To save space, we convert them into a binary format as head.tif and brain.tif, respectively, where the brain is simple a segmented subregions from the head scan.

The raw binary head image is shown as upload:head_orig.png

Reading the head binary image into matlab and plot a slice, you will find there are lots of gaps and openings in the volume. As the dimensions of these gaps are close to the sale of a few voxel, it will give iso2mesh a really bad time trying to mesh with these gaps and holes. Therefore, we first perform a binary volume repairing process: deisland3d, to fill the holes and remove the isolated regions form the image. The cleaned up version of the head image can be found at here. The brain volume is clean by itself, so, we do not need to do hole-filling and island-removing.

The head image after hole-filling/island-removing is upload:head_fill.png

In this case, we will make a 2-level surface meshes: we will make 3D FEM mesh to fill both the regions between the head surface and brain surface, as well as the volume within the brain surface, and we want to make sure that the surface representation of the brain appears exactly in the resulting volumetric mesh (i.e., there is a conformal surface layer in the FEM mesh to match the brain surface). To do this, we first add the brain binary image on top of the head image, to create a two level image: voxels with value of 0 represent air, voxels of 1 represent the head region outside the brain, and voxels of value 2 represent the brain region. We pass this 3-valued array, fullimg, to vol2mesh, vol2mesh will first extract the iso-surface the interface between 0 and 1, performing mesh repairing and simplification; and then the interfaces between 1 and 2. After both surfaces complete, we will produce volumetric mesh for this mesh cluster. The resulting mesh will have a exterior surface conforms to the head contour, and an internal surface to conform to the brain surface.

The 3-level head and segmented brain image upload:head_seg.png
The brain surface mesh is upload:head_brain_meshed.png
The head surface mesh is upload:head_surface.png
The produced volumetric mesh cross-cut view 1 upload:head_mesh_cut1.png
The produced volumetric mesh cross-cut view 2 upload:head_mesh_cut2.png"

MATLAB Central - How to create 3D mesh model?

MATLAB Central - How to create 3D mesh model?

MATLAB Central - Newsreader - How to create 3D mesh model?: "Thread Subject: How to create 3D mesh model?

Subject: How to create 3D mesh model?

From: Tong

Date: 14 Jul, 2009 19:55:03

Message: 1 of 6
Reply to this message
Add author to My Watch List
View original format
Flag as spam

I have segmented meniscus images from MRI that is created in about 3mm slices. How would I combine these slices together to create a 3D model of the meniscus?

Subject: How to create 3D mesh model?

From: Luigi Giaccari

Date: 14 Jul, 2009 20:49:03

Message: 2 of 6
Reply to this message
Add author to My Watch List
View original format
Flag as spam

Please send me that models of yours, I am plannig to build a surface recostructor for sliced cloud. Send to : giaccariluigi@msn.com

In the mean time look for:

http://www.mathworks.com/matlabcentral/fileexchange/22185
http://giaccariluigi.altervista.org/blog/

and related

Subject: How to create 3D mesh model?

From: Brad Henrie

Date: 17 Jul, 2009 21:45:18

Message: 3 of 6
Reply to this message
Add author to My Watch List
View original format
Flag as spam

'Tong ' <celticbaseball06@gmail.com> wrote in message <h3inqn$ni5$1@fred.mathworks.com>...
> I have segmented meniscus images from MRI that is created in about 3mm slices. How would I combine these slices together to create a 3D model of the meniscus?

First place all of your slices into a 3-d matrix. This will give you a cube of data. You can then view it from multiple planes by using this format variable(:,:,a) where a is the slice position in a direction directly into your displayed image. Using the same format you can display other planes variable(:,a,:). Converting your image to greyscale will allow you to display it using implay.

I'm sure that since you are working with MRI you have access to the image processing toolbox.

While viewing images in a plane where the pixels are not square you need to scale your image. (if you have a 3x3x5 voxel and display the 3x5 pixel representation) Also remember your slice separation if you don't have 3-d k-space.

Subject: How to create 3D mesh model?

From: Image Analyst

Date: 18 Jul, 2009 04:02:35

Message: 4 of 6
Reply to this message
Add author to My Watch List
View original format
Flag as spam

'Tong ' <celticbaseball06@gmail.com> wrote in message <h3inqn$ni5$1@fred.mathworks.com>...
> I have segmented meniscus images from MRI that is created in about 3mm slices. How would I combine these slices together to create a 3D model of the meniscus?
----------------------------------------
I'm not sure what you mean by 'model,' but you can combine 2D images together to form a 3D image by using the cat(3, slice1, slice2, slice3, slice4, slice5,......) function.

Subject: How to create 3D mesh model?

From: Tong

Date: 20 Jul, 2009 18:36:02

Message: 5 of 6
Reply to this message
Add author to My Watch List
View original format
Flag as spam

'Image Analyst' <imageanalyst@mailinator.com> wrote in message <h3rhgr$of5$1@fred.mathworks.com>...
> 'Tong ' <celticbaseball06@gmail.com> wrote in message <h3inqn$ni5$1@fred.mathworks.com>...
> > I have segmented meniscus images from MRI that is created in about 3mm slices. How would I combine these slices together to create a 3D model of the meniscus?
> ----------------------------------------
> I'm not sure what you mean by 'model,' but you can combine 2D images together to form a 3D image by using the cat(3, slice1, slice2, slice3, slice4, slice5,......) function.

What about when I am using regions of interest, not images?

Subject: How to create 3D mesh model?

From: fabio freschi

Date: 20 Jul, 2009 21:06:01

Message: 6 of 6
Reply to this message
Add author to My Watch List
View original format
Flag as spam

you can try iso2mesh in FE
fabio"

matlab .m .ppt course : introduction to the theory of manifolds; mesh generation ; mesh processing

take care : too many bugs...

This course is an introduction to the  theory of manifolds. Manifold models arise in various area of mathematics, image processing, data mining or computer science. Surfaces of arbitrary dimension can be used to model non-linear datasets that one encounters in modern data processing. Numerical methods allow to exploit this geometric non-linear prior in order to extract relevant information from the data. These methods include in particular local differential computations (related to the Laplacian operator and its variants) and global distance methods (related to geodesic computations). In this course, you will learn how to perform differential and geodesic computations on images, volumes, surfaces and high dimensional graphs.

The course includes a set of Matlab experiments. These experiments give an overview of various tasks in computer vision, image processing, learning theory and mesh processing. This includes computation of shortest paths, voronoi segmentations, geodesic delaunay triangulations, surface flattening, dimensionality reduction and mesh processing.

Ref.

http://www.mathworks.de/matlabcentral/fileexchange/loadFile.do?objectId=13464&objectType=file



---

Lecture 0 - Basic Matlab Instructions

Abstract : Learn the basic features of Matlab, and learn how to load and visualize signals and images.

Basic Matlab commands.

• Create a variable and an array

% this is a comment a = 1; a = 2+1i; % real and complex numbers b = [1 2 3 4]; % row vector c = [1; 2; 3; 4]; % column vector d = 1:2:7; % here one has d=[1 3 5 7] A = eye(4); B = ones(4); C = rand(4); % identity, 1 and random matrices c = b'; % transpose

• Modification of vectors and matrices

A(2,2) = B(1,1) + b(1); % to access an entry in a vector or matrix b(1:3) = 0; % to access a set of indices in a matrix b(end-2:end) = 1: % to access the last entries b = b(end:-1:1); % reverse a vector b = sort(b); % sort values b = b .* (b>2); % set to zeros (threshold) the values below 2 b(3) = []; % suppress the 3rd entry of a vector B = [b; b]; % create a matrix of size 2x4 c = B(:,2); % to access 2nd column

• Advanced instructions

a = cos(b); a = sqrt(b); % usual function help perform_wavelet_transform; % print the help a = abs(b); a = real(b); a = imag(b); a = angle(b); % norm, real part and imaginary part of a complex disp('Hello'); % display a text disp( sprintf('Value of x=%.2f', x) ); % print a values with 2 digits A(A==Inf) = 3; % replace Inf values by 3 A(:); % flatten a matrix into a column vector max(A(:)); % max of a matrix M = M .* (abs(M)>T); % threshold to 0 values below T.

• Display

plot( 1:10, (1:10).^2 ); % display a 1D function title('My title'); % title xlabel('variable x'); ylabel('variable y'); % axis subplot(2, 2, 1); % divide the screen in 2x2 and select 1st quadrant

• Programmation

for i=1:4 % repeat the loop for i=1, i=2, i=3 et i=4 disp(i); % make here something end i = 4; while i>0 % while syntax disp(i); % do smth i = i-1; end

Load and visualize signals and images

• Load and plot a signal: (function load_signal.m should be in the toolbox of each course)

f = load_signal('Piece-Regular', n); % signal of size n plot(f);

• Load and display an image (download function load_image.m should be in the toolboxes)

I = load_image('lena'); imagesc(I); axis image; colormap gray(256);

Copyright © 2006 Gabriel Peyré

Lecture 1 - Active Contours and Level Sets

Abstract : The goals of this lecture is to use the level set framework in order to do curve evolution. The mean curvature motion is the basic tool, and it can be extended into edge-based (geodesic active contours) and region-based (Chan-Vese) snakes.

Setting up Matlab.

• First download the Matlab toolbox toolbox_fast_marching.zip. Unzip it into your working directory. You should have a directory toolbox_fast_marching/ in your path.
• The first thing to do is to install this toolbox in your path.

path(path, 'toolbox_fast_marching/'); path(path, 'toolbox_fast_marching/data/'); path(path, 'toolbox_fast_marching/toolbox/');

• Recompile the mex file for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_fast_marching/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

cd toolbox_fast_marching compile_mex; cd ..

Managing level set function.

• In order to perform curve evolution, we will deal with a distance function stored in a 2D image D. The curve will be embedded in the level set D=0. This curve will be evolved by modifying the image D. A curve evolution ODE can be replaced by an PDE on D. This allows to deal with topological changes when a curve split or two curves merge.

n = 200; % size of the image % load a distance function D0 = compute_levelset_shape('circlerect2', n); % type 'help compute_levelset_shape' to see other % basic curve you can load. % display the curve clf; hold on; imagesc(D); axis image; axis off; axis([1 n 1 n]); [c,h] = contour(D,[0 0], 'r'); set(h, 'LineWidth', 2); hold off; colormap gray(256); % do the union of two curves options.center = [0.15 0.15]*n; options.radius = 0.1*n; D1 = compute_levelset_shape('circle', n,options); imagesc(min(D0,D1)<0);

• During the curve evolution, the image D might become far from being a distance function. In order to stabilize the algorithm, one needs to re-compute this distance function.

% here we simulate a modification of the distance function [Y,X] = meshgrid(1:n,1:n); D = (D0.^3) .* (X+n/3); D1 = perform_redistancing(D); % display both the original and the new, % redistanced, curve (should be very close) ...

Mean Curvature Motion.

• In order to compute differential quantities (tangent, normal, curvature, etc) on the curve, you can compute derivatives of the image D.

% the gradient g0 = divgrad(D); % display the gradient (as arrow field with 'quiver', ...) ... % the normalized gradient d = max(eps, sqrt(sum(g0.^2,3)) ); g = g0 ./ repmat( d, [1 1 2] ); % display ... % the curvature K = d .* divgrad( g ); % display ...

• The mean curvature motion of the level sets of some image u is driven be the following equation.
  
  Implement this evolution explicitly in time using finite differences.

Tmax = 1000; % maximum time of evolution dt = 0.4; % time step (should be small) niter = round(Tmax/dt); % number of iterations D = D0; % initialization for i=1:niter % compute the right hand size of the PDE ... % update the distance field D = ...; % redistance the function from time to time if mod(i,30)==0 D = perform_redistancing(D); end % display from time to time if mod(i,30)=1 % display here ... end end

Curve evolution under the mean curvature motion (the background is the distance function D).

Edge-based Segmentation with Geodesic Active Contour (snakes + level set).

• Given a background image M to segment, one needs to compute an edge-stopping function E. It should be small in area of high gradient, and high in area of large gradient.

% load an image name = 'brain'; M = rescale( sum( load_image(name, n), 3) ); % display it ... % compute a smoothed gradient sigma = 4; % blurring size G = divgrad( perform_blurring(M,sigma) ); % compute the norm of the gradient d = ... % compute the edge-stopping function E = ... % rescale it so that it is in realistic ranges E = rescale(E,0.3,1);

• The geodesic active contour evolution is a mean curvature motion modulated by the edge-stopping function:
  
  Implement this evolution explicitly in time using finite differences.

...

Segmentation with geodesic active contours.

Region-based Segmentation with Chan-Vese (Mumord-Shah + level sets).

• The geodesic active contour uses an edge-based energy. It has lots of local minima and is very sensitive to initialization. In order to circumvent these drawbacks, one can use a region based energy like the Mumford-Shah functional. Re-casted into the level set framework, it reads:
  
  The corresponding gradient descent is the Chan-Vese active contour method:
  
  Implement this evolution explicitly in time using finite differences, when c1 and c2 are known in advanced.

% initialize with a complex distance function D0 = compute_levelset_shape('small-disks', n); % set parameters lambda = 0.8; c1 = ...; % black c2 = ...; % gray ...

Segmentation with Chan-Vese active contour without edges.

• In the case that one does not know in advance the constants c0 and c1, how to update them automatically during the evolution ? Implement this method.

Copyright © 2006 Gabriel Peyré

Lecture 2 - Front propagation
in 2D and 3D

Abstract : The goals of this lecture is to manipulate the fast marching algorithm in 2D and 3D. Application to shortest path extraction (e.g. road tracking and tubular structure extraction in medical images) and Voronoi cell segmentation is presented.

Setting up Matlab.

• First download the Matlab toolbox toolbox_fast_marching.zip. Unzip it into your working directory. You should have a directory toolbox_fast_marching/ in your path.
• The first thing to do is to install this toolbox in your path.

path(path, 'toolbox_fast_marching/'); path(path, 'toolbox_fast_marching/data/'); path(path, 'toolbox_fast_marching/toolbox/');

• Recompile the mex file for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_fast_marching/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

cd toolbox_fast_marching compile_mex; cd ..

Front propagation in 2D.

• We can now load an image and build a speed propagation function W. Note that area with W values near 0 (black pixels) are those in which the front propagate slowly, so that W is the inverse of the potential that weight the metric.

n = 128; name = 'cavern'; % other possibilities are 'mountain' and 'road2' W = load_image(name, n); W = rescale(W, 0.01, 1); % set up a reasonable range for the potential % display the weighting function clf; imagesc(W); colormap gray(256);

• After asking to the user the starting and ending points, we proceed to the front propagation. Note that the propagation terminates when the front reaches the ending point.

[start_points,end_points] = pick_start_end_point(W); options.end_points = end_points; [D,S] = perform_fast_marching_2d(W, start_points, options); % display the distance function clf; imagesc(D); D(I==Inf) = 0; % remove Inf values that make contour crash figure; contour(D,50);

• The shortest path is extracted by performing a gradient descent of D, starting from the ending point.

p = extract_path_2d(D,end_points, options); % display the path clf; plot_fast_marching_2d(W,S,p,start_points,end_points);

• Now that you know the basic commands, you can try other speed functions loaded from other files, and you can try more complicated potential than just the value of the image. You can try a potential based on the gradient of the image computed using grad(W).

Left : distance function for 3 different W maps. Center : the corresponding level set maps. Right : shortest path together with the explored area (in red).

3D Volumetric Shortest Paths.

• 
  
  Firs you need to load a 3D array of data M. Such a volumetric dataset is more difficult to visualize than a standard 2D image. You can render slices like M(:,:,i) or use a volumetric renderer such a vol3d. In order to do so, you need to set up a correct alpha mapping to make transparent some parts of the volume. If you rotate the data, then you need to re-render the volume using vol3d(h). You can download the 3D data set here.

% load the whole volume load brain1-crop-256.mat % crop to retain only the central part n = 100; M = rescale( crop(M,n) ); % display some horizontal slices imageplot(M(:,:,50)); % same thing in the other directions ...

Cross sections of the data-set.

• 
  
  Another, more efficient way, to render volumetric data is to display a semi-transparent volumetric image. This can be achieved using the function vol3d than render semi-transparent slices of the data orthogonal to the viewing direction. In order to do so, you need to set up a correct alpha mapping to make transparent some parts of the volume. If you rotate the data, then you need to re-render the volume using vol3d(h).

clf; h = vol3d('cdata',M,'texture','2D'); view(3); axis off; % set up a colormap colormap bone(256); % set up an alpha map options.center = ...; % here a value in [0,1] options.sigma = .08; % control the width of the non-transparent region a = compute_alpha_map('gaussian', options); % you can plot(a) to see the alphamap colormap bone(256); % refresh the rendering vol3d(h); % try with other alphamapping and colormapping ...

Volumetric rendering with different alpha-mapping. Each time the options.center value is increased.

• Geodesic distances can be computed on a 3D volume using the Fast Marching. The important point here is to define the correct potential field W that should be large in the region where you want the front to move fast. It means that geodesic will follow these regions. In order to do so, we will ask the user to click on a starting point in a given horizontal slice W(:,:,delta). The potential is then computed in order to be large for value of M that are close to the value at this starting point.

% ask to the user for some input point delta = 5; clf; imageplot(M(:,:,delta)); title('Pick starting point'); start_point = round( ginput(1) ); start_point = [start_point(2); start_point(1); delta]; % compute a potential that is high only very close % to the value of M at the selected point W = ...; W = rescale(W,.001,1); % perform the front propagation options.nb_iter_max = Inf; [D,S] = perform_fast_marching_3d(W, start_point, options); % display the results using vol3d ...

Exploration of the distance function using vol3d.

• In order to extract a geodesic, we need to select an ending point and perform a descent of the distance function D from this starting point. The selection is done by choosing a point of low distance value in the slice D(:,:,end-delta).

% extract a slice d = D(:,:,n-delta); % select the point (x,y) of minimum value of d % hint: use functions 'min' and 'ind2sub' ... [x,y] = ...; end_point = [x;y;n-delta]; ... % extract the geodesic by discrete descent path = compute_discrete_geodesic(D,end_point); % draw the path Dend = D(end_point(1),end_point(2),end_point(3)); D1 = double( D<=Dend ); clf; plot_fast_marching_3d(M,D1,path,start_point,end_point);

Exploration of the distance function using vol3d. The red surface indicates the region of the volume that has been explored before hitting the ending point.

Voronoi segmentation and geodesic Delaunay triangulation in 2D.

• With the same code as above, one can use multiple starting points. The function perform_fast_marching_2d returns a segmentation map Q that contains the Voronoi segmentation.

n=300; name = 'constant'; % other possibility is 'mountain' W = load_potential_map(name, n); m = 20; % number of center points % compute the starting point at random. start_points = floor( rand(2,m)*(n-1) ) +1; [D,Z,Q] = perform_fast_marching_2d(W, start_points); % display the sampling with the distance function clf; hold on; imagesc(D'); axis image; axis off; plot(start_points(1,:), start_points(2,:), '.'); hold off; colormap gray(256); % display the segmentation figure; clf; hold on; imagesc(Q'); axis image; axis off; plot(start_points(1,:), start_points(2,:), '.'); hold off; colormap gray(256);

Example of Voronoi cells (distance functions) obtained with a constant speed W (left) and the mountain map (right). Note how the cell on the left have polygonal boundaries whereas cells on the right have curvy boundaries.

• A geodesic Delaunay triangulation is obtained by linking starting points whose Voronoi cells touch. This is the dual of the original Voronoi segmentation.

faces = compute_voronoi_triangulation(Q); hold on; imagesc(Q'); axis image; axis off; plot(start_points(1,:), start_points(2,:), 'b.', 'MarkerSize', 20); plot_edges(compute_edges(faces), start_points, 'k'); hold off; axis tight; axis image; axis off; colormap jet(256);

Examples of triangulations. Notice how the canonical euclidean Delaunay triangulation (left) differs from the geodesic one (right) when the metric is not constant.

Farthest point sampling.

• We are now back to computations in 2D. The function farthest_point_sampling iteratively compute the distance to the already computed set of points, and add to the list (initialized as empty) the farthest point. You should read the function so that you fully understand what it is doing. You can also try different speed functions to see the resulting sampling.

n=300; name = 'mountain'; W = load_potential_map(name, n); % perform sampling landmark = []; landmark = farthest_point_sampling( W, landmark, nbr_landmarks-size(landmark,2) ); % display hold on; imagesc(M'); plot(landmark(1,:), landmark(2,:), 'b.', 'MarkerSize', 20); hold off; axis tight; axis image; axis off; colormap gray(256); % try with other metric W, like 'constant' ...

Farthest point sampling with 50, 100, 150, 200, 250 and 300 points respectively.

• Now you can compute the corresponding triangulation using the already given code.

Farthest point triangulations.

Heuristically driven front propagation.

• The function perform_fast_marching_2d can be used with a fourth argument that gives an estimation of the distance to the end point. It helps the algorithm to reduce the number of explored pixels. This remaining distance H should be estimated quickly (hence the name "heuristic"). Here we propose to cheat and use directly the true remaining distance using a classical propagation. You have to test this code with various values of weight.

[start_points,end_points] = pick_start_end_point(W); % compute the heuristic [H,S] = perform_fast_marching_2d(W, end_points); % perform the propagation options.end_points = end_points; weight = 0.5; % should be between 0 and 1. [D,S] = perform_fast_marching_2d(W, start_points, options, H*weight); % compute the path p = extract_path_2d(D,end_points, options); % display clf; plot_fast_marching_2d(W,S,p,start_points,end_points); colormap jet(256); saveas(gcf, [rep name '-heuristic' num2str(weight) '.jpg'], 'jpg');

• As a final question, try to devise a fast way to estimate the remaining distance function. If you have no idea, you can use the function perform_fmstar_2d which implement two methods to achieve this.

Heuristic front propagation with a weight parameter of 0, 0.5 and 0.9 respectively. Notice how the explored area shrinks around the true path.

Copyright © 2006 Gabriel Peyré

Lecture 3 - Geodesic Computations
on 3D Meshes

Abstract : The goals of this lecture is to manipulate the fast marching algorithm on triangulated meshes. Application to geodesic extraction, Voronoi segmentation, remeshing and bending invariants are presented.

Setting up Matlab.

• First download the Matlab toolbox toolbox_fast_marching.zip and toolbox_graph.zip. Unzip it into your working directory. You should have directories toolbox_fast_marching/ and toolbox_graph/ in your path.
• The first thing to do is to install these toolboxes in your path.

path(path, 'toolbox_fast_marching/'); path(path, 'toolbox_fast_marching/toolbox/'); path(path, 'toolbox_graph/'); path(path, 'toolbox_graph/off/');

• Recompile the mex files for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_fast_marching/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

cd toolbox_fast_marching compile_mex; cd ..

Distance Computation on Triangulated Surface.

• Using the fast marching on a triangulated surface, one can compute the distance from a set of input points. This function also returns the segmentation of the surface into geodesic Voronoi cells.

name = 'elephant'; % other choices includes 'skull' or 'bunny' [vertex,faces] = read_mesh([name '.off']); nverts = max(size(vertex)); % number of vertices nstart = 8; % number of starting points start_points = floor(rand(nstart,1)*nverts)+1; options.end_points = end_points(:); % perform the front propagation [D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points, options); % display the distance function col = D; col(col==Inf) = 0; % the color of the vertices clf; hold on; options.face_vertex_color = col; plot_mesh(vertex, faces, options); h = plot3(vertex(1,start_points),vertex(2,start_points), vertex(3,start_points), 'r.'); set(h, 'MarkerSize', 25); hold off; colormap jet(256); axis tight; shading interp; % you can now display the segmentation function Q ...

Example of distance function from a set of random starting points (upper row), and the corresponding Voronoi segmentation (bottom row).

• You can even stop the propagation after a fixed number of steps and see the resulting partially propagated front.

options.nb_iter_max = ...; % trying with a varying number of iterations [D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points, options); % display the distance ...

Example of front propagation.

Geodesic Remeshing.

• A regular sampling of the surface can be performed by seeding in a greedy farthest fashion samples. These points can be linked according to the Voronoi cells adjacency which gives a powerful but yet simple remeshing scheme.

nbr_landmarks = ...; % number of points, eg. 400 W = ones(nverts,1); % the speed function, for now constant speed to perform uniform remeshing % perform the sampling of the surface landmark = farthest_point_sampling_mesh( vertex,faces, [], nbr_landmarks, options ); % compute the associated triangulation [D,Z,Q] = perform_fast_marching_mesh(vertex, faces, landmark); [vertex_voronoi,faces_voronoi] = compute_voronoi_triangulation_mesh(Q,vertex,faces); % display the distance function (same as before) ... % display the remeshed triangulation (same but with vertex_voronoi and faces_voronoi) ...

Farthest point sampling with an increasing number of points (upper row) and corresponding remeshing (bottom row).

• One can use a non-constant speed function in order to have an adapted remeshing. The sampling will be denser in area where the speed function is low.

% first kind of speed: low on the left, high on the right v = rescale(vertex(1,:)); options.W = rescale(v>0.5, 1, 3); options.W = options.W(:); % do the remeshing ... % second kind of speed: continuously increasing v = rescale(-vertex(1,:),1,8); options.W = v(:); % do the remeshing ...

Rows 1&2: uniform sampling and remeshing. Rows 3&4: adapted (split left/right) remeshing. Rows 5&6: adapted (continuously increasing) remeshing.

Bending Invariants of a Surface.

• One can use the Isomap procedure in order to modify the surface to get bending invariant signature, useful to perform articulation-invariant recognition of 3D surfaces. One first need to compute the pairwise geodesic distances between points on the surfaces. One then look for new 3D positions of the vertices so that the new Euclidean distances matches closely the geodesic distances. In order to speed up computation, the geodesic distances are computed only on a reduced set of landmarks points. The 3D new locations are then interpolated.

nlandmarks = 300; % number of landmarks (low to speed up) landmarks = floor(rand(nlandmarks,1)*nverts)+1; % samples landmarks at random % compute the distance between landmarks and the rest of the vertices Dland = zeros(nverts,nlandmarks); for i=1:nlandmarks fprintf('.'); [d,S,Q] = perform_fast_marching_mesh(vertex, faces, landmarks(i)); Dland(:,i) = d(:); end fprintf('\n'); % perform isomap on the reduced set of points D1 = Dland(landmarks,:); % reduced pairwise distances D1 = (D1+D1')/2; % force symmetry J = eye(nlandmarks) - ones(nlandmarks)/nlandmarks; % centering matrix K = -1/2 * J*D1*J; % inner product matrix % compute the rank-3 approximation of the inner product to compute embedding opt.disp = 0; [xy, val] = eigs(K, 3, 'LR', opt); xy = xy .* repmat(1./sqrt(diag(val))', [nlandmarks 1]);% interpolation on the full set of points % extend the embeding using geodesic interpolation vertex1 = zeros(nverts,3); deltan = mean(Dland,1); for x=1:nverts deltax = Dland(x,:); vertex1(x,:) = 1/2 * ( xy' * ( deltan-deltax )' )'; end % display both the original mesh and the embedding

Original mesh (left) and bending invariant (right).

Copyright © 2006 Gabriel Peyré

Lecture 4 - Mesh Processing

Abstract : The goal of this lecture is to manipulate a 3D mesh. This includes the loading and display of a 3D mesh and then the processing of the mesh. This processing is based on computations involving various kinds of Laplacians. These Laplacians are extensions of the classical second order derivatives to 3D meshes. They can be used to perform heat diffusion (smoothing), compression and parameterization of the mesh.

Setting up Matlab.

• First download the Matlab toolbox toolbox_graph.zip. Unzip it into your working directory. You should have a directory toolbox_graph/ in your path. Download also the set of additional meshes in .off format: meshes.zip and unzip this file into your directory.
• The first thing to do is to install this toolbox in your path.

path(path, 'toolbox_graph/'); path(path, 'toolbox_graph/off/'); path(path, 'toolbox_graph/toolbox/'); path(path, 'meshes/');

• Recompile the mex file for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available in toolbox_graph/ for MacOs and Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

cd toolbox_graph compile_mex; cd ..

Mesh Loading and Displaying.

• You can load a mesh from a file and display it. Remember that a mesh is given by two matrix: vertex store the 3D vertex positions and face store 3-tuples giving the number of the vertices forming faces.

name = 'mushroom'; % other possibilities include 'venus' %load from file [vertex,face] = read_mesh([name '.off']); % display the mesh clf; plot_mesh(vertex,face); % remove the display of the triangle shading interp; camlight;

Two examples of rendered triangle meshes.

• A mesh can also be built by starting from a coarse triangulation and then subdividing it.

name = 'sphere'; % another choice could be 'L1' j = 2; % number of subdivision levels [vertex,face] = gen_base_mesh(name, 0); % the initial mesh [vertex,face] = gen_base_mesh(name, j); clf; plot_mesh(vertex,face);

Two example of meshes computed by regular subdivision.

Heat Diffusion on 3D Meshes.

• A Laplacian on a 3D mesh is a (n,n) matrix L, where n is the number of vertices, that generalize the classical laplacian of image processing to a 3D mesh. There are several forms of laplacian depending whether it is symmetric (L'=L) and normalized (1 on the diagonal) :
  

  	L0 = D + W	(symmetric, normalized)
  	L1 = D^{-1}*L0 = Id + D^{-1}*W	(non-symmetric, non-normalized)
  	L2 = D^{-1/2}*L0*D^{-1/2} = Id + D^{-1/2}*W*D^{-1/2}	(symmetric, normalized)

  
  Where W is a weight matrix, W(i,j)=0 if (i,j) is not an edge of the graph. There are several kinds of such weights
  

  	W(i,j)=1	(combinatorial)
  	W(i,j)=1/|vi-vj|^2	(distance)
  	W(i,j)=cot(alpha_ij)+cot(beta_ij)	(harmonic)

  
  where {vi}_i are the vertex of the mesh, i.e. vi=vertex(:,i), and (alpha_ij,beta_ij) is the pair of adjacent angles to the edge (vi,vj). A gradient matrix G associated to the laplacian L is an (m,n) matrix where m is the number of edges in the mesh, that satisfies L=G'*G. It can be computed as G((i,j),k)=+sqrt(W(i,j)) if k=j and G((i,j),k)=-sqrt(W(i,j)) if k=i, and G((i,j),k)=0 otherwise.
• In the following, you can compute gradient, weights and laplacian using the compute_mesh_xxx functions.

% kind of laplacian, can be 'combinatorial', 'distance' or 'conformal' (slow) laplacian_type = ...; % load two different kind of laplacian and check the gradient factorization options.symmetrize = 1; options.normalize = 0; L0 = compute_mesh_laplacian(vertex,face,laplacian_type,options); G0 = compute_mesh_gradient(vertex,face,laplacian_type,options); disp(['Error (should be 0): ' num2str(norm(L0-G0'*G0, 'fro')) '.']); options.normalize = 1; L1 = compute_mesh_laplacian(vertex,face,laplacian_type,options); G1 = compute_mesh_gradient(vertex,face,laplacian_type,options); disp(['Error (should be 0): ' num2str(norm(L1-G1'*G1, 'fro')) '.']); % these matrices are stored as sparse matrix spy(L0);

• In order to smooth a vector f of size n (i.e. a function defined on each vertex of the mesh), one can perform a heat diffusion by solving the following PDE
  d F / dt = -L*f       with       F(x,t=0)=f
  until some stopping time t.
  When this diffusion is applied to each component of the positions of the vertices f=vertex(i,:), this smoothes the 3D mesh. Implement this PDE using an explicit discretization in time.

% the time step should be small enough dt = 0.1; % stopping time Tmax = 10; % number of steps niter = round(Tmax/dt); % initialize the 3 vectors at time t=0 vertex1 = vertex; % solve the diffusion for i=1:niter % update the position by solving the PDE vertex1 = vertex1 + ...; end

Heat diffusion on a 3D mesh, at times t=0, t=10, t=40, t=200.

• Another way to smooth a mesh is to perform the following quadratic regularization for each componnent f=vertex(i,:)
  F = argmin_g |f-g|^2 + t * |G*f|^2
  The solution of this optimization is given in closed form using the Laplacian L=G'*G as the solution of the following linear system:
  (Id+t*L)*F = f
  Solve this problem for various t on a 3D mesh. You can use the operator \ to solve a system. How does this method compares with the heat diffusion ?

% solve the equation vertex1 = ...; % display ...

Combinatorial Laplacian: Spectral Decomposition and Compression.

• The combinatorial laplacian is a linear operator (thus a NxN matrix where N is the number of vertices). It depends only on the connectivity of the mesh, thus on face only. The eigenvector of this matrix (which is symmetric, thus can be decomposed by SVD) forms an orthogonal basis of the vector space of signal of NxN values (one real value per vertex). Those functions are the extension of the Fourier oscillating functions to surfaces.

% combinatorial laplacian computation options.symmetrize = 1; options.normalize = 0; L0 = compute_mesh_laplacian(vertex,face,'combinatorial',options); %% Performing eigendecomposition [U,S,V] = svd(L0); % extract one of the eigenvectors c = U(:,end-10); % you can try with other vector with higher frequencies % assign a color to each vertex options.face_vertex_color = rescale(c, 0,255); % display clf; plot_mesh(vertex,face, options); shading interp; lighting none;

Eigenvectors of the combinatorial laplacian with increasing frequencies from left to right.

• Like the Fourier basis, the laplacian eigenvector basis can be used to perform an orthogonal decomposition of a function. In order to perform mesh compression, we decompose each coordinate X/Y/Z of the mesh into this basis. Once this decomposition has been performed, a compression is achieved by keeping only the biggest coefficients (in magnitude).

% Warning : vertex should be of size 3 x nvert keep = 5; % number of percents of coefficients kept vertex2 = (U'*vertex')'; % projection of the vector in the laplacian basis % set threshold to remove coefficients vnorm = sum(vertex2.^2, 1); vnorms = sort(vnorm); vnorms = vnorms(end:-1:1); nvert = size(vertex,2); thresh = vnorms( round(keep/100*nvert) ); % remove small coefs by thresholding vertex2 = vertex2 .* repmat( vnorm>=thresh, [3 1] ); % reconstruction vertex2 = (U*vertex2')'; % display clf; plot_mesh(vertex2,face); shading interp; camlight; axis tight;

Spectral mesh compression performed by decomposition on the eigenvectors of the laplacian.

Mesh Flattening and Parameterization.

• The eigenvector of the combinatorial laplacian can also be used to perform mesh flattening. Flattening means finding a 2D location for each vertex of the mesh. These two coordinates are composed by the eigenvectors n°2 and n°3 of the laplacian.

% load a mesh name = 'nefertiti'; [vertex,face] = read_mesh([name '.off']); A = triangulation2adjacency(face); % perform the embedding using the combinatorial eigendecomposition xy = perform_spectral_embedding(2,A); % display the flattened mesh gplot(A,xy,'k.-'); axis tight; axis square; axis off;

• This combinatorial flattening does not use geometric information since it only use the connectivity of the mesh. So any mesh with the same connectivity will have the same 2D embedding. In order to improve the quality of the embedding, one can use a conformal laplacian, who approximate the Laplace Beltrami operator of the continuous underlying surface.

% this time, we use the information from vertex to compute flattening xy = perform_spectral_embedding(2,A,vertex,face); % display gplot(A,xy,'k.-'); axis tight; axis square; axis off;

• Another way to compute the flattening is to use the Isomap algorithm. This algorithm is not based on local differential operator such as laplacian. Instead, the geodesic distance between points on the mesh graph is first computed (see the course 4 for example of geodesic computations on graphs). Then the 2D layout of point is computed in order to match the geodesic distance with the distance in the plane.

% the embedding is now computed with isomap xy = perform_spectral_embedding(2,A,vertex); % display gplot(A,xy,'k.-'); axis tight; axis square; axis off;

Comparison of the flattening obtained with the combinatorial laplacian, the conformal laplacian and isomap.

• Mesh parameterization is similar to flattening, except that we fix the boundary of the mesh to be flattened onto some given convex 2D curve. The flattening can then be proved to be 1:1 (no triangle flip) and long as the curve is convex and Id+laplacian is positive. While flattening involve spectral computation (eigenvectors extraction), which is very slow, mesh parameterization involve the solution of a sparse linear system, which is quite fast (even if the Laplacian matrix is ill-conditioned).

% pre-compute 1-ring, i.e. the neighbors of each triangle. ring = compute_1_ring(face); lap_type = 'combinatorial'; % the other choice is 'conformal' boundary_type = 'circle'; % the other choices are 'square' and 'triangle' % compute the parameterization by solving a linear system xy = compute_parametrization(vertex,face,lap_type,boundary_type,ring); % display clf; gplot(A,xy,'k.-'); axis tight; axis square; axis off;

The first row shows parameterization using the combinatorial laplacian (with various boundary conditions). This assumes that the edge lengths is 1. To take into account the geometry of the mesh, the second row uses the conformal laplacian.

Copyright © 2006 Gabriel Peyré

Lecture 5 - Graph-based
data processing

Abstract : The goal of this lecture is to manipulate data in arbitrary dimensions using graph-based method. The points is the data are linked together in a graph structure. Geodesic computations can be performed to compute distance on the dataset.

Setting up Matlab.

• First download the Matlab toolbox toolbox_dimreduc.zip and toolbox_graph.zip. Unzip it into your working directory. You should have directories toolbox_dimreduc/ and toolbox_graph/ in your path.
• The first thing to do is to install this toolbox in your path.

path(path, 'toolbox_dimreduc'); path(path, 'toolbox_dimreduc/toolbox/'); path(path, 'toolbox_dimreduc/data/'); path(path, 'toolbox_graph');

• Recompile the mex file for your machine (this can produce some warning). If it does not work, either use the already compiled mex file (they should be available for MacOs and/or Unix) or try to set up matlab with a C compiler (e.g. gcc) using 'mex -setup'.

cd toolbox_graph compile_mex; cd ..

Distance Computation on Graphs.

• You can load synthetic and real datasets (only frey_rawface is included in the distribution however)

name = 'swissroll'; % you can also try with 'scurve' n = 800; % number of points [X,col] = load_points_set( name, n ); clf; plot_scattered(X,col); axis equal; axis off;

• From points in arbitrary dimension, you can create a graph data structure using the nearest-neighbor rule (you can also try the alternative epsilon rule).

options.use_nntools = 0; % set to 1 if you have compile TStool for your machine, this will speed up computations options.nn_nbr = 7; % number of nearest neighbor % compute NN graph and the distance between nodes [D,A] = compute_nn_graph(X,options); % display the graph clf; plot_graph(A,X, col); axis tight; axis off;

• You can now compute the geodesic distance on this graph using the Dijkstra algorithm.

% set up the length of the edges (0 for no edges) D(D==Inf) = 0; % weight on the graph % find some cool location for starting point [tmp,start_point] = min( abs(col(:)-mean(col(:)))); % starting point % compute the geodesic distance [d,S] = perform_dijkstra(D, start_point); % display the graph with distance colors clf; hold on; plot_scattered(X, d); h = plot3(X(1,start_point), X(2,start_point), X(3,start_point), 'k.'); set(h, 'MarkerSize', 25); axis tight; axis off; hold off;

Two examples of 3D point clouds (left) ; corresponding NN-graph (center) ; geodesic distance to the point in black (right), blue means close.

PCA and Isomap Flattening.

• The principal component analysis realize a linear dimensionality reduction by projecting the data set on the axes of largest variance.

% dimension reduction using PCA [Y,xy] = pca(X,2); % display clf; plot_scattered(xy,col); axis equal; axis off;

• In order to better capture the manifold geometry of the data, Isomap compute the geodesic distance between pair of points, and find a low-dimensional layout that best respect these geodesic distance.

% dimension reduction using Isomap options.nn_nbr = 7; % number of NN for the graph xy = isomap(X,2, options); % display clf; plot_scattered(xy,col); axis equal; axis off;

Original 3D data set (left) ; 2D flattening using PCA (center) and using Isomap (right). Note how the PCA does not recover the simple 2D parameterization of the manifold since it is a linear process. In contrast, Isomap is able to "unfold" the curvy surface.

Dimension Reduction for Image Libraries.

• We can use the same process of flattening for data of arbitrary dimension. We first use a simple library of translating disks.

name = 'disks'; options.nbr = 1000; % Read database M = load_images_dataset(name, options); % turn it into a set of points a = size(M,1);b = size(M,2);n = size(M,3); X = reshape(M, a*b, n); % perform isomap options.nn_nbr = 7; options.use_nntools = 0; xy = isomap(X,2, options); k = 30; clf; plot_flattened_dataset(xy,M,k); % perform pca [tmp,xy] = pca(X,2); clf; plot_flattened_dataset(xy,M,k);

• We can do the same computation on a more complex library of faces.

name = 'frey_rawface'; options.nbr = 2000; % Read database M = load_images_dataset(name, options); % subsample at random n = 1000; sel = randperm(size(M,3)); M = M(:,:,sel(1:n)); M = permute(M,[2 1 3]); % fix x/y orientation of the faces % turn it into a set of points a = size(M,1);b = size(M,2);n = size(M,3); X = reshape(M, a*b, n); % perform isomap options.nn_nbr = 7; options.use_nntools = 0; xy = isomap(X,2, options); k = 30; clf; plot_flattened_dataset(xy,M,k); % perform pca [tmp,xy] = pca(X,2); clf; plot_flattened_dataset(xy,M,k);

Dimension reduction for two different libraries of images. Left: translating disks, right: face images. Although the disks data set depends on 2D translation, this is not a flat (euclidean) manifold (it is a bit curvy due to the disk shape). This is why the PCA mapping does not recover exactly a 2D square. The face database exhibits a more complex embedding.

Copyright © 2006 Gabriel Peyré

---

conclusion

0 - Basic Matlab instructions.	Matlab
1 - Active contour and level sets.	Matlab	PDF	PPT
2 - Front propagation in 2D and 3D.	Matlab	PDF	PPT
3 - Geodesic computation on 3D meshes.	Matlab	PDF	PPT
4 - Differential Calculus on 3D meshes.	Matlab	PDF	PPT
5 - High Dimensional Data Processing.	Matlab	PDF	PPT

Data Processing Using
Manifold Methods

This course is an introduction to the computational theory of manifolds. Manifold models arise in various area of mathematics, image processing, data mining or computer science. Surfaces of arbitrary dimension can be used to model non-linear datasets that one encounters in modern data processing. Numerical methods allow to exploit this geometric non-linear prior in order to extract relevant information from the data. These methods include in particular local differential computations (related to the Laplacian operator and its variants) and global distance methods (related to geodesic computations). In this course, you will learn how to perform differential and geodesic computations on images, volumes, surfaces and high dimensional graphs.

The course includes a set of Matlab experiments. These experiments give an overview of various tasks in computer vision, image processing, learning theory and mesh processing. This includes computation of shortest paths, Voronoi segmentations, geodesic Delaunay triangulations, surface flattening, dimensionality reduction and mesh processing.

One should copy/paste the provided code into a file named e.g. tp.m, and launch the script directly from Matlab command line > tp;. Some of the scripts contain "holes" that you should try to fill on your own.