colorization.rar_图形图像处理_matlab_

#include <fstream.h> #include <math.h> #include <stdlib.h> #include "tensor3d.h" Tensor3d::Tensor3d(int lx, int ly, int lz) { data = NULL ; set(lx, ly, lz) ; } Tensor3d::Tensor3d(const Tensor3d& other) { data = NULL ; *this = other ; } real Tensor3d::cube(double x, double y, double z) { int i,j; int ix = (int) x ; int iy = (int) y ; int iz = (int) z ; double dx = x - ix ; double dy = y - iy ; double dz = z - iz ; ix -- ; iy -- ; iz -- ; if(ix < 0 || iy < 0 || iz < 0 || ix+3 >= len_x || iy+3 >= len_y || iz+3 >= len_z) return 0 ; double A = dz * (dz * (-0.5 * dz + 1) - 0.5) ; double B = dz * (dz * (1.5 * dz - 2.5)) + 1 ; double C = dz * (dz * (-1.5 * dz + 2) + 0.5) ; double D = dz * (dz * (0.5 * dz - 0.5)) ; double p[16] ; for( i = 0 ; i < 4 ; i++) for( j = 0 ; j < 4 ; j++) p[j+i*4] = operator()(ix + i, iy + j, iz) * A + operator()(ix + i, iy + j, iz+1) * B + operator()(ix + i, iy + j, iz+2) * C + operator()(ix + i, iy + j, iz+3) * D ; A = dy * (dy * (-0.5 * dy + 1) - 0.5) ; B = dy * (dy * (1.5 * dy - 2.5)) + 1 ; C = dy * (dy * (-1.5 * dy + 2) + 0.5) ; D = dy * (dy * (0.5 * dy - 0.5)) ; double l[4] ; for( i = 0 ; i < 4 ; i++) l[i] = p[i*4] * A + p[i*4+1] * B + p[i*4+2] * C + p[i*4+3] * D ; A = dx * (dx * (-0.5 * dx + 1) - 0.5) ; B = dx * (dx * (1.5 * dx - 2.5)) + 1 ; C = dx * (dx * (-1.5 * dx + 2) + 0.5) ; D = dx * (dx * (0.5 * dx - 0.5)) ; return l[0] * A + l[1] * B + l[2] * C + l[3] * D ; } void Tensor3d::swap(Tensor3d& other) { int tt = other.len_x ; other.len_x = len_x ; len_x = tt ; tt = other.len_y ; other.len_y = len_y ; len_y = tt ; tt = other.plane ; other.plane = plane ; plane = tt ; tt = other.length ; other.length = length ; length = tt ; real *t = data ; data = other.data ; other.data = t ; } Tensor3d& Tensor3d::operator=(const Tensor3d& other) { len_x = other.len_x ; len_y = other.len_y ; len_z = other.len_z ; length = len_x * len_y * len_z ; plane = len_x * len_y ; if(data) delete[] data ; data = new real[length] ; if (!data) { cerr << "Allocation error of Tensor3d. (size req - " << length << ")" << endl ; exit(1) ; } for(int i = 0 ; i < length ; i++) data[i] = other.data[i] ; return *this ; } void Tensor3d::set(int lx, int ly, int lz) { if(data && lx == len_x && ly == len_y && lz == len_z) return ; len_x = lx ; len_y = ly ; len_z = lz ; length = len_x * len_y * len_z ; plane = len_x * len_y ; if(data) delete[] data ; data = new real[length] ; if (!data) { cerr << "Allocation error of Tensor3d. (size req - " << length << ")" << endl ; exit(1) ; } } void Tensor3d::clear(real val) { for(int i = 0 ; i < length ; i++) data[i] = val ; } void Tensor3d::restrict_full(Tensor3d& r) { r.set((len_x - 1) / 2 + 1, (len_y - 1) / 2 + 1, (len_z - 1) / 2 + 1) ; int i, j,x,y,z ; for( z = 1 ; z < r.len_z - 1 ; z++) { for(y = 1 ; y < r.len_y - 1 ; y++) { for( x = 1 ; x < r.len_x - 1 ; x++) { i = x + y * r.len_x + z * r.plane ; j = 2 * x + 2 * y * len_x + 2 * z * plane ; r.data[i] = 0.125 * data[j] + 0.0625 * (data[j-1] + data[j+1] + data[j-len_x] + data[j+len_x] + data[j-plane] + data[j+plane]) + 0.03125 * (data[j-1-len_x] + data[j+1-len_x] + data[j-1+len_x] + data[j+1+len_x] + data[j-1-plane] + data[j+1-plane] + data[j-1+plane] + data[j+1+plane] + data[j-len_x-plane] + data[j+len_x-plane] + data[j-len_x+plane] + data[j+len_x+plane]) + 0.015625 * (data[j-1-len_x-plane] + data[j-1-len_x+plane] + data[j+1-len_x-plane] + data[j+1-len_x+plane] + data[j-1+len_x-plane] + data[j-1+len_x+plane] + data[j+1+len_x-plane] + data[j+1+len_x+plane]) ; } } } /* Z-Plane (XY) */ for( y = 1 ; y < r.len_y - 1 ; y++) { for( x = 1 ; x < r.len_x - 1 ; x++) { i = x + y * r.len_x ; j = 2 * x + 2 * y * len_x ; r.data[i] = 0.1666667 * data[j] + 0.0833333 * (data[j-1] + data[j+1] + data[j-len_x] + data[j+len_x] + data[j+plane]) + 0.0416667 * (data[j-1+len_x] + data[j+1+len_x] + data[j-1-len_x] + data[j+1-len_x] + data[j-1+plane] + data[j+1+plane] + data[j+plane-len_x] + data[j+plane+len_x]) + 0.0208333 * (data[j-1-len_x+plane] + data[j+1-len_x+plane] + data[j-1+len_x+plane] + data[j+1+len_x+plane]) ; i += (r.len_z - 1) * r.plane ; j += (len_z - 1) * plane ; r.data[i] = 0.1666667 * data[j] + 0.0833333 * (data[j-1] + data[j+1] + data[j-len_x] + data[j+len_x] + data[j-plane]) + 0.0416667 * (data[j-1+len_x] + data[j+1+len_x] + data[j-1-len_x] + data[j+1-len_x] + data[j-1-plane] + data[j+1-plane] + data[j-plane-len_x] + data[j-plane+len_x]) + 0.0208333 * (data[j-1-len_x-plane] + data[j+1-len_x-plane] + data[j-1+len_x-plane] + data[j+1+len_x-plane]) ; } } /* Y-Plane (ZX) */ for( z = 1 ; z < r.len_z - 1 ; z++) { for( x = 1 ; x < r.len_x - 1 ; x++) { i = x + z * r.plane ; j = 2 * x + 2 * z * plane ; r.data[i] = 0.1666667 * data[j] + 0.0833333 * (data[j-1] + data[j+1] + data[j-plane] + data[j+plane] + data[j+len_x]) + 0.0416667 * (data[j-1+plane] + data[j+1+plane] + data[j-1-plane] + data[j+1-plane] + data[j-1+len_x] + data[j+1+len_x] + data[j+len_x-plane] + data[j+len_x+plane]) + 0.0208333 * (data[j-1-plane+len_x] + data[j+1-plane+len_x] + data[j-1+plane+len_x] + data[j+1+plane+len_x]) ; i += (r.len_y-1) * r.len_x ; j += (len_y-1) * len_x ; r.data[i] = 0.1666667 * data[j] + 0.0833333 * (data[j-1] + data[j+1] + data[j-plane] + data[j+plane] + data[j-len_x]) + 0.0416667 * (data[j-1+plane] + data[j+1+plane] + data[j-1-plane] + data[j+1-plane] + data[j-1-len_x] + data[j+1-len_x] + data[j-len_x-plane] + data[j-len_x+plane]) + 0.0208333 * (data[j-1-plane-len_x] + data[j+1-plane-len_x] + data[j-1+plane-len_x] + data[j+1+plane-len_x]) ; } } /* X-Plane (ZY) */ for( z = 1 ; z < r.len_z - 1 ; z++) { for( y = 1 ; y < r.len_y - 1 ; y++) { i = y * r.len_x + z * r.plane ; j = 2 * y * len_x + 2 * z * plane ; r.data[i] = 0.1666667 * data[j] + 0.0833333 * (data[j-len_x] + data[j+len_x] + data[j-plane] + data[j+plane] + data[j+1]) + 0.0416667 * (data[j-len_x+plane] + data[j+len_x+plane] + data[j-len_x-plane] + data[j+len_x-plane] + data[j-len_x+1] + data[j+len_x+1] + data[j+1-plane] + data[j+1+plane]) + 0.0208333 * (data[j-len_x-plane+1] + data[j+len_x-plane+1] + data[j-len_x+plane+1] + data[j+len_x+plane+1]) ; i += r.len_x - 1 ; j += len_x - 1 ; r.data[i] = 0.1666667 * data[j] + 0.0833333 * (data[j-len_x] + data[j+len_x] + data[j-plane] + data[j+plane] + data[j-1]) + 0.0416667 * (data[j-len_x+plane] + data[j+len_x+plane] + data[j-len_x-plane] + data[j+len_x-plane] + data[j-len_x-1] + data[j+len_x-1] + data[j-1-plane] + data[j-1+plane]) + 0.0208333 * (data[j-len_x-plane-1] + data[j+len_x-plane-1] + data[j-len_x+plane-1] + data[j+len_x+plane-1]) ; } } /* X-Axis */ for( x = 1 ; x < r.len_x - 1 ; x++) { i = x ; j = 2 * x ; r.data[i] = 0.222222 * data[j] + 0.111111 * (data[j+1] + data[j-1] + data[j+len_x] + data[j+plane]) + 0.055556 * (data[j+len_x+1] + data[j+len_x-1] + data[j+plane+1] + data[j+plane-1] + data[j+plane+len_x]) + 0.027778 * (data[j+len_x+1+plane] + data[j+len_x-1+plane]) ; i = x + (r.len_y-1) * r.len_x ; j = 2 * x + (len_y-1) * len_x ; r.data[i] = 0.2222