4. Vesselness-Filter

Vesselness-Filter

  • J

    Heyho,

    ich versuche gerade auf einem 3D-Datensatz eine Vesselness-Filterung vor zu nehmen.
    Meine Recherchen ergaben ganz schön viel Mathematik, die ich vom Grund-Ansatz her auch verstehe.
    Aber ich bin nicht in der Lage diese dann auf meine Bild-Daten an zu wenden.
    Im Anschluss an die Vesselness-Filterung möchte ich noch eien Level-Set Segmentierung durchführen.
    Ob ich nun Vesselness nach Sat et al. oder Frangi et al. mache ist mir erstmal ziemlich egal.
    Solang ich hinkriege es an zu wenden.

    Als Ursprung der Hesse-Matrix nehme ich eine Gaußsche Normalverteilung.

    f(x,y,z) = 1/(sqrt(2pi)sigma)³ * e-(x²+y²+z²)/(2sigma²)

    Die zieh ich dann nen bisl auseinander um später 1D-Filtermasken über mein Bild ziehen zu können.

    f(x,y,z) = 1/(sqrt(2pi)sigma) * e-(x²)/(2sigma²) * 1/(sqrt(2pi)sigma) * e-(y²)/(2sigma²) * 1/(sqrt(2pi)sigma) * e-(z²)/(2sigma²)

    Dann Bilde ich jeweils die erste und zweite Ableitung

    d/dx f(x,y,z) = x/(sqrt(2pi)sigma³) * e-(x²)/(2sigma²) * 1/(sqrt(2pi)sigma) * e-(y²)/(2sigma²) * 1/(sqrt(2pi)sigma) * e-(z²)/(2sigma²)

    d/dx(d/dx f(x,y,z)) = (sigma²-x²)/(sqrt(2pi)sigma5) * e-(x²)/(2sigma²) * 1/(sqrt(2pi)sigma) * e-(y²)/(2sigma²) * 1/(sqrt(2pi)sigma) * e-(z²)/(2sigma²)

    usw...

    Jetzt ziehe ich die Formeln auseinander und erhalte je Dimension eine Maske
    Also berechne ich beispielsweise Ixx meiner Hesse-Matrix wie folgt.

    #define MATH_PI 3.1415926f

VesselFilter::VesselFilter(double sigma) {
	this->sigma = sigma;
	this->maskSize = 2*sigma + 1;
	this->hessiansGenerated = false;
	this->maskGauss = new double[maskSize];
	this->maskFirstDerivativeOfGauss = new double[maskSize];
	this->maskSecDerivativeOfGauss = new double[maskSize];
}

void VesselFilter::calcMaskGauss() {
	for(int i = -maskSize/2; i < maskSize/2; i ++)
		maskGauss[i+maskSize/2] = exp(-1 * (pow(i, 2)/(2*pow(sigma,2))) ) / (sqrt(2 * MATH_PI) * sigma);
}
void VesselFilter::calcMaskFirstDerivativeOfGauss() {
	for(int i = -maskSize/2; i < maskSize/2; i ++)
		maskFirstDerivativeOfGauss[i+maskSize/2] = (i * exp(-1 * (pow(i, 2)/(2*pow(sigma,2))) )) / (sqrt(2 * MATH_PI) * pow(sigma,3));
}
void VesselFilter::calcMaskSecDerivativeOfGauss() {
	for(int i = -maskSize/2; i < maskSize/2; i ++)
		maskSecDerivativeOfGauss[i+maskSize/2] = (pow(sigma,2) - pow(i,2)) / (sqrt(2 * MATH_PI) * pow(sigma,3)) * exp(-1 * (pow(i, 2)/(2*pow(sigma,2))));
}

void VesselFilter::generadeHessianImages() {
	hessiansGenerated = true;
	gaussXX = new unsigned short **[depth];
	gaussXY = new unsigned short **[depth];
	gaussXZ = new unsigned short **[depth];
	gaussYX = gaussXY;
	gaussYY = new unsigned short **[depth];
	gaussYZ = new unsigned short **[depth];
	gaussZX = gaussXZ;
	gaussZY = gaussYZ;
	gaussZZ = new unsigned short **[depth];
	for(int z = 0; z < depth; z ++) {
		gaussXX[z] = new unsigned short *[height];
		gaussXY[z] = new unsigned short *[height];
		gaussXZ[z] = new unsigned short *[height];
		gaussYY[z] = new unsigned short *[height];
		gaussYZ[z] = new unsigned short *[height];
		gaussZZ[z] = new unsigned short *[height];
		for(int y = 0; y < height; y ++) {
			gaussXX[z][y] = new unsigned short[width];
			gaussXY[z][y] = new unsigned short[width];
			gaussXZ[z][y] = new unsigned short[width];
			gaussYY[z][y] = new unsigned short[width];
			gaussYZ[z][y] = new unsigned short[width];
			gaussZZ[z][y] = new unsigned short[width];
		}
	}
}
void VesselFilter::deleteHessianImages() {
	if(!hessiansGenerated)
		return;
	hessiansGenerated = false;

	for(int z = 0; z < depth; z ++) {
		for(int y = 0; y < height; y ++) {
			delete gaussXX[z][y];
			delete gaussXY[z][y];
			delete gaussXZ[z][y];
			delete gaussYY[z][y];
			delete gaussYZ[z][y];
			delete gaussZZ[z][y];
		}
		delete gaussXX[z];
		delete gaussXY[z];
		delete gaussXZ[z];
		delete gaussYY[z];
		delete gaussYZ[z];
		delete gaussZZ[z];
	}
	delete gaussXX;
	delete gaussXY;
	delete gaussXZ;
	delete gaussYY;
	delete gaussYZ;
	delete gaussZZ;
}

void VesselFilter::calcHessianImages3D(unsigned short ***imageData, int width, int height, int depth) {
	this->width = width;
	this->height = height;
	this->depth = depth;

	generadeHessianImages();

	unsigned short tmpXX = 0, tmpXY = 0, tmpXZ = 0, tmpYY = 0, tmpYZ = 0, tmpZZ = 0; 
	// 3D-Filter - First Step - Filter everything in X-Direktion based on input image
	for(int z = 0; z < depth; z ++) {
		for(int y = 0; y < height; y ++) {
			for(int x = maskSize/2; x < width-maskSize/2; x ++) {
				tmpXX = 0;
				tmpXY = 0;
				tmpXZ = 0;
				tmpYY = 0;
				tmpYZ = 0;
				tmpZZ = 0; 
				for(int m = -maskSize/2; m < maskSize/2; m++) {
					tmpXX += imageData[z][y][x+m] * maskSecDerivativeOfGauss[m];
					tmpXY += imageData[z][y][x+m] * maskFirstDerivativeOfGauss[m];
					tmpYY += imageData[z][y][x+m] * maskGauss[m];
				}
				tmpXZ = tmpXY;
				tmpYZ = tmpYY;
				tmpZZ = tmpYY;

				gaussXX[z][y][x] = tmpXX;
				gaussXX[z][y][x] = tmpXY;
				gaussXX[z][y][x] = tmpXZ;
				gaussXX[z][y][x] = tmpYY;
				gaussXX[z][y][x] = tmpYZ;
				gaussXX[z][y][x] = tmpZZ;
			}
		}
	}
	// 3D-Filter - Second Step - Filter everything in Y-Direktion based on X-Filtered-Image
	for(int z = 0; z < depth; z ++) {
		for(int y = maskSize/2; y < height-maskSize/2; y ++) {
			for(int x = 0; x < width; x ++) {
				tmpXX = 0;
				tmpXY = 0;
				tmpXZ = 0;
				tmpYY = 0;
				tmpYZ = 0;
				tmpZZ = 0; 
				for(int m = -maskSize/2; m < maskSize/2; m++) {
					tmpXX += gaussXX[z][y+m][x] * maskGauss[m];
					tmpXY += gaussXY[z][y+m][x] * maskFirstDerivativeOfGauss[m];
					tmpXZ += gaussXZ[z][y+m][x] * maskGauss[m];
					tmpYY += gaussYY[z][y+m][x] * maskSecDerivativeOfGauss[m];
					tmpYZ += gaussYZ[z][y+m][x] * maskFirstDerivativeOfGauss[m];
					tmpZZ += gaussXZ[z][y+m][x] * maskGauss[m];

				}
				gaussXX[z][y][x] = tmpXX;
				gaussXX[z][y][x] = tmpXY;
				gaussXX[z][y][x] = tmpXZ;
				gaussXX[z][y][x] = tmpYY;
				gaussXX[z][y][x] = tmpYZ;
				gaussXX[z][y][x] = tmpZZ;
			}
		}
	}
	// 3D-Filter - Third Step - Filter everything in Z-Direktion based on XY-Filtered-Image
	for(int z = maskSize/2; z < depth-maskSize/2; z ++) {
		for(int y = 0; y < height; y ++) {
			for(int x = 0; x < width; x ++) {
				tmpXX = 0;
				tmpXY = 0;
				tmpXZ = 0;
				tmpYY = 0;
				tmpYZ = 0;
				tmpZZ = 0; 
				for(int m = -maskSize/2; m < maskSize/2; m++) {
					tmpXX += gaussXX[z+m][y][x] * maskGauss[m];
					tmpXY += gaussXY[z+m][y][x] * maskGauss[m];
					tmpXZ += gaussXZ[z+m][y][x] * maskFirstDerivativeOfGauss[m];
					tmpYY += gaussYY[z+m][y][x] * maskGauss[m];
					tmpYZ += gaussYZ[z+m][y][x] * maskFirstDerivativeOfGauss[m];
					tmpZZ += gaussXZ[z+m][y][x] * maskSecDerivativeOfGauss[m];

				}
				gaussXX[z][y][x] = tmpXX;
				gaussXX[z][y][x] = tmpXY;
				gaussXX[z][y][x] = tmpXZ;
				gaussXX[z][y][x] = tmpYY;
				gaussXX[z][y][x] = tmpYZ;
				gaussXX[z][y][x] = tmpZZ;
			}
		}
	}
}

VesselFilter::~VesselFilter() {
	delete maskGauss;
	delete maskFirstDerivativeOfGauss;
	delete maskSecDerivativeOfGauss;
	deleteHessianImages();
}

    Nun ziehe ich meine Masken über mein 3D-Bild
    und erhalte je Maske ein "zwischenBild"

    So wie ich das verstanden habe muss ich jetzt noch die Eigenwerte und Eigenvektoren bestimmen
    Und damit dann irgendwas machen.... 😕 😕

    Aber ab hier komm ich nicht mehr weiter...

    Kann mir da jemand helfen?

    Lieben Gruß

    Jango

    0
Anmelden zum Antworten