postprocessing = modification of the image = graphic algorithms = image processing
Algorithms by graphic type
Five classes of image processing algorithms:
Postprocessing
Algorithms:
The image processing operators ( algorithms) can be classified into the 4 categories
Dense image[11][12][13][14][15]
"the denser the area, the more heavy the anti-aliasing have to be in order to make it look good." knighty[19]
"the details are smaller than pixel spacing, so all that remains is the bands of colour shift from period-doubling of features making it even denser" claude[20]
path finding
Algorithms[21]
Methods:
/* How to tell whether a point is to the right or left side of a line ? http://stackoverflow.com/questions/1560492/how-to-tell-whether-a-point-is-to-the-right-or-left-side-of-a-line a, b = points line = ab pont to check = z position = sign((Bx - Ax) * (Y - Ay) - (By - Ay) * (X - Ax)) It is 0 on the line, and +1 on one side, -1 on the other side.*/doubleCheckSide(doubleZx,doubleZy,doubleAx,doubleAy,doubleBx,doubleBy){return((Bx-Ax)*(Zy-Ay)-(By-Ay)*(Zx-Ax));}
/*c console programgcc t.c -Wall./a.out*/#include<stdio.h>// 3 points define triangledoubleZax=-0.250000000000000;doubleZay=0.433012701892219;// left when ydoubleZlx=-0.112538773749444;doubleZly=0.436719687479814;doubleZrx=-0.335875821657728;doubleZry=0.316782798339332;// points to test// = inside triangledoubleZx=-0.209881783739630;doubleZy=+0.4;// outside triangledoubleZxo=-0.193503885412548;doubleZyo=0.521747636163664;doubleZxo2=-0.338750000000000;doubleZyo2=+0.440690927838329;// ============ http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-2d-triangle// In general, the simplest (and quite optimal) algorithm is checking on which side of the half-plane created by the edges the point is.doublesign(doublex1,doubley1,doublex2,doubley2,doublex3,doubley3){return(x1-x3)*(y2-y3)-(x2-x3)*(y1-y3);}intPointInTriangle(doublex,doubley,doublex1,doubley1,doublex2,doubley2,doublex3,doubley3){doubleb1,b2,b3;b1=sign(x,y,x1,y1,x2,y2)<0.0;b2=sign(x,y,x2,y2,x3,y3)<0.0;b3=sign(x,y,x3,y3,x1,y1)<0.0;return((b1==b2)&&(b2==b3));}intDescribe_Position(doubleZx,doubleZy){if(PointInTriangle(Zx,Zy,Zax,Zay,Zlx,Zly,Zrx,Zry))printf(" Z is inside\n");elseprintf(" Z is outside\n");return0;}// ======================================intmain(void){Describe_Position(Zx,Zy);Describe_Position(Zxo,Zyo);Describe_Position(Zxo2,Zyo2);return0;}
Orientation and area of the triangle : how to do it ?
// gcc t.c -Wall// ./a.out#include<stdio.h>// http://ncalculators.com/geometry/triangle-area-by-3-points.htmdoubleGiveTriangleArea(doublexa,doubleya,doublexb,doubleyb,doublexc,doubleyc){return((xb*ya-xa*yb)+(xc*yb-xb*yc)+(xa*yc-xc*ya))/2.0;}/*wiki Curve_orientation[http://mathoverflow.net/questions/44096/detecting-whether-directed-cycle-is-clockwise-or-counterclockwise]The orientation of a triangle (clockwise/counterclockwise) is the sign of the determinantmatrix = { {1 , x1, y1}, {1 ,x2, y2} , {1, x3, y3}}where(x_1,y_1), (x_2,y_2), (x_3,y_3)$are the Cartesian coordinates of the three vertices of the triangle.:<math>\mathbf{O} = \begin{bmatrix}1 & x_{A} & y_{A} \\1 & x_{B} & y_{B} \\1 & x_{C} & y_{C}\end{bmatrix}.</math>A formula for its determinant may be obtained, e.g., using the method of [[cofactor expansion]]::<math>\begin{align}\det(O) &= 1\begin{vmatrix}x_{B}&y_{B}\\x_{C}&y_{C}\end{vmatrix}-x_{A}\begin{vmatrix}1&y_{B}\\1&y_{C}\end{vmatrix}+y_{A}\begin{vmatrix}1&x_{B}\\1&x_{C}\end{vmatrix} \\&= x_{B}y_{C}-y_{B}x_{C}-x_{A}y_{C}+x_{A}y_{B}+y_{A}x_{C}-y_{A}x_{B} \\&= (x_{B}y_{C}+x_{A}y_{B}+y_{A}x_{C})-(y_{A}x_{B}+y_{B}x_{C}+x_{A}y_{C}).\end{align}</math>If the determinant is negative, then the polygon is oriented clockwise. If the determinant is positive, the polygon is oriented counterclockwise. The determinant is non-zero if points A, B, and C are non-[[collinear]]. In the above example, with points ordered A, B, C, etc., the determinant is negative, and therefore the polygon is clockwise.*/doubleIsTriangleCounterclockwise(doublexa,doubleya,doublexb,doubleyb,doublexc,doubleyc){return((xb*yc+xa*yb+ya*xc)-(ya*xb+yb*xc+xa*yc));}intDescribeTriangle(doublexa,doubleya,doublexb,doubleyb,doublexc,doubleyc){doublet=IsTriangleCounterclockwise(xa,ya,xb,yb,xc,yc);doublea=GiveTriangleArea(xa,ya,xb,yb,xc,yc);if(t>0)printf("this triangle is oriented counterclockwise, determinent = %f ; area = %f\n",t,a);if(t<0)printf("this triangle is oriented clockwise, determinent = %f; area = %f\n",t,a);if(t==0)printf("this triangle is degenerate: colinear or identical points, determinent = %f; area = %f\n",t,a);return0;}intmain(){// clockwise oriented trianglesDescribeTriangle(-94,0,92,68,400,180);// https://www-sop.inria.fr/prisme/fiches/Arithmetique/index.html.enDescribeTriangle(4.0,1.0,0.0,9.0,8.0,3.0);// clockwise orientation https://people.sc.fsu.edu/~jburkardt/datasets/triangles/tex5.txt// counterclockwise oriented trianglesDescribeTriangle(-50.00,0.00,50.00,0.00,0.00,0.02);// a "cap" triangle. This example has an area of 1.DescribeTriangle(0.0,0.0,3.0,0.0,0.0,4.0);// a right triangle. This example has an area of (?? 3 ??) = 6DescribeTriangle(4.0,1.0,8.0,3.0,0.0,9.0);// https://people.sc.fsu.edu/~jburkardt/datasets/triangles/tex1.txtDescribeTriangle(-0.5,0.0,0.5,0.0,0.0,0.866025403784439);// an equilateral triangle. This triangle has an area of sqrt(3)/4.// degenerate trianglesDescribeTriangle(1.0,0.0,2.0,2.0,3.0,4.0);// This triangle is degenerate: 3 colinear points. https://people.sc.fsu.edu/~jburkardt/datasets/triangles/tex6.txtDescribeTriangle(4.0,1.0,0.0,9.0,4.0,1.0);//2 identical pointsDescribeTriangle(2.0,3.0,2.0,3.0,2.0,3.0);// 3 identical pointsreturn0;}
/*gcc p.c -Wall./a.out----------- git --------------------cd existing_foldergit initgit remote add origin git@gitlab.com:adammajewski/PointInPolygonTest_c.gitgit add .git commitgit push -u origin master*/#include<stdio.h>#define LENGTH 6/*ArgumentMeaningnvertNumber of vertices in the polygon. Whether to repeat the first vertex at the end is discussed below.vertx, vertyArrays containing the x- and y-coordinates of the polygon's vertices.testx, testyX- and y-coordinate of the test point.https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.htmlPNPOLY - Point Inclusion in Polygon TestW. Randolph Franklin (WRF)I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point,and count how many edges it crosses.At each crossing, the ray switches between inside and outside.This is called the Jordan curve theorem.The case of the ray going thru a vertex is handled correctly via a careful selection of inequalities.Don't mess with this code unless you're familiar with the idea of Simulation of Simplicity.This pretends to shift the ray infinitesimally down so that it either clearly intersects, or clearly doesn't touch.Since this is merely a conceptual, infinitesimal, shift, it never creates an intersection that didn't exist before,and never destroys an intersection that clearly existed before.The ray is tested against each edge thus:Is the point in the half-plane to the left of the extended edge? andIs the point's Y coordinate within the edge's Y-range?Handling endpoints here is tricky.I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point,and count how many edges it crosses. At each crossing,the ray switches between inside and outside. This is called the Jordan curve theorem.The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge.So basically it's keeping track of whether the number of edges crossed are even or odd.0 means even and 1 means odd.*/intpnpoly(intnvert,double*vertx,double*verty,doubletestx,doubletesty){inti,j,c=0;for(i=0,j=nvert-1;i<nvert;j=i++){if(((verty[i]>testy)!=(verty[j]>testy))&&(testx<(vertx[j]-vertx[i])*(testy-verty[i])/(verty[j]-verty[i])+vertx[i]))c=!c;}returnc;}voidCheckPoint(intnvert,double*vertx,double*verty,doubletestx,doubletesty){intflag;flag=pnpoly(nvert,vertx,verty,testx,testy);switch(flag){case0:printf("outside\n");break;case1:printf("inside\n");break;default:printf(" ???\n");}}intmain(){// values from http://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon// number from 0 to (LENGTH-1)doublezzx[LENGTH]={13.5,6.0,13.5,42.5,39.5,42.5};doublezzy[LENGTH]={100.0,70.5,41.5,56.5,69.5,84.5};CheckPoint(LENGTH,zzx,zzy,zzx[4]-0.001,zzy[4]);CheckPoint(LENGTH,zzx,zzy,zzx[4]+0.001,zzy[4]);return0;}
types:
methods
Examples
reduce the number of points, but still keep the overall shape of the curve = PolyLine Reduction
thick line
neighborhood variants of 2D algorithms:
Field line[39]
Tracing curve[40]
Methods
Task: draw a 2D curve when
trace an external ray for angle t in turns, which means ( description by Claude Heiland-Allen)[42]
this algorithm is O(period^2), which means doubling the period takes 4x as long, so it's only really feasible for periods up to a few 1000.
Three Curve Tracing Models[43]
Images
Problems:
Examples
Ray can be parametrised with radius (r)
Simple closed curve ("a connected curve that does not cross itself and ends at the same point where it begins"[46] = having no endpoints) can be parametrized with angle (t).
Sobel filter G consist of 2 kernels (masks):
The Sobel kernel contains weights for each pixel from the 8-point neighbourhood of a tested pixel. These are 3x3 kernels.
There are 2 Sobel kernels, one for computing horizontal changes and other for computing vertical changes. Notice that a large horizontal change may indicate a vertical border, and a large vertical change may indicate a horizontal border. The x-coordinate is here defined as increasing in the "right"-direction, and the y-coordinate is defined as increasing in the "down"-direction.
The Sobel kernel for computing horizontal changes is:
The Sobel kernel for computing vertical changes is:
Note that:
Pixel kernel A containing central pixel with its 3x3 neighbourhood:
Other notations for pixel kernel:
where:[54]
unsignedcharul,// upper leftunsignedcharum,// upper middleunsignedcharur,// upper rightunsignedcharml,// middle leftunsignedcharmm,// middle = central pixelunsignedcharmr,// middle rightunsignedcharll,// lower leftunsignedcharlm,// lower middleunsignedcharlr,// lower right
In array notation it is:[55]
In geographic notation usede in cellular aotomats[check spelling] it is central pixel of Moore neighbourhood.
So central (tested) pixel is:
Compute Sobel filters (where here denotes the 2-dimensional convolution operation not matrix multiplication). It is a sum of products of pixel and its weights:
Because 3 weights in each kernal are zero so there are only 6 products.[56]
shortGh=ur+2*mr+lr-ul-2*ml-ll;shortGv=ul+2*um+ur-ll-2*lm-lr;
Result is computed (magnitude of gradient):
It is a color of tested pixel.
One can also approximate result by sum of 2 magnitudes:
which is much faster to compute.[57]
Lets take array of 8-bit colors (image) called data. To find borders in this image simply do:
for(iY=1;iY<iYmax-1;++iY){for(iX=1;iX<iXmax-1;++iX){Gv=-data[iY-1][iX-1]-2*data[iY-1][iX]-data[iY-1][iX+1]+data[iY+1][iX-1]+2*data[iY+1][iX]+data[iY+1][iX+1];Gh=-data[iY+1][iX-1]+data[iY-1][iX+1]-2*data[iY][iX-1]+2*data[iY][iX+1]-data[iY-1][iX-1]+data[iY+1][iX+1];G=sqrt(Gh*Gh+Gv*Gv);if(G==0){edge[iY][iX]=255;}/* background */else{edge[iY][iX]=0;}/* boundary */}}
Note that here points on borders of array (iY= 0, iY = iYmax, iX=0, iX=iXmax) are skipped
Result is saved to another array called edge (with the same size).
One can save edge array to file showing only borders, or merge 2 arrays:
for(iY=1;iY<iYmax-1;++iY){for(iX=1;iX<iXmax-1;++iX){if(edge[iY][iX]==0)data[iY][iX]=0;}}
to have new image with marked borders.
Above example is for 8-bit or indexed color. For higher bit colors "the formula is applied to all three color channels separately" (from RoboRealm doc).
Other implementations:
Edge position:
In ImageMagick as "you can see, the edge is added only to areas with a color gradient that is more than 50% white! I don't know if this is a bug or intentional, but it means that the edge in the above is located almost completely in the white parts of the original mask image. This fact can be extremely important when making use of the results of the "-edge" operator."[58]
The result is:
See also new operators from 6 version of ImageMagick: EdgeIn and EdgeOut fromMorphology[60]
convert$tmp0-convolve"1,1,1,1,1,1,1,1,1"-threshold0$outfile
/* distance between 2 points z1 = x1 + y1*I z2 = x2 + y2*I en.wikipedia.org/wiki/Distance#Geometry*/doubleGiveDistance(intx1,inty1,intx2,inty2){returnsqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));}/*mutually and externally tangent circlesmathworld.wolfram.com/TangentCircles.htmlTwo circles are mutually and externally tangent if distance between their centers is equal to the sum of their radii*/doubleTestTangency(intx1,inty1,intr1,intx2,inty2,intr2){doubledistance;distance=GiveDistance(x1,y1,x2,y2);return(distance-(r1+r2));// return should be zero}
Stereographic projection is a map projection obtained by projecting points P on the surface of sphere from the sphere's north pole N to point P' in a plane tangent to the south pole S.[64]
take the north pole N to be the standard unit vector (0, 0, 1) and the center of the sphere to be the origin, so that the tangent plane at the south pole is has equation z = 1. Given a point P = (x, y, z) on the unit sphere which is not the north pole, its image is equal to[65]
Cylindrical projections in general have an increased vertical stretching as one moves towards either of the poles. Indeed, the poles themselves can't be represented (except at infinity). This stretching is reduced in the Mercator projection by the natural logarithm scaling.[66]
Two equivalent constructions of standard Mercator projection:
complex logarithm f(z) = ln(z). Effectively, a complex point z with polar coordinates (ρ,θ) is mapped to f(z) = ln(ρ) + iθ. This yields a horizontal strip,
the standard vertical Mercator mapping by rotating 90° afterwards, or equivalently f(z) can be defined as f(z) = i*ln(z).
| Normal Mercator | Transverse Mercator | |||
|---|---|---|---|---|
 | \ |  |
Cylindrical projections in general have an increased vertical stretching as one moves towards either of the poles. Indeed, the poles themselves can't be represented (except at infinity). This stretching is reduced in the Mercator projection by the natural logarithm scaling.[67]
Cylinder is nor finity so Mercator projections give ininite stripes with increasing distorion at poles of sphere. Thus the coordinate y of the Mercator projection becomes infinite at the poles and the map must be truncated ( cropped) at some latitude less than ninety degrees at both ends .
This need not be done symmetrically:
earth as sphere
Coordinate systems for the Earth in geography
One step method:
Consider a point on the globe of radiusR with longitudeλ and latitudeφ. The valueλ0 is the longitude of an arbitrary central meridian that is usually, but not always, that of Greenwich (i.e., zero). The anglesλ andφ are expressed in radians.
Map from point P on the unit sphere to point on the Cartesian plane
Two step method
The spherical form of the stereographic projection is usually expressed inpolar coordinates:
where is the radius of the sphere, and and are the latitude and longitude, respectively.
rectangular coordinates
Assumption: A, B, C, and D are real numbers such that.
Definitions
for every.
Cylindrical coordinates
Cylindrical coordinates (axialradiusρ,azimuthφ,elevationz) may be converted into spherical coordinates (central radiusr,inclinationθ,azimuthφ), by the formulas
Conversely, the spherical coordinates may be converted into cylindrical coordinates by the formulae
These formulae assume that the two systems have the same origin and same reference plane, measure the azimuth angleφ in the same senses from the same axis, and that the spherical angleθ is inclination from the cylindricalz axis.
Code
FilterLinear continuous-time filters
