00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00047
00048
00049
00050 #include "DiscreteApproximations.h"
00051
00052 namespace lsseg {
00053
00054
00055 void compute_squared_gradient_sum_2D(const Image<double>& image, Image<double>& result)
00056
00057 {
00058 assert(image.dimz() == 1);
00059 result.resize(image.dimx(), image.dimy(), 1, 1);
00060 result = 0;
00061 double Inc, Ipc, Icn, Icp;
00062 const int X = image.dimx();
00063 const int Y = image.dimy();
00064 int xp, xn, yp, yn;
00065
00066 for (int c = 0; c < image.numChannels(); ++c) {
00067 for (int y = 0; y < Y; ++y) {
00068 yp = (y > 0) ? y - 1 : y;
00069 yn = (y+1 < Y) ? y + 1 : y;
00070 double y_div = (yp + 2 == yn) ? 0.5 : (yp + 1 == yn) ? 1 : 0;
00071 for (int x = 0; x < X; ++x) {
00072 xp = (x > 0) ? x - 1 : x;
00073 xn = (x+1 < X) ? x + 1 : x;
00074 double x_div = (xp + 2 == xn) ? 0.5 : (xp + 1 == xn) ? 1 : 0;
00075
00076 Ipc = image(xp, y, 0, c);
00077 Inc = image(xn, y, 0, c);
00078 Icp = image( x, yp, 0, c);
00079 Icn = image( x, yn, 0, c);
00080
00081 double dx = (Inc - Ipc) * x_div;
00082 double dy = (Icn - Icp) * y_div;
00083
00084 result(x, y) += (dx * dx) + (dy * dy);
00085 }
00086 }
00087 }
00088 }
00089
00090
00091 void compute_squared_gradient_sum_3D(const Image<double>& image, Image<double>& result)
00092
00093 {
00094 result.resize(image.dimx(), image.dimy(), image.dimz(), 1);
00095 result = 0;
00096 double Incc, Ipcc, Icnc, Icpc, Iccn, Iccp;
00097 const int X = image.dimx();
00098 const int Y = image.dimy();
00099 const int Z = image.dimz();
00100 int xp, xn, yp, yn, zp, zn;
00101
00102 for (int c = 0; c < image.numChannels(); ++c) {
00103 for (int z = 0; z < Z; ++z) {
00104 zp = (z > 0) ? z - 1 : z;
00105 zn = (z+1 < Z) ? z + 1 : z;
00106 double z_div = (zp + 2 == zn) ? 0.5 : (zp + 1 == zn) ? 1 : 0;
00107 for (int y = 0; y < Y; ++y) {
00108 yp = (y > 0) ? y - 1 : y;
00109 yn = (y+1 < Y) ? y + 1 : y;
00110 double y_div = (yp + 2 == yn) ? 0.5 : (yp + 1 == yn) ? 1 : 0;
00111 for (int x = 0; x < X; ++x) {
00112 xp = (x > 0) ? x - 1 : x;
00113 xn = (x+1 < X) ? x + 1 : x;
00114 double x_div = (xp + 2 == xn) ? 0.5 : (xp + 1 == xn) ? 1 : 0;
00115
00116 Ipcc = image(xp, y, z, c);
00117 Incc = image(xn, y, z, c);
00118 Icpc = image( x, yp, z, c);
00119 Icnc = image( x, yn, z, c);
00120 Iccp = image( x, y, zp, c);
00121 Iccn = image( x, y, zn, c);
00122
00123 double dx = (Incc - Ipcc) * x_div;
00124 double dy = (Icnc - Icpc) * y_div;
00125 double dz = (Iccn - Iccp) * z_div;
00126
00127 result(x, y, z) += (dx * dx) + (dy * dy) + (dz * dz);
00128 }
00129 }
00130 }
00131 }
00132 }
00133
00134
00135
00136
00137
00138 };