}
void
-lmrcImageHoughTransformLines(mrcImage* out, mrcImage* in, lmrcImageHoughTransformInfo* linfo, int mode)
+lmrcImageHoughTransformInit(mrcImage* out, mrcImage* in, lmrcImageHoughTransformInfo* linfo, int mode)
{
- int iTheta, iRho, iRho0, iRho1, iRhoIndex, index, index2;
- double theta, rho;
- mrcImageParaTypeInteger x, y;
- double data, data2, weight, weight0, weight1, Rho;
-
mrcDefaultHeaderValueSet(out);
out->HeaderN.x = (int)(M_PI/linfo->deltaTheta);
out->HeaderN.y = 2*(int)(sqrt(SQR(in->HeaderN.x)+SQR(in->HeaderN.y))/linfo->deltaRho)+1;
- out->HeaderN.z = 1;
+ out->HeaderN.z = (int)((linfo->widthMax - linfo->widthMin)/linfo->deltaRho) + 1;
out->HeaderMode = mrcFloatImage;
out->HeaderLength.x = linfo->deltaTheta*DEGREE;
out->HeaderLength.y = linfo->deltaRho*in->HeaderLength.x;
out->HeaderStartN.y = -(out->HeaderN.y - 1)/2;
out->HeaderStartN.z = 0;
mrcInit(out, NULL);
+}
+
+void
+lmrcImageHoughTransformInverseInit(mrcImage* out, mrcImage* in, lmrcImageHoughTransformInfo* linfo, int mode)
+{
+ mrcDefaultHeaderValueSet(out);
+ out->HeaderN.x = linfo->Nx;
+ out->HeaderN.y = linfo->Ny;
+ out->HeaderN.z = 1;
+ out->HeaderMode = mrcFloatImage;
+ out->HeaderLength.x = linfo->dx;
+ out->HeaderLength.y = linfo->dy;
+ out->HeaderLength.z = linfo->dy;
+ out->HeaderStartN.x = 0;
+ out->HeaderStartN.y = 0;
+ out->HeaderStartN.z = 0;
+ mrcInit(out, NULL);
+}
+
+void
+lmrcImageHoughTransformLines(mrcImage* out, mrcImage* in, lmrcImageHoughTransformInfo* linfo, int mode)
+{
+ int iTheta, iRho, iRho0, iRho1, iRhoIndex, index, index2;
+ double theta, rho;
+ mrcImageParaTypeInteger x, y;
+ double data, data2, weight, weight0, weight1, Rho;
+
+
+ lmrcImageHoughTransformInit(out, in, linfo, mode);
for(y=0; y<in->HeaderN.y; y++) {
for(x=0; x<in->HeaderN.x; x++) {
mrcImageParaTypeReal x, y;
double data, data2, dx, dy, del, delta, thres;
- mrcDefaultHeaderValueSet(out);
- out->HeaderN.x = linfo->Nx;
- out->HeaderN.y = linfo->Ny;
- out->HeaderN.z = 1;
- out->HeaderMode = mrcFloatImage;
- out->HeaderLength.x = linfo->dx;
- out->HeaderLength.y = linfo->dy;
- out->HeaderLength.z = linfo->dy;
- out->HeaderStartN.x = 0;
- out->HeaderStartN.y = 0;
- out->HeaderStartN.z = 0;
- mrcInit(out, NULL);
+ lmrcImageHoughTransformInverseInit(out, in, linfo, mode);
delta = sqrt(SQR(out->HeaderLength.x) + SQR(out->HeaderLength.y))/2.0*in->HeaderLength.x*RADIAN;
thres = 1.0/sqrt(2);
mrcImageParaTypeReal x, y, srcx1, srcy1, srcx2, srcy2;
double data0, data1, data2, weight0, weight1;
- mrcDefaultHeaderValueSet(out);
- out->HeaderN.x = (int)(M_PI/linfo->deltaTheta);
- out->HeaderN.y = (int)(2*sqrt(SQR(in->HeaderN.x)+SQR(in->HeaderN.y))/linfo->deltaRho)+1;
- out->HeaderN.z = (int)((linfo->widthMax - linfo->widthMin)/linfo->deltaRho)+1;
- out->HeaderMode = mrcFloatImage;
- out->HeaderLength.x = linfo->deltaTheta*DEGREE;
- out->HeaderLength.y = linfo->deltaRho*in->HeaderLength.x;
- out->HeaderLength.z = linfo->deltaRho*in->HeaderLength.x;
- out->HeaderStartN.x = 0;
- out->HeaderStartN.y = -(out->HeaderN.y - 1)/2;
- out->HeaderStartN.z = 0;
- mrcInit(out, NULL);
+
+ lmrcImageHoughTransformInit(out, in, linfo, mode);
for(iWidth=0; iWidth<out->HeaderN.z; iWidth++) {
width = (iWidth*linfo->deltaRho + linfo->widthMin)/2.0;
mrcImageParaTypeReal x, y;
double data, data2, dx, dy, del, delta, thres;
- mrcDefaultHeaderValueSet(out);
- out->HeaderN.x = linfo->Nx;
- out->HeaderN.y = linfo->Ny;
- out->HeaderN.z = 1;
- out->HeaderMode = mrcFloatImage;
- out->HeaderLength.x = linfo->dx;
- out->HeaderLength.y = linfo->dy;
- out->HeaderLength.z = linfo->dy;
- out->HeaderStartN.x = 0;
- out->HeaderStartN.y = 0;
- out->HeaderStartN.z = 0;
- mrcInit(out, NULL);
+
+ lmrcImageHoughTransformInverseInit(out, in, linfo, mode);
thres = 1.0/sqrt(2.0);
delta = in->HeaderLength.y*out->HeaderLength.x*RADIAN;
}
}
+
+
+void
+lmrcImageHoughTransformFilledRectangles(mrcImage* out, mrcImage* in, lmrcImageHoughTransformInfo* linfo, int mode)
+{
+ int iTheta, iRho, iWidth, index, iRho0, iRho1;
+ double theta, rho, rho0, rho1, rhoCentre, Rho, width, w, widthOut;
+ mrcImageParaTypeReal x, y, srcx1, srcy1, srcx2, srcy2, srcx, srcy;
+ double data0, data1, data2, weight0, weight1, data, norm;
+
+
+ lmrcImageHoughTransformInit(out, in, linfo, mode);
+
+ for(iWidth=0; iWidth<out->HeaderN.z; iWidth++) {
+ width = (iWidth*linfo->deltaRho + linfo->widthMin)/2.0;
+ widthOut = linfo->widthOut+width;
+ norm = 1.0/out->HeaderN.x/(2.0*widthOut);
+ //DEBUGPRINT2("width: %f %f\n", width, widthOut);
+ for(y=0; y<in->HeaderN.y; y++) {
+ for(x=0; x<in->HeaderN.x; x++) {
+ for(iTheta=0; iTheta<out->HeaderN.x; iTheta++) {
+ theta = iTheta*linfo->deltaTheta;
+ rhoCentre = x*cos(theta) + y*sin(theta);
+
+ switch(mode) {
+ case 0: {
+ if(0<=rhoCentre) {
+ iRho = (int)(rhoCentre/linfo->deltaRho+0.5);
+ } else {
+ iRho = (int)(rhoCentre/linfo->deltaRho-0.5);
+ }
+ iRho -= out->HeaderStartN.y;
+
+ if(0<=iRho && iRho <out->HeaderN.y) {
+ for(w=-width; w<=width; w+=linfo->deltaRho) {
+ srcx = +w*cos(theta) + x;
+ srcy = +w*sin(theta) + y;
+ mrcPixelDataGet(in, srcx, srcy, 0, &data, mrcPixelRePart, linfo->InterpolationMode);
+ index = iTheta + iRho*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] += (data*norm);
+ }
+ for(w=-widthOut; w<-width; w+=linfo->deltaRho) {
+ srcx = +w*cos(theta) + x;
+ srcy = +w*sin(theta) + y;
+ mrcPixelDataGet(in, srcx, srcy, 0, &data, mrcPixelRePart, linfo->InterpolationMode);
+ index = iTheta + iRho*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] -= (data*norm);
+ }
+ for(w= widthOut; width<w; w-=linfo->deltaRho) {
+ srcx = +w*cos(theta) + x;
+ srcy = +w*sin(theta) + y;
+ mrcPixelDataGet(in, srcx, srcy, 0, &data, mrcPixelRePart, linfo->InterpolationMode);
+ index = iTheta + iRho*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] -= (data*norm);
+ }
+ }
+ break;
+ }
+ case 1: {
+ Rho = rhoCentre/linfo->deltaRho;
+ iRho = (int)Rho;
+ if(0<=rhoCentre) {
+ iRho0 = iRho;
+ iRho1 = iRho + 1;
+ } else {
+ iRho0 = iRho - 1;
+ iRho1 = iRho;
+ }
+ weight0 = iRho1 - Rho;
+ weight1 = Rho - iRho0;
+
+ //iRho -= out->HeaderStartN.y;
+ iRho0 -= out->HeaderStartN.y;
+ iRho1 -= out->HeaderStartN.y;
+
+ if(0<=iRho0 && iRho1 <out->HeaderN.y) {
+ for(w= -width; w<=width; w+=linfo->deltaRho) {
+ srcx = +w*cos(theta) + x;
+ srcy = +w*sin(theta) + y;
+ mrcPixelDataGet(in, srcx, srcy, 0, &data, mrcPixelRePart, linfo->InterpolationMode);
+ index = iTheta + iRho0*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] += (data*norm*weight0);
+ index = iTheta + iRho1*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] += (data*norm*weight1);
+ }
+ for(w= -widthOut; w < -width; w+=linfo->deltaRho) {
+ srcx = +w*cos(theta) + x;
+ srcy = +w*sin(theta) + y;
+ mrcPixelDataGet(in, srcx, srcy, 0, &data, mrcPixelRePart, linfo->InterpolationMode);
+ index = iTheta + iRho0*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] -= (data*norm*weight0);
+ index = iTheta + iRho1*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] -= (data*norm*weight1);
+ }
+ for(w= widthOut; width < w; w-=linfo->deltaRho) {
+ srcx = +w*cos(theta) + x;
+ srcy = +w*sin(theta) + y;
+ mrcPixelDataGet(in, srcx, srcy, 0, &data, mrcPixelRePart, linfo->InterpolationMode);
+ index = iTheta + iRho0*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] -= (data*norm*weight0);
+ index = iTheta + iRho1*out->HeaderN.x + iWidth*out->HeaderN.x*out->HeaderN.y;
+ out->ImageFloatImage[index] -= (data*norm*weight1);
+ }
+ }
+ break;
+ }
+ dafault: {
+ fprintf(stderr, "Not supported Mode: %d\n", mode);
+ exit(EXIT_FAILURE);
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+}
OSDN Git Service
