// DLL source file pascal.cpp #include "pascal.h" #include #include #include /* used for EXIT_SUCCESS */ #include LPSAFEARRAY lpsa1; LPSAFEARRAY lpsa2; #define FALSE 0 #define TRUE 1 // double Fourn(double *Data, int *NN, int NDIM, int ISIGN); double __stdcall MyFunc_atan2(double a, double b) { return atan2(a,b); } double __stdcall MyFunc_atan(double a) { return atan(a); } double __stdcall MyFunc_floor(double a) { return floor(a); } double __stdcall MyFunc_asin(double a) //dimensiot { return asin(a); } double __stdcall MyFunc_acos(double a) //dimensiot { return acos(a); } int __stdcall MyFunc_mktime(struct tm *time_s ) { int i; i = mktime(time_s); return 1; } double __stdcall MyFunc_ccor(unsigned char *akuva, short na, short ma,unsigned char *bkuva, short nb, short mb, short x_a, short y_a, short x_b, short y_b, short w, short h) // double __stdcall MyFunc_ccor(unsigned char *akuva, short na,short ma) // akuva = Bytearray size na * ma, position x_a, y_a, window height h, width w // bkuva = Bytearray size nb * mb, position x_b, y_b, window height h, width w { // char buffer[100]; unsigned char temp_a[50][50],temp_b[50][50]; // size of akuva [0...na*nb-1] int i,j; int l,m, start_row, start_col, start_index, elements, lkm_a, lkm_b; double os, nim1, nim2, meanA, meanB, cor; long jump; //j = sprintf(buffer,"Rajat %d %d %d %d \n", na, ma, nb, mb); //MessageBox ( NULL, buffer,"rajat",MB_OK); //j = sprintf(buffer,"Rajat %d %d %d %d \n", x_a, y_a, x_b, y_b); //MessageBox ( NULL, buffer,"rajat",MB_OK); //j = sprintf(buffer,"Rajat %d %d \n", w, h); //MessageBox ( NULL, buffer,"rajat",MB_OK); // fill temp_a and temp_b templates (1 to w, 1 to h) m=0; l=0; jump=0; start_row = y_a - h/2; start_col = x_a - w/2; start_index = (start_row)*(na)+start_col; elements =0; lkm_a=0; for (l=1;l<=w;l++) { jump = long(na)*long(l-1); for (m=1;m<=h;m++) { lkm_a=lkm_a+1; if ((start_index+m+jump) < 0 || (start_index+m+jump) >= (na*ma-1)) return double(elements); temp_a[l][m]=akuva[start_index+m+jump]; elements = elements + 1; } } m=0; l=0; jump=0; start_row = y_b - h/2; start_col = x_b - w/2; start_index = (start_row)*(nb)+start_col; elements =0; lkm_b=0; for (l=1;l<=w;l++) { //jump = nb-w-1; jump = long(nb)*long(l-1); for (m=1;m<=h;m++) { lkm_b=lkm_b+1; if ((start_index+m+jump) < 0 || (start_index+m+jump) >= (nb*mb-1)) return double(elements); temp_b[l][m]=bkuva[start_index+m+jump]; elements = elements + 1; } } os = 0.0; nim1 = 0.0; nim2 = 0.0; meanA = 0.0; meanB = 0.0; for (i=1;i<=w;i++) { for (j=1;j<=h;j++) { meanA = meanA + double(temp_a[i][j]); meanB = meanB + double(temp_b[i][j]); } } meanA = meanA / double(h*w); meanB = meanB / double(h*w); for (i=1;i<=w;i++) { for (j=1;j<=h;j++) { os = os + (double(temp_a[i][j]) - meanA) * (double(temp_b[i][j]) - meanB); nim1 = nim1 + (double(temp_a[i][j]) - meanA) * (double(temp_a[i][j]) - meanA); nim2 = nim2 + (double(temp_b[i][j]) - meanB) * (double(temp_b[i][j]) - meanB); } } if ((nim1 * nim2) < 0.002 ) return -2.0; cor = os / pow((nim1 * nim2),0.5); return cor; } double __stdcall MyFunc_corima(unsigned char *akuva, short width,short height, unsigned char *temp, short wt, short ht, short c_col, short c_row, double *rkuva) { /* This routine computes [-1,1] correlation image using cross-correlation. Correlation is computed for image akuva[] and template[] for the common area. BW-image BYTE (uchar) array *akuva[], declared as akuva(0 to width-1, 0 to height-1) in VB Template is the BYTE (uchar) array *temp[], declared also as temp(0 to wt-1, 0 to ht-1) in VB The center of the template (tree top's location) is given in (c_col, c_row) The output goes to the correlation double-array, declared same size as akuva-image-array. */ const int MAXDIM = 100; char buffer[MAXDIM]; unsigned char temp_a[MAXDIM][MAXDIM],temp_b[MAXDIM][MAXDIM]; double nom_term1[MAXDIM][MAXDIM]; int i,j; int l,m, startrow, startcol, startindex, lkm_a, lkm_b; double nominator, denominator1, denominator2, meanA, meanB; // j = sprintf(buffer,"Passed arguments width: %d height: %d wt: %d ht: %d \n", width, height, wt, ht); // MessageBox ( NULL, buffer,"MyFunc_corima",MB_OK); lkm_a=0; meanA=0.0; // Loop the template area, calculate mean of template, omit all zero values (outside ! for (l=0;l<=wt-1;l++) { // cols for (m=0;m<=ht-1;m++) { // rows temp_a[l][m]=temp[m*wt+l]; // minimum is zero, maximum is (ht-1)*wt+wt-1 <=> wt((ht-1)+1)-1= wt*ht-1, in VB-declaration wt*ht elemnets if (temp_a[l][m]!=0) { meanA += double(temp_a[l][m]); lkm_a++; } } } meanA /= double(lkm_a); // Too be sure we do not travel outside the image area, main loop is restricted to the area // with template width & height removed from top&bottom, left & right. // Constants denominator1=0.0; for (l=0;l<=wt-1;l++) { // loop template columns for (m=0;m<=ht-1;m++) { // loop template rows if (temp_a[l][m]!=0) // it's a valid element { denominator1 += (double(temp_a[l][m]) - meanA) * (double(temp_a[l][m]) - meanA); nom_term1[l][m]= double(temp_a[l][m]) - meanA; } } } for (i=wt;i<=(width-1)-wt;i++) { // i loops main image columns for (j=ht;j<=(height-1)-ht;j++) { // j loops main image rows // Cross-correlation is now computed for point (i,j), => rkuva[j*width+i] // c_col and c_row give the location of the HOTSPOT-point of the template. startrow = j; // -(c_row); // the location of the upper left corner of the template on the image startcol = i; // -(c_col); // " startindex = startrow*width+startcol; // " index in akuva[] lkm_b=0; meanB=0.0; for (l=0;l<=wt-1;l++) { // loop template columns for (m=0;m<=ht-1;m++) { // loop template rows if (temp_a[l][m]!=0) // it's a valid element { lkm_b++; temp_b[l][m]=akuva[startindex+(m*width)+l]; meanB += double(temp_b[l][m]); } } } meanB /= double(lkm_b); nominator=0.0; denominator2=0.0; for (l=0;l<=wt-1;l++) { // loop template cols for (m=0;m<=ht-1;m++) { // loop template rows if (temp_a[l][m]!=0) { nominator += nom_term1[l][m] * (double(temp_b[l][m]) - meanB); denominator2 += (double(temp_b[l][m]) - meanB) * (double(temp_b[l][m]) - meanB); } } } if ((denominator1 * denominator2) < 0.000002 ) { j = sprintf(buffer,"Trying to divide by zero in Myfunc_corima"); MessageBox ( NULL, buffer,"MyFunc_corima",MB_OK); return -1.0; } // Here we put the results to the location of the hotspot rkuva[(j+(c_row))*width+(i+(c_col))] = nominator/pow((denominator1 * denominator2),0.5); } // main loop rows } // main loop cols return 1.0; } double __stdcall MyFunc_Corima_Fourier(unsigned char *akuva, short na,short ma, unsigned char *temp, short wt, short ht,double *rkuva) { // akuva width na, height ma, indexing starts from [0], passed as kuva(1,1) in VBasic // temp passed temp[0] as temp(1,1), width wt, height ht // answer in double rkuva(1,1) in VB passed as rkuva[0] (C) // Compute crosscorr-image in rkuva[] na* nb // Image akuva[] na x ma // Template temp[] wt x ht // store values in FFTt[] (template) // store values in FFTima[] (image) // w_padd & H_apdd tell the numebr of zero's to padd FFTima with // double FFTt[524288] ; // max size 512x512 x 2 (complex array) /* double FFTima[long(524288)+1]; // max size 512x512 x 2 double FFTres[long(524288)+1] ; // max size 512x512 x 2 (stores the FFT'd product of temp x image* int i,j, i_re,i_im,ix,jx, lkm_c; int l,m, start_row, start_col, start_index, jump, lkm_a, lkm_b; double os, nim1, nim2, meanA, meanB; int NN[3]; int NDIM, ISIGN; char buffer[100]; j = sprintf(buffer,"Rajat %d %d %d %d \n", na, ma, wt, ht); MessageBox ( NULL, buffer,"rajat",MB_OK); for (l=1;l<=long(2)*long(na)*long(ma);l=l++) // Fill FFTt-array with template's values & zeros { FFTt[l]=double(temp[l-1]); } NN[1]=int(na); NN[2]=int(ma); // os = Fourn(&FFTt[0], &NN[0], 2, 1); for (l=1;l<=int(2)*int(na)*int(ma);l=l++) // Fill FFTt-array with template's values & zeros { FFTima[l]=double(akuva[l-1]); } // os = Fourn(&FFTima[0], &NN[0], 2, 1); for (ix = 1;ix<=wt;ix++) { for (jx = 1;jx<=ht;jx++) { // Indeces i_im = (ix * wt * 2) - wt * 2 + 2 * jx; i_re = i_im - 1; // real part: FFTres[i_re] = FFTima[i_re] * FFTt[i_re] + FFTima[i_im] * FFTt[i_im]; // imaginary part: FFTres[i_im] = FFTima[i_re] * FFTt[i_im] - FFTt[i_re] * FFTima[i_im]; } } // Calculate the inverse Fourier spectrum // os = Fourn(&FFTres[0], &NN[0], 2, -1); rkuva[0]=0.0; for (l=1;l<=long(2)*long(na)*long(ma);l=l++) // Fill FFTt-array with template's values & zeros { rkuva[l-1] = (FFTres[l] / double(double(na) * double(ma))); } // i = sprintf(buffer,"Rajat %f \n", FFTres[1]); // MessageBox ( NULL, buffer,"rajat",MB_OK); */ return 1.0; } /* double __stdcall MyFunc_Comp_Conj_Mult(double *FFTima, double *FFTt, int dummy, double *FFTres, int WCW, int HCW) { int ix, jx, i_im, i_re, i; char buffer[100]; i = sprintf(buffer,"Rajat %d %d \n", WCW, HCW); MessageBox ( NULL, buffer,"rajat",MB_OK); for (ix=1;ix<=WCW;ix++) { for (jx=1;jx<=HCW;jx++) { // real part is at index (ix * WCW * 2) - WCW * 2 + 2 * jx - 1 // imag part is at index (ix * WCW * 2) - WCW * 2 + 2 * jx i_im = (ix * WCW * 2) - WCW * 2 + 2 * jx; i_re = i_im - 1; // real part: FFTres[i_re] = FFTima[i_re] * FFTt[i_re] + FFTima[i_im] * FFTt[i_im]; // imaginary part: FFTres[i_im] = FFTima[i_re] * FFTt[i_im] - FFTt[i_re] * FFTima[i_im]; } } return 1.0; } */ double __stdcall MyFunc_Fourn(double *Data, int *NN, int NDIM, int ISIGN) { double WR, WI, WPR, WPI, WTEMP, THETA; long NTOT, IDIM, NPREV, N, NREM, IP1, IP2, IP3; long I2REV, I2, I1, I3, I3REV, IBIT; long IFP1, K1, K2, IFP2; double TEMPR, TEMPI; NTOT = 1; for (IDIM=1; IDIM<=NDIM;IDIM++) { // loop thru dimensions (2), compute total number of complex values NTOT = NTOT * NN[IDIM]; // to be stored / computed into data() complex -array, in C arrays start at [0] } NPREV = 1; for (IDIM=1;IDIM<=NDIM;IDIM++) // Main loop over the dimensions NDIM { N = NN[IDIM]; NREM = (NTOT) / (N * NPREV); IP1 = 2 * NPREV; IP2 = IP1 * N; IP3 = IP2 * NREM; I2REV = 1; for (I2 = 1;I2<=IP2; I2=I2+IP1) { // Bit reversal section of the routine if (I2 < I2REV) { for (I1=I2; I1<=(I2 + IP1 - 2);I1=I1+2) { for (I3=I1;I3<=IP3;I3=I3+IP2) { I3REV = I2REV + I3 - I2; TEMPR = Data[I3]; TEMPI = Data[I3 + 1]; Data[I3] = Data[I3REV]; Data[I3 + 1]= Data[I3REV + 1]; Data[I3REV] = TEMPR; Data[I3REV + 1] = TEMPI; } } } IBIT = IP2 / 2; line1: if ((IBIT >= IP1) && (I2REV > IBIT)) { I2REV = I2REV - IBIT; IBIT = IBIT / 2; goto line1; } I2REV = I2REV + IBIT; } // bit reversal section stops here. IFP1 = IP1; line2: if (IFP1 < IP2) { IFP2 = 2 * IFP1; THETA = ISIGN * 2.0 * PI / (IFP2 / IP1) ; // initialize for troigonometric recurrence WPR = -2.0 * sin(0.5 * THETA) * sin(0.5 * THETA) ; WPI = sin(THETA); WR = 1.0; WI = 0.0; for (I3 = 1;I3<=IFP1;I3=I3+IP1) { for (I1 = I3;I1<=(I3 + IP1 - 2);I1=I1+2) { for (I2 = I1;I2<=IP3;I2=I2+IFP2) { K1 = I2; K2 = K1 + IFP1; TEMPR = WR * Data[K2] - WI * Data[K2 + 1]; TEMPI = WR * Data[K2 + 1] + WI * Data[K2]; Data[K2] = Data[K1] - TEMPR; Data[K2 + 1] = Data[K1 + 1] - TEMPI; Data[K1] = Data[K1] + TEMPR; Data[K1 + 1] = Data[K1 + 1] + TEMPI; } } WTEMP = WR; WR = WR * WPR - WI * WPI + WR; WI = WI * WPR + WTEMP * WPI + WI; } IFP1 = IFP2; goto line2; } NPREV = N * NPREV; } return -1.0; } long __stdcall MyFunc_ComputeBkProj(int image,int Npoints, double *af, double *ori, double *cimainfo, float *atau, double *SSpace, double *cor_arr) { // Nima holds the number of images involved // Npoints holds the number of points in the 3D search volume // af[] holds the coefficients of affine transformation // ori[] holds the exterior orientation information // cimainfo[] holds the image height, and origin col,row values of the correlation image pathces // cor_arr[] holds the correlation image // FILE *f; double a_, b_, c_, d_, e_, f_; double omega, phi, kappa, Xo, Yo, Zo, c; double a[9], ck, p_x, p_y, Xa, Ya, Za, ima_x, ima_y, o_col, o_row, ima_height; int i,j, width; // char buffer[100]; // f = fopen("c:\\test.txt","w"); // Check argument flow / traffic width = (int)cimainfo[3]; ima_height = cimainfo[0]; o_col = cimainfo[1]; o_row = cimainfo[2]; // Assign values to locals and form the 3D-rotation matrix a[] omega = ori[0]; phi = ori[1]; kappa = ori[2]; Xo = ori[3]; Yo = ori[4]; Zo = ori[5]; c=ori[6]; a_=af[0];b_=af[1];c_=af[2];d_=af[3];e_=af[4];f_=af[5]; a[0] = cos(phi) * cos(kappa); a[1] = -cos(phi) * sin(kappa); a[2] = sin(phi); a[3] = cos(omega) * sin(kappa) + sin(omega) * sin(phi) * cos(kappa); a[4] = cos(omega) * cos(kappa) - sin(omega) * sin(phi) * sin(kappa); a[5] = -sin(omega) * cos(phi); a[6] = sin(omega) * sin(kappa) - cos(omega) * sin(phi) * cos(kappa); a[7] = sin(omega) * cos(kappa) + cos(omega) * sin(phi) * sin(kappa), a[8] = cos(omega) * cos(phi); int col, row; for (i = 0; i<=Npoints-1; i++) // Loop thru all points in the srach space { // Collinear equations for photo-x and photo-y j = i*3; Xa = SSpace[j] - Xo; Ya = SSpace[j+1] - Yo; Za = SSpace[j+2] - Zo; ck = -c / (a[2] * Xa + a[5] * Ya + a[8] * Za); p_x = ck * (a[0] * Xa + a[3] * Ya + a[6] * Za) ; // photo-x p_y = ck * (a[1] * Xa + a[4] * Ya + a[7] * Za) ; // photo-y ima_x = a_ * p_x + b_ * p_y + c_ ; // image-x ima_y = d_ * p_x + e_ * p_y + f_ ; // image-y // ima_x is in the main image coordinates the column-value. // ima_y is in the main image coordinates the row-value from bottom // cor_arr represents a sub-image of the main image! // o_col is in right-handed main image coords // o-row is in right-handed main image coords: Upper-left pixel's row value in the main image col = (int)(ima_x - o_col +0.5); // Col gives the column in cor_arr // ((cimainfo[0] - 1) - ima_y) == pixel coords from top // cimainfo[2] == row coordinates from top for the origin // cor_arr(row = 0) == image file structure row 0 up row = (int)(o_row - ima_y +.5); atau[i] += (float)cor_arr[9*width*row+col*9+image]; // Round-off errors!! (int) truncates!!, here we should perhaps sample a 3x3 window!, get weigths for each! } return 1; } long __stdcall MyFunc_Make3dpoints(double origoX, double origoY, double origoZ, double gridXextent, double gridYextent, double gridZextent, double gridXtess, double gridYtess, double gridZtess, double XYrotangle, char *filename, long stringlength) { FILE *f; double XYrotangle_sin; double XYrotangle_cos; double x1;double y1; double z1; double x, y, z; int l; XYrotangle_sin = sin(XYrotangle * 3.1415926 / 180); XYrotangle_cos = cos(XYrotangle * 3.1415926 / 180); f = fopen(filename,"w"); if (f==NULL) return 0; l=0; for (x=-gridXextent/2.0;x<=gridXextent/2.0;x=x+gridXtess) { for (y=-gridYextent/2.0;y<=gridYextent/2.0;y=y+gridYtess) { for (z=-gridZextent/2.0;z<=gridZextent/2.0;z=z+gridZtess) { l = l+1; x1 = XYrotangle_cos * x - XYrotangle_sin * y + origoX; y1 = XYrotangle_sin * x + XYrotangle_cos * y + origoY; z1 = origoZ+z; fprintf(f, "%.3f,%.3f,%.3f\n",x1,y1,z1); } } } fclose(f); return l; } long __stdcall MyFuncSplittoRGB(long kork, long lev, char *filename, char *filenameRed,char *filenameBlue,char *filenameGreen) { // 2.Dec.2001, Ilkka Korpela // Routine to read a RAW (rgb) -file of the size width: lev, height: kork and output // its channel split in three binary files. The file names are passed as strings (NUL-ended). struct RGBtriplet { unsigned char r; unsigned char g; unsigned char b; }; FILE *fpIn, *fpRed, *fpGreen, *fpBlue; // file pointers char buffer[100]; int i,j; RGBtriplet ImageRow[16384]; unsigned char RedRow[16384],GreenRow[16384],BlueRow[16384]; // Let the user know if the dimensions were correctly passed. i = sprintf(buffer,"Split: arvo kork %d \n", kork); MessageBox ( NULL, buffer,"Slpit",MB_OK); i = sprintf(buffer,"Split: arvo lev %d \n", lev); MessageBox ( NULL, buffer,"Slpit",MB_OK); // Open files, return if failure fpIn = fopen(filename, "rb"); if (fpIn==NULL) { i = sprintf(buffer,"Split: luettavan tiedoston avaus ei onnistunut! \n"); MessageBox ( NULL, buffer,"Slpit",MB_OK); return 0; } fpRed = fopen(filenameRed, "wb"); if (fpRed==NULL) { i = sprintf(buffer,"Split:kirjoitettavan (Red) tiedoston avaus ei onnistunut! \n"); MessageBox ( NULL, buffer,"Slpit",MB_OK); fclose(fpIn); return 0; } fpGreen = fopen(filenameGreen, "wb"); if (fpGreen==NULL) { i = sprintf(buffer,"Split:kirjoitettavan (Green) tiedoston avaus ei onnistunut! \n"); MessageBox ( NULL, buffer,"Split",MB_OK); fclose(fpIn); fclose(fpRed); return 0; } fpBlue = fopen(filenameBlue, "wb"); if (fpBlue==NULL) { i = sprintf(buffer,"Split:kirjoitettavan (Blue) tiedoston avaus ei onnistunut! \n"); MessageBox ( NULL, buffer,"Slpit",MB_OK); fclose(fpIn); fclose(fpRed); fclose(fpGreen); return 0; } // Start reading the input file row by row for (i=0; i<=int(kork)-1; i++) { fread(&ImageRow[0], sizeof(ImageRow[0])*int(lev), 1, fpIn); for (j=0; j<=int(lev)-1;j++) { RedRow[j]=ImageRow[j].r; GreenRow[j]=ImageRow[j].g; BlueRow[j]=ImageRow[j].b; } fwrite(&RedRow[0], sizeof(unsigned char)*int(lev), 1, fpRed ); fwrite(&GreenRow[0], sizeof(unsigned char)*int(lev), 1, fpGreen ); fwrite(&BlueRow[0], sizeof(unsigned char)*int(lev), 1, fpBlue ); } fclose(fpIn); fclose(fpRed); fclose(fpGreen); fclose(fpBlue); return 1; } double __stdcall MyFuncReadBinaryFiles(unsigned char *bkuva, int lstartrow, int rstartrow, int lstartcol, int rstartcol, int lendcol, int nwidth) { int i,m; char buffer[100]; unsigned char lrow[1200*2], rrow[1200*2]; char nimiL[16+1]="d:\\test_left.raw"; char nimiR[17+1]="d:\\test_right.raw"; FILE *fpInL; FILE *fpInR; MessageBox ( NULL, nimiL,"Msgbox from withn pascall.DLL",MB_OK); MessageBox ( NULL, nimiR,"Msgbox from withn pascall.DLL",MB_OK); fpInL = fopen("d:\\test_left.raw", "rb"); if (fpInL==NULL) { i = sprintf(buffer,"Normalized: Opening the RAW input file failed left, \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); return 1; } fpInR = fopen(nimiR, "rb"); if (fpInR==NULL) { i = sprintf(buffer,"Normalized: Opening the RAW input file failed right, \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpInL); return 1; } int k,j; long rl, rr, paikkal, paikkar; k = -2; for(i = 0;i<=1199;i++) { k = k + 2; rl = lstartrow + k; rr = rstartrow + k; paikkal = long(11736 - rl) * long(nwidth) + long(lstartcol) ; paikkar = long(11736 - rr) * long(nwidth) + long(rstartcol) ; fseek(fpInL, paikkal, SEEK_SET); fread(&lrow[0], sizeof(unsigned char)*1200*2, 1, fpInL); fseek(fpInR, paikkal, SEEK_SET); fread(&rrow[0], sizeof(unsigned char)*1200*2, 1, fpInR); j = -1; for (m = 0; m <= (lendcol - lstartcol);m=m+2) { j = j + 1; if (lrow[j] * 1.4 < 255) lrow[j] = lrow[j] * 1.4; bkuva[j+i*1200*3]= rrow[m]; bkuva[j+i*1200*3+1]=lrow[m]; bkuva[j+i*1200*3+2]=0; } } fclose(fpInL); fclose(fpInR); return 1; } double __stdcall MyFuncCreateNormalizedImage(double c, int Nwidth, int Nheight, int Owidth, int Oheight, char *filename_in, char *filename_out, double *pixelsize, double *af, double *R) { // af [] coeffiecients of the interrior orientation // R [] coeffiecienst of the 3x3 normalizing matrice struct RGBtriplet { unsigned char r; unsigned char g; unsigned char b; }; int i, row, col; // long maxpaikka; long paikka; double a_, b_, c_, d_, e_, f_, x_n, y_n, ima_x, ima_y, p_x, p_y; // double ec, fb, cd, af2 // double rowb, rowa, rowc; // double a1, a2, a3, a4; // double jako; FILE *fpIn; FILE *fpOut; char buffer[150]; unsigned char *norm_row; unsigned char *InImage[20000]; // double r4r8, r1r8, r2r7, r5r7, r4r6, r3r7, r0r7, r1r6, r3r8, r2r6, r0r8, r5r6; double ymax, xmin, pixelwidth; int color; int InRow, InCol; int j; unsigned int numcount; xmin=pixelsize[0]; ymax=pixelsize[1]; pixelwidth=pixelsize[2]; color=int(pixelsize[3]); if (color < 1) goto BwImage; i = _heapmin(); fpIn = fopen(filename_in, "rb"); if (fpIn==NULL) { i = sprintf(buffer,"Normalized: Opening the RAW input file failed, \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpIn); return 0; } fpOut = fopen(filename_out, "wb"); if (fpOut==NULL) { i = sprintf(buffer,"Opening the RAW output file failed\n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpOut); return 0; } norm_row = (unsigned char *) calloc(Nwidth*3, sizeof(unsigned char)); if (norm_row==NULL) { i = sprintf(buffer,"Normalized: Could not allocate memory: norm_row \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpIn); fclose(fpOut); return 0; } for (i=0;i<=Oheight-1;i++) { InImage[i] = (unsigned char *) calloc(Owidth*3,sizeof(unsigned char)); if (InImage[i]==NULL) { j = sprintf(buffer,"Normalized: Could not allocate enough memory, row: InImage %d \n",i); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpIn); fclose(fpOut); for (j=0;j<=i;j++) { free(InImage[j]); } free(norm_row); return 0; } } a_=af[0]; b_=af[1]; c_=af[2]; d_=af[3]; e_=af[4]; f_=af[5]; // Read the whole image by rows numcount = 0; paikka = 0; for (i=0;i<=Oheight-1;i++) { paikka=i*Owidth*3; fseek(fpIn, paikka, SEEK_SET); numcount += fread(&InImage[i][0], sizeof(unsigned char)*Owidth*3, 1, fpIn); } // i = sprintf(buffer,"InImage READ Bytes: %d \n", (numcount) ); // MessageBox ( NULL, buffer,"rajat",MB_OK); // Now compute OutImage for (row = 0;row<=Nheight-1;row++) { // y-coordinate in the normalized coordinates y_n = ymax - double(row)*pixelwidth + pixelwidth/2; for (col = 0;col<=Nwidth*3-1;col=col+3) { // Using 14 um pixel size, 1st row (0) y = ymin, row (height-1) y = ymax x_n = xmin + double(col/3)*pixelwidth + pixelwidth/2.0; // These are the original camera coordinates ima_x=c*(-c*R[3]*R[7]+c*R[4]*R[6]-R[3]*y_n*R[8]-x_n*R[5]*R[7]+R[4]*x_n*R[8]+y_n*R[5]*R[6])/(-c*R[3]*R[1]+c*R[4]*R[0]-R[3]*y_n*R[2]+R[4]*x_n*R[2]+y_n*R[5]*R[0]-x_n*R[5]*R[1]); ima_y=(x_n*R[2]*R[7]-x_n*R[8]*R[1]-R[6]*c*R[1]-R[2]*y_n*R[6]+c*R[0]*R[7]+R[8]*y_n*R[0])*c/(-c*R[3]*R[1]+c*R[4]*R[0]-R[3]*y_n*R[2]+R[4]*x_n*R[2]+y_n*R[5]*R[0]-x_n*R[5]*R[1]); // These are the pixel coordinates in the original image, whete to sample p_x = a_ * ima_x + b_ * ima_y + c_; p_y = d_ * ima_x + e_ * ima_y + f_; if ( (p_x < 0) || (p_x > (Owidth-1)) || (p_y < 0) || (p_y > (Oheight-1)) ) { norm_row[col] = 0; norm_row[col+1] = 0; norm_row[col+2] = 0; goto skip; } // This is the position in the file, note inverted row-direction InRow = long(Oheight-1-p_y); InCol = long(p_x)*3; norm_row[col] = InImage[InRow][InCol]; norm_row[col+1] = InImage[InRow][InCol+1]; norm_row[col+2] = InImage[InRow][InCol+2]; skip: ; } // Print a row of data in the normalized image fseek(fpOut, long(Nheight-1-row) * long(Nwidth)*3, SEEK_SET); fwrite(&norm_row[0], sizeof(unsigned char)*Nwidth*3, 1, fpOut); } free (norm_row); for (i=0;i<=Oheight-1;i++) { free (InImage[i]); } fclose(fpOut); fclose(fpIn); i = _heapmin(); return 1; // Here starts the section for BW-images (grayscale) BwImage: // i = sprintf(buffer,"InImage BWFile: %d \n", (color) ); // MessageBox ( NULL, buffer,"rajat",MB_OK); i = _heapmin(); fpIn = fopen(filename_in, "rb"); if (fpIn==NULL) { i = sprintf(buffer,"Normalized: Opening the RAW input file failed, \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpIn); return 0; } fpOut = fopen(filename_out, "wb"); if (fpOut==NULL) { i = sprintf(buffer,"Opening the RAW output file failed\n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpOut); return 0; } norm_row = (unsigned char *) calloc(Nwidth, sizeof(unsigned char)); if (norm_row==NULL) { i = sprintf(buffer,"Normalized: Could not allocate memory: norm_row \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpIn); fclose(fpOut); return 0; } for (i=0;i<=Oheight-1;i++) { InImage[i] = (unsigned char *) calloc(Owidth, sizeof(unsigned char)); if (InImage[i]==NULL) { j = sprintf(buffer,"Normalized: Could not allocate enough memory, row: InImage %d \n",i); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpIn); fclose(fpOut); for (j=0;j<=i;j++) { free(InImage[j]); } free(norm_row); return 0; } } a_=af[0]; b_=af[1]; c_=af[2]; d_=af[3]; e_=af[4]; f_=af[5]; // Read the whole image by rows numcount = 0; paikka = 0; for (i=0;i<=Oheight-1;i++) { paikka=i*Owidth; fseek(fpIn, paikka, SEEK_SET); numcount += fread(&InImage[i][0], sizeof(unsigned char)*Owidth, 1, fpIn); } // Now compute OutImage for (row = 0;row<=Nheight-1;row++) { // y-coordinate in the normalized coordinates y_n = ymax - double(row)*pixelwidth + pixelwidth/2; for (col = 0;col<=Nwidth-1;col++) { // Using 14 um pixel size, 1st row (0) y = ymin, row (height-1) y = ymax x_n = xmin + double(col)*pixelwidth + pixelwidth/2.0; // These are the original camera coordinates ima_x=c*(-c*R[3]*R[7]+c*R[4]*R[6]-R[3]*y_n*R[8]-x_n*R[5]*R[7]+R[4]*x_n*R[8]+y_n*R[5]*R[6])/(-c*R[3]*R[1]+c*R[4]*R[0]-R[3]*y_n*R[2]+R[4]*x_n*R[2]+y_n*R[5]*R[0]-x_n*R[5]*R[1]); ima_y=(x_n*R[2]*R[7]-x_n*R[8]*R[1]-R[6]*c*R[1]-R[2]*y_n*R[6]+c*R[0]*R[7]+R[8]*y_n*R[0])*c/(-c*R[3]*R[1]+c*R[4]*R[0]-R[3]*y_n*R[2]+R[4]*x_n*R[2]+y_n*R[5]*R[0]-x_n*R[5]*R[1]); // These are the pixel coordinates in the original image, whete to sample p_x = a_ * ima_x + b_ * ima_y + c_; p_y = d_ * ima_x + e_ * ima_y + f_; if ( (p_x < 0) || (p_x > (Owidth-1)) || (p_y < 0) || (p_y > (Oheight-1)) ) { norm_row[col] = 0; goto skip2; } // This is the position in the file, note inverted row-direction InRow = long(Oheight-1-p_y); InCol = long(p_x); norm_row[col] = InImage[InRow][InCol]; skip2: ; } // Print a row of data in the normalized image fseek(fpOut, long(Nheight-1-row) * long(Nwidth), SEEK_SET); fwrite(&norm_row[0], sizeof(unsigned char)*Nwidth, 1, fpOut); } // i = sprintf(buffer,"Now freeing memory: %d \n", (numcount) ); // MessageBox ( NULL, buffer,"rajat",MB_OK); free (norm_row); for (i=0;i<=Oheight-1;i++) { free (InImage[i]); } fclose(fpOut); fclose(fpIn); i = _heapmin(); return 1; } double __stdcall MyFuncCreateBmp(char *lfile, char *rfile, char *ofile, int index, int lstartcol, int rstartcol, int lstartrow, int rstartrow, int lendcol, int lendrow, int rendcol, int rendrow, double pan_x, double pan_y, int nwidth, int nheight, int stereoheight, int stereowidth, int colormodel) { // This routine reads a RAW-files and makes it into a BMP-file (Anaglyph) struct BITMAPFILEHEADER { short bfType; int bfSize; short bfReserved1; short bfReserved2; int bfOffBits; }; struct BITMAPINFOHEADER { int biSize; int biWidth; int biHeight; short biPlanes; short biBitCount; int biCompression; int biSizeImage; int biXPelsPerMeter; int biYPelsPerMeter; int biClrUsed; int biClrImportant; }; struct RGBtriplet { unsigned char r; unsigned char g; unsigned char b; }; struct BGRQ { unsigned char b; unsigned char g; unsigned char r; }; int i,j, w_byte_array, h_byte_array; int lisa, padding; int shift; FILE *fpInL; FILE *fpInR; FILE *fpOut; BITMAPFILEHEADER filehederi; BITMAPINFOHEADER hederi; char buffer[100]; unsigned char *p_lrow[2000]; //unsigned char *p_lrowbyte; //unsigned char *p_rrowbyte; unsigned char *p_rrow[2000]; unsigned char *p_bkuva[2000]; h_byte_array = stereoheight; w_byte_array = stereowidth; shift = 100; i = _heapmin(); lisa = (w_byte_array * 3) % 4; if (lisa == 0) padding = 0; else padding = 4 - lisa; int rl, rr, paikkal, paikkar, colorm, numcount; int k, m; if (colormodel==1) colorm=3; else colorm=1; // Allocate storage for (i=0;i<=stereoheight-1;i++) { p_bkuva[i] = (unsigned char *) calloc(stereowidth*3,sizeof(unsigned char)); if (p_bkuva[i]==NULL) { j = sprintf(buffer,"Normalized: Could not allocate memory, row: p_bkuva \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); for (j=0;j<=i;j++) { free(p_bkuva[j]); } return 0; } } fpInL = fopen(lfile, "rb"); if (fpInL==NULL) { i = sprintf(buffer,"Opening the RAW input file failed, \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); return 0; } fpInR = fopen(rfile, "rb"); if (fpInR==NULL) { i = sprintf(buffer,"Opening the RAW input file failed, \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpInL); return 0; } fpOut = fopen(ofile, "wb+"); if (fpOut==NULL) { i = sprintf(buffer,"Opening the BMP output file failed\n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpInL); fclose(fpInR); return 0; } //filehederi.bfType = 0x4D42; "BM" filehederi.bfType = 0x4D42; //2 filehederi.bfSize = (w_byte_array + padding) * h_byte_array * 3; //4 filehederi.bfReserved1 = 0; //2 filehederi.bfReserved2 = 0; //2 // Offset from the BitMapFileheader to the Bitmap Bits filehederi.bfOffBits = 54; // 55 (54 byte offset, read 55th as first data-byte) //24 // Rem akind of Dib - bitmap , header is defined here hederi.biSize = 40; // 4 bytes hederi.biWidth = w_byte_array; // 4 hederi.biHeight = -h_byte_array; // 4 // Rem By making Height negative we assume (WINDOWS programs assume) it's a BMP with Low-left corner on first scan-line!! hederi.biPlanes = 1; // 2 bytes hederi.biBitCount = 24; // 2 hederi.biCompression = 0; // 4 hederi.biSizeImage = 0; // 4 hederi.biXPelsPerMeter = 0; // 4 hederi.biYPelsPerMeter = 0; // 4 hederi.biClrUsed = 0; // 4 hederi.biClrImportant = 0; // 4 // Write the headers of the BMP-file (first 54 bytes) fseek(fpOut, 0, SEEK_SET); fwrite(&filehederi.bfType, 2, 1, fpOut); fseek(fpOut, 2, SEEK_SET); fwrite(&filehederi.bfSize, 4, 1, fpOut); fseek(fpOut, 2+4, SEEK_SET); fwrite(&filehederi.bfReserved1, 2, 1, fpOut); fseek(fpOut, 2+4+2, SEEK_SET); fwrite(&filehederi.bfReserved2, 2, 1, fpOut); fseek(fpOut, 2+4+2+2, SEEK_SET); fwrite(&filehederi.bfOffBits, 4, 1, fpOut); fseek(fpOut, 14, SEEK_SET); fwrite(&hederi.biSize, 4, 1, fpOut); fseek(fpOut, 14+4, SEEK_SET); fwrite(&hederi.biWidth, 4, 1, fpOut); fseek(fpOut, 18+4, SEEK_SET); fwrite(&hederi.biHeight, 4, 1, fpOut); fseek(fpOut, 22+4, SEEK_SET); fwrite(&hederi.biPlanes, 2, 1, fpOut); fseek(fpOut, 26+2, SEEK_SET); fwrite(&hederi.biBitCount, 2, 1, fpOut); fseek(fpOut, 28+2, SEEK_SET); fwrite(&hederi.biCompression, 4, 1, fpOut); fseek(fpOut, 30+4, SEEK_SET); fwrite(&hederi.biSizeImage, 4, 1, fpOut); fseek(fpOut, 34+4, SEEK_SET); fwrite(&hederi.biXPelsPerMeter, 4, 1, fpOut); fseek(fpOut, 38+4, SEEK_SET); fwrite(&hederi.biYPelsPerMeter, 4, 1, fpOut); fseek(fpOut, 42+4, SEEK_SET); fwrite(&hederi.biClrUsed, 4, 1, fpOut); fseek(fpOut, 46+4, SEEK_SET); fwrite(&hederi.biClrImportant, 4, 1, fpOut); numcount=0; int addcolor; if (colorm==1) addcolor=0; else addcolor=2; for (i=0;i<=stereoheight-1;i++) { p_lrow[i] = (unsigned char *) calloc(nwidth*colorm,sizeof(unsigned char)); if (p_lrow[i]==NULL) { j = sprintf(buffer,"Normalized: Could not allocate memory, row: p_lrow \n"); MessageBox ( NULL, buffer,"Msgbox from withn pascall.DLL",MB_OK); fclose(fpInL); fclose(fpInR); fclose(fpOut); return 0; } } for (i=0;i<=stereoheight-1;i++) { p_rrow[i] = (unsigned char *) calloc(nwidth*colorm,sizeof(unsigned char)); if (p_rrow[i]==NULL) { j = sprintf(buffer,"Normalized: Could not allocate memory, row: p_rrow \n"); MessageBox ( NULL, buffer,"Msgbox from within pascall.DLL",MB_OK); fclose(fpInL); fclose(fpInR); fclose(fpOut); return 0; } } int c1, c2, c3, c4; c1=(nwidth) * colorm; c2=(lstartcol) * colorm; c3=(rstartcol) * colorm; c4=sizeof(unsigned char)*nwidth*colorm; for (i = 0; i <= stereoheight-1;i++) { rl = lstartrow + i; paikkal = (nheight - rl) * c1 + c2; if ((paikkal < 0) || (paikkal > nheight * nwidth * colorm)) goto skipleft; fseek(fpInL,paikkal+1,SEEK_SET); numcount = fread(&p_lrow[i][0], c4, 1, fpInL); skipleft:; } for (i = 0; i <= stereoheight-1;i++) { rr = rstartrow + i; paikkar = (nheight - rr) * c1 + c3; // Read a row from the left image, colorm tells if grayscale or rgb if ((paikkar < 0) || (paikkar > nheight * nwidth * colorm)) goto skiprigth; fseek(fpInR,paikkar+1,SEEK_SET); numcount = fread(&p_rrow[i][0], c4, 1, fpInR); skiprigth:; } for (i = 0; i <= stereoheight-1;i++) { k = -colorm; for (m = 0; m <= (lendcol - lstartcol -1); m++) { k = k + colorm; // Compute average if colormodel is RGB else take graylevels switch (colormodel) { case 0: p_bkuva[i][m*3] = 0 ; // blue-cmoponenet p_bkuva[i][m*3+1] = p_rrow[i][k + addcolor]; // green p_bkuva[i][m*3+2] = p_lrow[i][k]; // red break; case 1: p_bkuva[i][m*3] = 0 ; // blue p_bkuva[i][m*3+1] = (p_rrow[i][k]+p_rrow[i][k+1]+p_rrow[i][k+2])/3; // green p_bkuva[i][m*3+2] = (p_lrow[i][k]+p_lrow[i][k+1]+p_lrow[i][k+2])/3; // red break; } } fseek(fpOut,54+i*stereowidth*3,SEEK_SET); fwrite(&p_bkuva[i][0], 3 * stereowidth, 1, fpOut); } // Write a row in the BMP-file // We're done, close files and free allocated storage & return. fclose(fpInL); fclose(fpInR); fclose(fpOut); for (j=0;j<=stereoheight-1;j++) { free(p_bkuva[j]); free(p_lrow[j]); free(p_rrow[j]); } return 1; // End of function } long __stdcall MyFuncReadTiffHeader(long *kork, long *lev, long *offset, long *konfiguraatio,int *kompressio, int *kanavia,char *filename, long *osoitteet ) { typedef unsigned char BYTE; typedef unsigned short WORD; typedef unsigned long DWORD; typedef struct Tag { WORD nNumber; WORD nType; DWORD dwLen; DWORD dwData; } TAG; typedef struct Hdr { BYTE cProcessor[2]; WORD nVersion; DWORD dwIFD; } HDR; typedef struct Img { DWORD ImageWidth; DWORD ImageLength; WORD Compression; DWORD StripOffset; WORD SamplesPerPixel; WORD PlanarConfig; } IMG; short int i,j,nTags, nNumTags, NTiles; DWORD dwIFD; HDR Hdr; IMG Img; FILE *fpIn; TAG Tag; char buffer[100]; fpIn = NULL; long TileWidth,RowsPerStrip,StripByteCounts,TileHeight,TileOffsets,TileByteCounts; short BitsPerSample, PhotometricInterpretation; i = sprintf(buffer,"Harmiton ilmoitus Dll:n sisältä\n"); MessageBox ( NULL, buffer,"Tageja",MB_OK); fpIn = fopen(filename, "rb"); if (fpIn==NULL) return 0; /* Initialize and read header from input file */ memset(&Hdr, 0, sizeof(Hdr)); fread(&Hdr, sizeof(Hdr), 1, fpIn); /* Read Offset to 0'th IFD*/ dwIFD = Hdr.dwIFD; i = sprintf(buffer,"Rajat %d \n", dwIFD); MessageBox ( NULL, buffer,"Tageja",MB_OK); /* Seek IFD address and read number of tags */ fseek(fpIn, dwIFD, SEEK_SET); fread(&nNumTags, sizeof(nNumTags), 1, fpIn); /* Read Tags */ i = sprintf(buffer,"Rajat %d \n", nNumTags); MessageBox ( NULL, buffer,"Tageja",MB_OK); nTags = nNumTags; for (j = 0; j < nTags; j++) { /* Read tag */ fread(&Tag, sizeof(TAG), 1, fpIn); switch (Tag.nNumber) { case 256: Img.ImageWidth = Tag.dwData; break; case 257: Img.ImageLength = Tag.dwData; break; case 258: BitsPerSample = (WORD)Tag.dwData;break; case 259: Img.Compression = (WORD)Tag.dwData;break; case 262: PhotometricInterpretation=(WORD)Tag.dwData;break; case 273: Img.StripOffset = Tag.dwData; break; case 277: Img.SamplesPerPixel = (WORD)Tag.dwData;break; case 278: RowsPerStrip=(WORD)Tag.dwData;break; case 279: StripByteCounts=Tag.dwData;break; case 284: Img.PlanarConfig = (WORD)Tag.dwData;break; case 322: TileWidth=Tag.dwData; i = sprintf(buffer,"Tiilitetty: TileWidth= %d \n",TileWidth); MessageBox ( NULL, buffer,"Tageja",MB_OK); break; case 323: TileHeight=Tag.dwData; i = sprintf(buffer,"Tiilitetty: TileHeight= %d \n",TileHeight); MessageBox ( NULL, buffer,"Tageja",MB_OK); break; case 324: TileOffsets=Tag.dwData; NTiles = Tag.dwLen; i = sprintf(buffer,"Tiilitetty: TileOffsets and Tiles= %d %d\n",TileOffsets,NTiles); MessageBox ( NULL, buffer,"Tageja",MB_OK); // There are Tag.dwlen tiles, the first offset (tile(0)) is at TileOffsets break; case 325: TileByteCounts=Tag.dwData; i = sprintf(buffer,"Tiilitetty: TileByteCounts= %d \n",TileByteCounts); MessageBox ( NULL, buffer,"Tageja",MB_OK); break; } i = sprintf(buffer,"Luettu Tagi %d %d %d %d \n", (WORD)Tag.nNumber, (WORD)Tag.nType ,(WORD)Tag.dwLen, (DWORD)Tag.dwData); MessageBox ( NULL, buffer,"rajat",MB_OK); } fseek(fpIn, TileOffsets, SEEK_SET); // The first tile offset is here // NOTE Here we explicitly say the number of tiles (its the value in Tag.dwLen of Tag 324!! for (i=0; i<=NTiles-1; i++) // loop thru all in order, fill osoitteet() array- (to be returned) { fread(&osoitteet[i], sizeof(osoitteet[i]), 1, fpIn); } i = sprintf(buffer,"Rajat %d %d %d %d %d %d\n", Img.ImageWidth, Img.ImageLength ,Img.Compression, Img.StripOffset,Img.SamplesPerPixel,Img.PlanarConfig); MessageBox ( NULL, buffer,"rajat",MB_OK); i = sprintf(buffer,"Rajat %d %d \n", RowsPerStrip ,PhotometricInterpretation); MessageBox ( NULL, buffer,"rajat",MB_OK); /* Check TIFF-format if (Img.Compression != 1) { fclose(fpIn); MessageBox (NULL,"Palaan. Kompressio vaara"," ",MB_OK); return(-1); } */ *lev = Img.ImageWidth; *kork = Img.ImageLength; *kanavia = Img.SamplesPerPixel; *offset=Img.StripOffset; *kompressio=Img.Compression; *konfiguraatio=Img.PlanarConfig; fcloseall(); return 1; /* OK */ } double __stdcall MyFuncMatlabCall(double b) { // This function computes x = inv(A'A)A'l char buffer[256]; FILE *fIn1; double *A, *y; int i, rows, cols; fIn1=fopen("C:\\apu\\matlab.txt","r"); i = fscanf(fIn1,"%d %d", &rows, &cols); fclose(fIn1); // i=_heapmin(); A=( double *) malloc(rows*cols* sizeof(double)); y=( double *) malloc(rows* sizeof(double)); // Inform the user about the size of the task (matrix dimensions of A) //i = sprintf(buffer,"Dim A ja y luettu: rows & cols %d %d \n", rows, cols); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); // Data in column order fIn1=fopen("C:\\apu\\A.txt","rb"); if (fIn1==NULL) { i = sprintf(buffer,"A:n avaus ei onnistunut \n"); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); } int numcount=0; fseek(fIn1,0,SEEK_SET); numcount = fread(&A[0],sizeof(double), rows*cols, fIn1); fclose(fIn1); // i = sprintf(buffer,"A:n numcount (alkioita) %d \n", numcount); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); fIn1=fopen("C:\\apu\\l.txt","rb"); if (fIn1==NULL) { i = sprintf(buffer,"l:n avaus ei onnistunut \n"); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); } numcount=0; fseek(fIn1,0,SEEK_SET); numcount += fread(&y[0],sizeof(double), rows, fIn1); // i = sprintf(buffer,"l:n numcount %d \n", numcount); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); fclose(fIn1); fIn1=fopen("C:\\apu\\test.txt","w"); for (i=0;i<=rows-1;i++) { fprintf(fIn1,"%12.8f \n", y[i]); } fclose(fIn1); // mxArray *AM=NULL; mxArray *A1=NULL; mxArray *AMT=NULL; mxArray *yM=NULL; mxArray *AMM=NULL; mxArray *AMI=NULL; mxArray *ATy=NULL; mxArray *ATAATy=NULL; mxArray *format; // String array(s) mxArray *count; // Return value mxArray *permission; // String array(s) mxArray *filename; // String array(s) mxArray *fid=NULL;; // Return value mxArray *status = NULL;/* Return value */ // mxArray *value = NULL; // mxArray *row = NULL; // mxArray *col = NULL; // mlfEnterNewContext(0,0); mxArray *AM mlfAssign(&AM,mlfDoubleMatrix(cols,rows,A,NULL)); mlfAssign(&A1,mlfDoubleMatrix(rows,cols,A,NULL)); mlfAssign(&AMT,mlfDoubleMatrix(cols,rows,A,NULL)); mlfAssign(&yM,mlfDoubleMatrix(rows,1,y,NULL)); mlfAssign(&AMM,mlfDoubleMatrix(cols,cols,A,NULL)); //i = sprintf(buffer,"Matlab taulukot muodostettu, aloitan laskennan, odota! \n"); //MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); // Transpose of A A1 = mlfTranspose(AM); // A1S = mlfSparse(A1,NULL,NULL,NULL,NULL,NULL); // A1 is the main matrix, now transpose it AMT = mlfTranspose(A1); // AMTS = mlfSparse(AMT,NULL,NULL,NULL,NULL,NULL); AMM = mlfMtimes(AMT,A1); //i = sprintf(buffer,"Kertolasku onnistui, odota! \n"); //MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); AMI = mlfInv(AMM); //i = sprintf(buffer,"Kääntö onnistui, odota! \n"); //MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ATy = mlfMtimes(AMT,yM); ATAATy = mlfMtimes(AMI,ATy); //i = sprintf(buffer,"Kertolasku kaksi onnistui, odota! \n"); //MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); mlfAssign(&format,mxCreateString("%12.8f\r")); mlfAssign(&filename,mxCreateString("c:\\apu\\x.txt")); mlfAssign(&permission,mxCreateString("w")); // Handle in fid mlfAssign(&fid, mlfFopen(NULL,NULL,filename,permission,NULL)); mlfAssign(&count, mlfFprintf(fid,format,ATAATy,NULL)); mlfAssign(&status, mlfFclose(fid)); char *apu; apu="c:\\apu\\x.txt"; i = sprintf(buffer,"Parannusvektori tiedostossa %s \n", apu); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); mxDestroyArray(AMM); mxDestroyArray(A1); mxDestroyArray(AM); mxDestroyArray(AMT); mxDestroyArray(yM); mxDestroyArray(AMI); mxDestroyArray(ATy); mxDestroyArray(ATAATy); mxDestroyArray(format); mxDestroyArray(fid); mxDestroyArray(count); mxDestroyArray(permission); mxDestroyArray(filename); mxDestroyArray(status); free (A); free (y); return 1.0; } double __stdcall MyFunc_FindClusters(long n_candidate_points, double xy_thin, double r_limit, double *plotdata, long *n_tau, float *a_tau, long *nclus, Point3D *SearchSpaceData, double *atausum, Point3D *cluscoords, float *rclus, char *filename_in, char *filename_out) { // This function does the clustering of volumetric correlation data in a_tau[] // // 20.Jan.2003 written from VB to speed up Monte Carlo -simulations // // struct TreeVect { double X; double Y; double Ztop; double Zbutt; double Zdem; double d13; double Height; int Num; int Species; int Status; double Vol; }; struct ClusMatchStruct { int Index; // Index to arrays holding clus-data unsigned char IsMatch; // True for a match with a tree unsigned char IsInside; // int MatchTreeNum; // int MatchTreeSpec; // int MatchTreeStat; double MatchTreed13; double MatchTreeh; double MatchTreeZbuttTach; double MatchTreeZbuttDem; double Xtree; // Coordinates of the tree top matched with double Ytree; double Ztree; unsigned char IsCommission; // True for a commision error (no match found for a candidate) double Dist3D; // 3D distance to the matched tree int NTreesInCylinder; // How many tree tops where in the match-cylinder int IndecesTreesInCyl[11]; double Dist2dToTreesInCyl[11]; double X; // Coordinates of the cluster double Y; double Z; double Rvalue; // 3D-correlation value int Npoints;} ; // Number of 3D-points that formed the cluster struct TreeMatchStruct { int Num; // Tree field number int Spec; // Species int Status; // Status double d13; // dbh double h; // height double X; // Treetops coordinates (using correct z-model) double Y; // double Z; // double ZButtTach; // Elevation for butt, obtained with the double ZbuttDEM; int NclusInCylinder; // Number of clusters in the match-cylinedr of the tree int IndecesClusInCyl[11]; double Dist2dToClusInCyl[11]; int IndexMatchClus; // Index of the matched cluster double Dist3D; // 3D distance betweem the top and the cluster unsigned char IsMatch; // True for a match unsigned char IsOmission; // True for a missed tree unsigned char IsInside; double Xclus; // clusters coordinates double Yclus; // double Zclus; double Vol;}; // struct PlotCenter { double X; double Y; double Z;}; struct PlotCenter Plotcenter; double plotradius, Match2Ddist, MatchZdist, ModelTreeX, ModelTreeY, ModelTreeZ; double TWidth, EllZasym, Zdiff, ZDepth, Zasym, Meshdist, Meanheight; int NimagesInMatch, ModelTreeNum; // These vars are needed for output / computations Plotcenter.X = plotdata[0]; Plotcenter.Y = plotdata[1]; Plotcenter.Z = plotdata[2]; plotradius = plotdata[3]; Match2Ddist = plotdata[4]; MatchZdist = plotdata[5]; NimagesInMatch = int(plotdata[6]); ModelTreeNum = int(plotdata[7]); ModelTreeX= plotdata[8]; ModelTreeY= plotdata[9]; ModelTreeZ= plotdata[10]; TWidth= plotdata[11]; EllZasym= plotdata[12]; Zdiff= plotdata[13]; ZDepth= plotdata[14]; Zasym= plotdata[15]; Meshdist= plotdata[16]; Meanheight= plotdata[17]; struct TreeVect Ftrees[400]; const N_CAND_MAX=2000; const N_CLUS_MAX=2000; float clus[N_CLUS_MAX+1]; int belongs,i,j,l, m, ii, jj; long mc; char buffer[200]; double xl, yl, dist, adist; FILE *fp; int ij; i = n_candidate_points; m = 0; mc = 0; // Before Euclidian distances are computed, compute BLOCK-metrics /////////////////// // ij = sprintf(buffer,"Data in\n"); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ////////////////// for (j=i; j>0 ;j--) { l = n_tau[j]; m=m+1; if (double(a_tau[l]) < r_limit) {goto piirto ;} if (m == 1) { clus[m] = 1; cluscoords[m].X = cluscoords[m].X + SearchSpaceData[l-1].X * double(a_tau[l]); cluscoords[m].Y = cluscoords[m].Y + SearchSpaceData[l-1].Y * double(a_tau[l]); cluscoords[m].Z = cluscoords[m].Z + SearchSpaceData[l-1].Z * double(a_tau[l]); atausum[m] = atausum[m] + double(a_tau[l]); mc = 1; } belongs = 0; for (ii=1; ii<=mc; ii++) { xl = cluscoords[ii].X / atausum[ii]; yl = cluscoords[ii].Y / atausum[ii]; // Do Not compute pow and sqrt when not needed if ((fabs(xl - SearchSpaceData[l-1].X) > xy_thin) || (fabs(yl - SearchSpaceData[l-1].Y) > xy_thin)) {goto NextClus;} adist = (xl - SearchSpaceData[l-1].X)*(xl - SearchSpaceData[l-1].X) + (yl - SearchSpaceData[l-1].Y)*(yl - SearchSpaceData[l-1].Y); dist = sqrt(adist); if (dist < xy_thin) { belongs= belongs + 1; if (clus[ii] > (N_CAND_MAX - 2)) { // ' overflow if > N_CAND_MAX jj = sprintf(buffer,"Virhe, liikaa kandeja! \n"); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); return -1.0;} // end if clus[ii] = clus[ii] + 1; cluscoords[ii].X = cluscoords[ii].X + SearchSpaceData[l-1].X * double(a_tau[l]); cluscoords[ii].Y = cluscoords[ii].Y + SearchSpaceData[l-1].Y * double(a_tau[l]); cluscoords[ii].Z = cluscoords[ii].Z + SearchSpaceData[l-1].Z * double(a_tau[l]); atausum[ii] = atausum[ii] + double(a_tau[l]); goto out; } // end if NextClus:; } // end for checking for existing out:; if (belongs == 0) // New cluster {mc = mc + 1; clus[mc] = 1; cluscoords[mc].X = cluscoords[mc].X + SearchSpaceData[l-1].X * double(a_tau[l]); cluscoords[mc].Y = cluscoords[mc].Y + SearchSpaceData[l-1].Y * double(a_tau[l]); cluscoords[mc].Z = cluscoords[mc].Z + SearchSpaceData[l-1].Z * double(a_tau[l]); atausum[mc] = atausum[mc] + double(a_tau[l]); } // end if // Go and pick a new point (if it still has correlation > r_limit) } // Next candidate point for loop piirto:; // Done with forming clusters. Update the label on the main -form. /////////////////// // ij = sprintf(buffer,"Clustering done\n"); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ////////////////// for (j = 1; j <= mc; j++) { rclus[j] = clus[j]; } if (mc == 0 ) { return double(mc);} // CHECK ACCURACY SECTION // // int numcount; double angle, X_origo, Y_origo, Z_origo, X_shift, Y_shift, Z_shift; fp = fopen(filename_in, "rb"); if (fp==NULL) { ij = sprintf(buffer,"BIN- Puutiedoston avaus ei onnistu mc: %d\n", mc); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); return double(mc); } int NtreesInPlot; // ij = sprintf(buffer,"Tietoja %d\n",sizeof(struct TreeVect)); // MessageBox (NULL, buffer,"MyFuncMatlabCall",MB_OK); numcount = fread(&NtreesInPlot,sizeof(int), 1, fp); /////////////////// // ij = sprintf(buffer,"Trees: %d\n", NtreesInPlot); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ////////////////// numcount = fread(&angle,sizeof(double), 1, fp); numcount = fread(&X_origo,sizeof(double), 1, fp); numcount = fread(&Y_origo,sizeof(double), 1, fp); numcount = fread(&Z_origo,sizeof(double), 1, fp); numcount = fread(&X_shift,sizeof(double), 1, fp); numcount = fread(&Y_shift,sizeof(double), 1, fp); numcount = fread(&Z_shift,sizeof(double), 1, fp); // WARNING use /Zp4 to have the size of Treevect at 76 bytes! (4 byte boundaries), use explicitly 76 bytes? numcount = fread(Ftrees,sizeof(struct TreeVect), NtreesInPlot, fp); numcount = fclose(fp); /////////////////// // ij = sprintf(buffer,"Treevect luettu \n"); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ////////////////// if (numcount!=0) { ij = sprintf(buffer,"d:\\temp\\kv\\MA3_Trees_H9.bin sulkeminen ei onnistu mc: %d\n", mc); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); return double(mc); } double MeanHeightField =0.0; double d2_dist, d3_dist; int N_tree_inside_border=0; int N_tree_inside=0; int kx, lx, jx; int N_trees_in_plot =0; struct TreeMatchStruct *TreeMatchStruct; struct ClusMatchStruct *ClusMatchStruct; TreeMatchStruct = ( struct TreeMatchStruct *) calloc(NtreesInPlot, sizeof(struct TreeMatchStruct)); ClusMatchStruct = ( struct ClusMatchStruct *) calloc(mc+1, sizeof(struct ClusMatchStruct)); /////////////////// // ij = sprintf(buffer,"struct malloc done \n"); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ////////////////// for (lx = 0; lx <= NtreesInPlot-1; lx++) // Loop thru tree-data { if ((fabs(Ftrees[lx].X - Plotcenter.X) > (plotradius + Match2Ddist)) || (fabs(Ftrees[lx].Y - Plotcenter.Y) > (plotradius + Match2Ddist))) {goto skiptree;} d2_dist = sqrt ( pow((Ftrees[lx].X - Plotcenter.X), 2) + pow((Ftrees[lx].Y - Plotcenter.Y), 2)); if (d2_dist <= (plotradius + Match2Ddist)) { N_tree_inside_border = N_tree_inside_border + 1; // N of field trees within plotradius // Running Mean kx = N_tree_inside_border ; // Number of trees inside the plot TreeMatchStruct[kx].Num = Ftrees[lx].Num; TreeMatchStruct[kx].Spec = Ftrees[lx].Species; TreeMatchStruct[kx].Status = Ftrees[lx].Status; TreeMatchStruct[kx].d13 = Ftrees[lx].d13; TreeMatchStruct[kx].h = Ftrees[lx].Height; TreeMatchStruct[kx].X = Ftrees[lx].X; TreeMatchStruct[kx].Y = Ftrees[lx].Y; TreeMatchStruct[kx].Z = Ftrees[lx].Ztop; TreeMatchStruct[kx].ZbuttDEM = Ftrees[lx].Zdem; TreeMatchStruct[kx].ZButtTach = Ftrees[lx].Zbutt; TreeMatchStruct[kx].NclusInCylinder = 0; TreeMatchStruct[kx].IndexMatchClus = 0; TreeMatchStruct[kx].Dist3D = -1; TreeMatchStruct[kx].IsMatch = FALSE; TreeMatchStruct[kx].IsOmission = FALSE; TreeMatchStruct[kx].Xclus = 0; TreeMatchStruct[kx].Yclus = 0; TreeMatchStruct[kx].Zclus = 0; TreeMatchStruct[kx].IsInside = FALSE; TreeMatchStruct[kx].Vol = Ftrees[lx].Vol; } // Check if the tree falls inside the plot + border-zone/buffer if (d2_dist < plotradius) { N_tree_inside = N_tree_inside + 1; MeanHeightField = MeanHeightField + Ftrees[lx].Ztop; TreeMatchStruct[kx].IsInside = TRUE; } skiptree:; } /////////////////// // ij = sprintf(buffer,"1. loop done \n"); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); ////////////////// N_trees_in_plot = N_tree_inside; double MeanHeightFoto,MatchXYrmse,MatchZrmse; MatchXYrmse = 0.0; MeanHeightFoto = 0.0; MatchZrmse = 0.0; double xp, yp, zp; MeanHeightField = MeanHeightField / double(N_tree_inside); int N_cand_inside=0; int N_cand_inside_border=0; for (i = 1; i <= mc; i++) { xp = cluscoords[i].X / (atausum[i]); yp = cluscoords[i].Y / (atausum[i]); zp = cluscoords[i].Z / (atausum[i]); d2_dist = sqrt ( pow((xp - Plotcenter.X), 2) + pow((yp - Plotcenter.Y),2)); if ((d2_dist > (plotradius + Match2Ddist))) {goto GetNextClus;} N_cand_inside_border = N_cand_inside_border + 1; j = N_cand_inside_border; MeanHeightFoto = MeanHeightFoto + zp; ClusMatchStruct[j].X = xp; ClusMatchStruct[j].Y = yp; ClusMatchStruct[j].Z = zp; ClusMatchStruct[j].Index = j; ClusMatchStruct[j].IsMatch = FALSE; ClusMatchStruct[j].IsCommission = TRUE; ClusMatchStruct[j].IsInside = FALSE; ClusMatchStruct[j].Npoints = int(clus[i]); ClusMatchStruct[j].Rvalue = atausum[i] / double(clus[i]); if (d2_dist < plotradius ) { N_cand_inside = N_cand_inside + 1; ClusMatchStruct[j].IsInside = TRUE; } GetNextClus:; } // Next cluster if (N_cand_inside != 0) {MeanHeightFoto = MeanHeightFoto / double(N_cand_inside);} double minXYZdist=50.0; int NTreeInClusCyl=0; // Fill ClusMatchStruct with Start-Up values for (j = 1; j <= N_cand_inside_border; j++) { minXYZdist = 50; NTreeInClusCyl = 0; for (jx = 1; jx <= N_tree_inside_border; jx++) { // if ((fabs(ClusMatchStruct[j].X - TreeMatchStruct[jx].X) > Match2Ddist) || (fabs(ClusMatchStruct[j].Y - TreeMatchStruct[jx].Y) > Match2Ddist)) {goto skipd2dist1;} d2_dist = sqrt ( pow((ClusMatchStruct[j].X - TreeMatchStruct[jx].X), 2) + pow((ClusMatchStruct[j].Y - TreeMatchStruct[jx].Y) , 2) ); // skipd2dist1:; d3_dist = sqrt ( pow((ClusMatchStruct[j].X - TreeMatchStruct[jx].X), 2) + pow((ClusMatchStruct[j].Y - TreeMatchStruct[jx].Y) , 2) + pow((ClusMatchStruct[j].Z - TreeMatchStruct[jx].Z) , 2)); if ((d2_dist < Match2Ddist) && (fabs(ClusMatchStruct[j].Z - TreeMatchStruct[jx].Z) < MatchZdist)) { NTreeInClusCyl = NTreeInClusCyl + 1; ClusMatchStruct[j].Dist2dToTreesInCyl[NTreeInClusCyl] = d3_dist; ClusMatchStruct[j].IndecesTreesInCyl[NTreeInClusCyl] = jx; TreeMatchStruct[jx].NclusInCylinder = TreeMatchStruct[jx].NclusInCylinder + 1; TreeMatchStruct[jx].IndecesClusInCyl[TreeMatchStruct[jx].NclusInCylinder] = j; TreeMatchStruct[jx].Dist2dToClusInCyl[TreeMatchStruct[jx].NclusInCylinder] = d3_dist; // m = sprintf(buffer,"In MyFunc_FindClusters %d %d %d %d %d \n", j, jx, ClusMatchStruct[j].NTreesInCylinder,TreeMatchStruct[jx].NclusInCylinder,TreeMatchStruct[jx].IndecesClusInCyl[TreeMatchStruct[jx].NclusInCylinder]); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); } if (d3_dist <= minXYZdist) { minXYZdist = d3_dist; // Rem At this point we store the minimum tree number and distance ClusMatchStruct[j].Dist3D = minXYZdist; ClusMatchStruct[j].MatchTreeNum = TreeMatchStruct[jx].Num; } // We have also use for the number of trees in a clusters match-cylinder ClusMatchStruct[j].NTreesInCylinder = NTreeInClusCyl; } // Next Tree jx } // Next cluster j int min_n; for (kx = 1; kx <= N_tree_inside_border; kx++) { if (TreeMatchStruct[kx].IsInside == TRUE) { if (TreeMatchStruct[kx].NclusInCylinder == 0) { // A SURE OMISSION ERROR TreeMatchStruct[kx].IsOmission = TRUE; goto NextTreeForCheck; } if ((TreeMatchStruct[kx].NclusInCylinder == 1) && (ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].NTreesInCylinder == 1)) { // A SURE MATCH goto markamatch; } goto NextTreeForCheck; markamatch:; TreeMatchStruct[kx].IndexMatchClus = TreeMatchStruct[kx].IndecesClusInCyl[1]; TreeMatchStruct[kx].Dist3D = TreeMatchStruct[kx].Dist2dToClusInCyl[1]; TreeMatchStruct[kx].IsMatch = TRUE; TreeMatchStruct[kx].Xclus = ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].X; TreeMatchStruct[kx].Yclus = ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].Y; TreeMatchStruct[kx].Zclus = ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].Z; // Rem Recall to mark the cluster ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].IsMatch = TRUE; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].IsCommission = FALSE; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreeNum = TreeMatchStruct[kx].Num; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreeSpec = TreeMatchStruct[kx].Spec; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreeStat = TreeMatchStruct[kx].Status; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreed13 = TreeMatchStruct[kx].d13; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreeh = TreeMatchStruct[kx].h; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreeZbuttDem = TreeMatchStruct[kx].ZbuttDEM; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].MatchTreeZbuttTach = TreeMatchStruct[kx].ZButtTach; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].Xtree = TreeMatchStruct[kx].X; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].Ytree = TreeMatchStruct[kx].Y; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].Ztree = TreeMatchStruct[kx].Z; ClusMatchStruct[TreeMatchStruct[kx].IndecesClusInCyl[1]].Dist3D = TreeMatchStruct[kx].Dist3D; } // End if Isinside NextTreeForCheck:; } // Next kx // m = sprintf(buffer,"In MyFunc_FindClusters N_match %d \n",N_match); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); long apu; int ind; for (jx = 1; jx <= N_tree_inside_border; jx++) // Loop Trees inside, allow a match from border! { if (TreeMatchStruct[jx].NclusInCylinder > 1) { // This is a tree that has several clusters in it's cylinder minXYZdist = 50; apu = 0; for (j = 1; j<= TreeMatchStruct[jx].NclusInCylinder; j++) { ind = TreeMatchStruct[jx].IndecesClusInCyl[j]; apu = apu + ClusMatchStruct[ind].NTreesInCylinder; if (TreeMatchStruct[jx].Dist2dToClusInCyl[j] < minXYZdist) { minXYZdist = TreeMatchStruct[jx].Dist2dToClusInCyl[j]; min_n = TreeMatchStruct[jx].IndecesClusInCyl[j]; } } if (apu == TreeMatchStruct[jx].NclusInCylinder) { // This tree has only clusters which on their part have just one tree, choose // the tree to attach the cluster with, it is the (min_n) -tree, make other trees missed! // There are N=apu clusters to choose from TreeMatchStruct[jx].IndexMatchClus = min_n; TreeMatchStruct[jx].Dist3D = minXYZdist; TreeMatchStruct[jx].IsMatch = TRUE; TreeMatchStruct[jx].Xclus = ClusMatchStruct[min_n].X; TreeMatchStruct[jx].Yclus = ClusMatchStruct[min_n].Y; TreeMatchStruct[jx].Zclus = ClusMatchStruct[min_n].Z; // Recall to mark the cluster ClusMatchStruct[min_n].IsMatch = TRUE; ClusMatchStruct[min_n].IsCommission = FALSE; ClusMatchStruct[min_n].MatchTreeNum = TreeMatchStruct[jx].Num; ClusMatchStruct[min_n].MatchTreeSpec = TreeMatchStruct[jx].Spec; ClusMatchStruct[min_n].MatchTreeStat = TreeMatchStruct[jx].Status; ClusMatchStruct[min_n].MatchTreed13 = TreeMatchStruct[jx].d13; ClusMatchStruct[min_n].MatchTreeh = TreeMatchStruct[jx].h; ClusMatchStruct[min_n].MatchTreeZbuttDem = TreeMatchStruct[jx].ZbuttDEM; ClusMatchStruct[min_n].MatchTreeZbuttTach = TreeMatchStruct[jx].ZButtTach; ClusMatchStruct[min_n].Xtree = TreeMatchStruct[jx].X; ClusMatchStruct[min_n].Ytree = TreeMatchStruct[jx].Y; ClusMatchStruct[min_n].Ztree = TreeMatchStruct[jx].Z; ClusMatchStruct[min_n].Dist3D = minXYZdist; goto NextMTreeToCheck; } // Check if This cluster is used // We have now NT trees with > 2 clus // if (ClusMatchStruct[min_n].IsMatch == TRUE) // { // m = sprintf(buffer,"In MyFunc_FindClusters Duplicate \n"); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); // } TreeMatchStruct[jx].IndexMatchClus = min_n; TreeMatchStruct[jx].Dist3D = minXYZdist; TreeMatchStruct[jx].IsMatch = TRUE; TreeMatchStruct[jx].Xclus = ClusMatchStruct[min_n].X; TreeMatchStruct[jx].Yclus = ClusMatchStruct[min_n].Y; TreeMatchStruct[jx].Zclus = ClusMatchStruct[min_n].Z; // Recall to mark the cluster ClusMatchStruct[min_n].IsMatch = TRUE; ClusMatchStruct[min_n].IsCommission = FALSE; ClusMatchStruct[min_n].MatchTreeNum = TreeMatchStruct[jx].Num; ClusMatchStruct[min_n].MatchTreeSpec = TreeMatchStruct[jx].Spec; ClusMatchStruct[min_n].MatchTreeStat = TreeMatchStruct[jx].Status; ClusMatchStruct[min_n].MatchTreed13 = TreeMatchStruct[jx].d13; ClusMatchStruct[min_n].MatchTreeh = TreeMatchStruct[jx].h; ClusMatchStruct[min_n].MatchTreeZbuttDem = TreeMatchStruct[jx].ZbuttDEM; ClusMatchStruct[min_n].MatchTreeZbuttTach = TreeMatchStruct[jx].ZButtTach; ClusMatchStruct[min_n].Xtree = TreeMatchStruct[jx].X; ClusMatchStruct[min_n].Ytree = TreeMatchStruct[jx].Y; ClusMatchStruct[min_n].Ztree = TreeMatchStruct[jx].Z; ClusMatchStruct[min_n].Dist3D = minXYZdist; } NextMTreeToCheck:; } // Next kx tree // Now remains to polish the Clusters-commission errors for which there lie's // a tree outside the circular plot area. These clusters must lie in the bor- // dering zone: (radius-match2ddist) - (radius) int NborderMatches = 0; for (j = 1; j <= N_cand_inside_border; j++) // loop clusters { if ((ClusMatchStruct[j].NTreesInCylinder > 0) && (ClusMatchStruct[j].IsMatch == FALSE)) { if (TreeMatchStruct[ClusMatchStruct[j].IndecesTreesInCyl[1]].NclusInCylinder > 0) { // A tree has been found for the border cluster ind = ClusMatchStruct[j].IndecesTreesInCyl[1]; if (TreeMatchStruct[ind].IsMatch == TRUE) {goto NextCand;} TreeMatchStruct[ind].IndexMatchClus = j; TreeMatchStruct[ind].Dist3D = ClusMatchStruct[j].Dist3D; TreeMatchStruct[ind].IsMatch = TRUE; TreeMatchStruct[ind].Xclus = ClusMatchStruct[j].X; TreeMatchStruct[ind].Yclus = ClusMatchStruct[j].Y; TreeMatchStruct[ind].Zclus = ClusMatchStruct[j].Z; ClusMatchStruct[j].IsMatch = TRUE; ClusMatchStruct[j].IsCommission = FALSE; ClusMatchStruct[j].MatchTreeNum = TreeMatchStruct[ind].Num; ClusMatchStruct[j].MatchTreeSpec = TreeMatchStruct[ind].Spec; ClusMatchStruct[j].MatchTreeStat = TreeMatchStruct[ind].Status; ClusMatchStruct[j].MatchTreed13 = TreeMatchStruct[ind].d13; ClusMatchStruct[j].MatchTreeh = TreeMatchStruct[ind].h; ClusMatchStruct[j].MatchTreeZbuttDem = TreeMatchStruct[ind].ZbuttDEM; ClusMatchStruct[j].MatchTreeZbuttTach = TreeMatchStruct[ind].ZButtTach; ClusMatchStruct[j].Xtree = TreeMatchStruct[ind].X; ClusMatchStruct[j].Ytree = TreeMatchStruct[ind].Y; ClusMatchStruct[j].Ztree = TreeMatchStruct[ind].Z; ClusMatchStruct[j].Dist3D = TreeMatchStruct[ind].Dist3D; NborderMatches = NborderMatches + 1; } } NextCand:; } // Next j double dXave, dYave , dZave; int CountMatches, Countomissions, CountCommissions; double ZmatchedmeanTrees, ZmatchedmeanClus; CountMatches = 0; Countomissions = 0; CountCommissions=0; ZmatchedmeanTrees = 0.0; ZmatchedmeanClus = 0.0; dXave = 0.0; dYave = 0.0; dZave = 0.0; //m = sprintf(buffer,"In MyFunc_FindClusters N_tree_inside_border = %d\n",N_tree_inside_border); //MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); for (j = 1; j <= N_tree_inside_border; j++) // Loop thru trees { if (TreeMatchStruct[j].IsInside == TRUE) { ZmatchedmeanTrees = ZmatchedmeanTrees + TreeMatchStruct[j].Z; switch (TreeMatchStruct[j].IsMatch) { case TRUE: CountMatches = CountMatches + 1; dXave = dXave + (TreeMatchStruct[j].X - TreeMatchStruct[j].Xclus); dYave = dYave + (TreeMatchStruct[j].Y - TreeMatchStruct[j].Yclus); dZave = dZave + (TreeMatchStruct[j].Z - TreeMatchStruct[j].Zclus); break; case FALSE: Countomissions = Countomissions + 1; break; default: m = sprintf(buffer,"In MyFunc_FindClusters this should not be\n"); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); break; } // Switch } // IsInside } // m = sprintf(buffer,"In MyFunc_FindClusters Countomissions %d\n",Countomissions); // MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); if (CountMatches != 0) { dXave = dXave / double(CountMatches); dYave = dYave / double(CountMatches); dZave = dZave / double(CountMatches); } if (CountMatches == 0) { dXave = -9.0; dYave = -9.0; dZave = -9.0; } for (j = 1; j <= N_cand_inside_border; j++) // ' Loop thru clusters { if (ClusMatchStruct[j].IsInside == TRUE) { ZmatchedmeanClus = ZmatchedmeanClus + ClusMatchStruct[j].Z;} if ((ClusMatchStruct[j].IsCommission == TRUE) && (ClusMatchStruct[j].IsInside == TRUE)) {CountCommissions = CountCommissions + 1;} } double Zbias; if (N_tree_inside != 0) {ZmatchedmeanTrees = ZmatchedmeanTrees / double(N_tree_inside);} if (N_cand_inside != 0) {ZmatchedmeanClus = ZmatchedmeanClus / double(N_cand_inside);} Zbias = ZmatchedmeanTrees - ZmatchedmeanClus; // ***************** RESULTS PRINTING SECTION **************************** // // Determine a linear LS-regression line for the Z-errors with respect to tree height double Zer, XYer, RMSEXY, RMSEZ, sx, sy, Ymean, VolMatch, VolMissed; RMSEXY =0.0; RMSEZ = 0.0; sx = 0.0; sy = 0.0; VolMatch = 0.0; VolMissed = 0.0, Ymean = 0.0; for (i = 1; i <= N_tree_inside_border; i++) // Loop thru trees inside the plot area { if ((TreeMatchStruct[i].IsMatch == TRUE) && (TreeMatchStruct[i].IsInside == TRUE)) { Zer = TreeMatchStruct[i].Z - TreeMatchStruct[i].Zclus; XYer = sqrt(pow((TreeMatchStruct[i].Y - TreeMatchStruct[i].Yclus) ,2) + pow((TreeMatchStruct[i].X - TreeMatchStruct[i].Xclus) , 2)); RMSEXY = RMSEXY + XYer; RMSEZ = RMSEZ + Zer * Zer; sx = sx + TreeMatchStruct[i].h; sy = sy + Zer; Ymean = Ymean + Zer; VolMatch = VolMatch + TreeMatchStruct[i].Vol; } if ((TreeMatchStruct[i].IsMatch == FALSE) && (TreeMatchStruct[i].IsInside == TRUE)) { VolMissed = VolMissed + TreeMatchStruct[i].Vol; } } if (CountMatches != 0) { RMSEXY = sqrt(RMSEXY / double(CountMatches)); RMSEZ = sqrt(RMSEZ / double(CountMatches)); } double sx0ss=0.0, t=0.0, st2, bb, aa=0.0, siga=0.0, sigb=0.0; st2 = 0.0; bb = 0.0; if (CountMatches != 0) { sx0ss = sx / double (CountMatches); Ymean = Ymean / double (CountMatches); } for (i = 1; i <= N_tree_inside_border; i++) // Loop thru trees inside the plot area { if ((TreeMatchStruct[i].IsMatch == TRUE) && (TreeMatchStruct[i].IsInside == TRUE)) { t = (TreeMatchStruct[i].h - sx0ss); st2 = st2 + t * t; bb = bb + t * (TreeMatchStruct[i].Z - TreeMatchStruct[i].Zclus); } } if ((CountMatches > 2) && st2 > (10E-19)) { bb = bb / st2; aa = (sy - sx * bb) / double (CountMatches); siga = sqrt ( (1.0 + sx * sx / (double (CountMatches) * st2)) / double (CountMatches) ); sigb = sqrt ( (1.0 / st2) / (double (CountMatches) - 2) ); } double ssres, ssreg; ssres = 0.0; ssreg = 0.0; for (i = 1; i <= N_tree_inside_border; i++) // Loop thru trees inside the plot area { if ((TreeMatchStruct[i].IsMatch == TRUE) && (TreeMatchStruct[i].IsInside == TRUE)) { ssres = ssres + pow ( (TreeMatchStruct[i].Z - TreeMatchStruct[i].Zclus) - (aa + bb * TreeMatchStruct[i].h) , 2); ssreg = ssreg + pow(((aa + bb * TreeMatchStruct[i].h) - Ymean) , 2); } } double r2=0.0, sres=0.0, seb=0.0; if ((CountMatches > 2) && st2 > (10E-19)) { r2 = ssreg / (ssreg + ssres); sres = sqrt ( ssres / (double (CountMatches) - 2.0) ); seb = sres / sqrt(st2); } fp = fopen(filename_out, "rb+"); if (fp==NULL) { ij = sprintf(buffer,"Output-tiedoston avaus ei onnistu\n"); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); free(TreeMatchStruct); free(ClusMatchStruct); return double(mc); } double MatchP = 0.0; MatchP = (double(CountMatches) / double(N_tree_inside))*100.0; time_t seconds; struct tm *p; time(&seconds); p = localtime(&seconds); p->tm_mon = p->tm_mon + 1; p->tm_year = (p->tm_year -100) + 2000; int a; a = fseek(fp, 0, SEEK_END); // Put the position at end + 1, this is where we write next unsigned char charapu; short shortapu; float floatapu; charapu = unsigned char(p->tm_mday); a = fwrite(&charapu, sizeof(unsigned char), 1, fp); charapu = unsigned char(p->tm_mon); a = fwrite(&charapu, sizeof(unsigned char), 1, fp); shortapu = short(p->tm_year); a = fwrite(&shortapu, sizeof(short), 1, fp); charapu = unsigned char(p->tm_hour); a = fwrite(&charapu, sizeof(unsigned char), 1, fp); charapu = unsigned char(p->tm_min); a = fwrite(&charapu, sizeof(unsigned char), 1, fp); charapu = unsigned char(p->tm_sec); a = fwrite(&charapu, sizeof(unsigned char), 1, fp); // 5 x 1 bytes, 1 x 2 bytes = 7 bytes charapu = unsigned char(NimagesInMatch); a = fwrite(&charapu, sizeof(unsigned char), 1, fp); shortapu = short(ModelTreeNum); a = fwrite(&shortapu, sizeof(short), 1, fp); a = fwrite(&ModelTreeX, sizeof(double), 1, fp); a = fwrite(&ModelTreeY, sizeof(double), 1, fp); a = fwrite(&ModelTreeZ, sizeof(double), 1, fp); // 1 x 1 bytes, 1 x 2 + 3 x 8 bytes = 27 bytes (Tot 34) floatapu= float(r_limit); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(xy_thin); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(Match2Ddist); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(MatchZdist); a = fwrite(&floatapu, sizeof(float), 1, fp); // 4 x 4 bytes = 16 (Total 50) floatapu= float(TWidth); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(EllZasym); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(Zdiff); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(ZDepth); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(Zasym); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(Meshdist); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu= float(Meanheight); a = fwrite(&floatapu, sizeof(float), 1, fp); // 7 x 4 bytes = 28 (Tot 78) floatapu= float(plotradius); a = fwrite(&floatapu, sizeof(float), 1, fp); a = fwrite(&Plotcenter, sizeof(Plotcenter), 1, fp); // 1 x 4 , 3 x 8 bytes = 28 (Tot 106) shortapu = short(N_tree_inside); a = fwrite(&shortapu, sizeof(short), 1, fp); shortapu = short(N_cand_inside); a = fwrite(&shortapu, sizeof(short), 1, fp); shortapu = short(CountMatches); a = fwrite(&shortapu, sizeof(short), 1, fp); shortapu = short(Countomissions); a = fwrite(&shortapu, sizeof(short), 1, fp); shortapu = short(CountCommissions); a = fwrite(&shortapu, sizeof(short), 1, fp); floatapu = float(MatchP); a = fwrite(&floatapu, sizeof(float), 1, fp); // 5 x 2 bytes, 1 x 4 bytes = 14 (Tot 120) floatapu = float(dZave); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(Zbias); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(dXave); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(dYave); a = fwrite(&floatapu, sizeof(float), 1, fp); // 4 x 4 bytes = 16 (Tot 136) floatapu = float(RMSEXY); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(RMSEZ); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(aa); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(bb); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(seb); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(r2); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(VolMatch); a = fwrite(&floatapu, sizeof(float), 1, fp); floatapu = float(VolMissed); a = fwrite(&floatapu, sizeof(float), 1, fp); // 8 x 4 bytes = 32 (Tot 158) // This section is used if ASCII output is desired, binary is used for speeding up i/o /* i = sprintf(&CO[0], "\"P3\",\"%02i.%02i.20%02i\",\"%02i:%02i:%02i\",\0", p->tm_mday, p->tm_mon+1, p->tm_year-100,p->tm_hour ,p->tm_min , p->tm_sec ); i = sprintf(&CO1[0], "%02i,%03i,%10.2f,%10.2f,%5.2f,\0", NimagesInMatch, ModelTreeNum,ModelTreeX,ModelTreeY,ModelTreeZ); strcat(CO,CO1); i = sprintf(&CO1[0], "%5.2f,%5.2f,%5.2f,%5.2f,\0",r_limit,xy_thin,Match2Ddist,MatchZdist); strcat(CO,CO1); i = sprintf(&CO1[0], "%4.2f,%4.2f,%4.2f,%4.2f,%4.2f,%4.2f,%4.2f,\0",TWidth,EllZasym,Zdiff,ZDepth,Zasym,Meshdist,Meanheight); strcat(CO,CO1); i = sprintf(&CO1[0], "%4.2f,%10.2f,%10.2f,%5.2f,\0",plotradius,Plotcenter.X,Plotcenter.Y,Plotcenter.Z); strcat(CO,CO1); i = sprintf (&CO1[0], "00000000,\0"); strcat(CO,CO1); i = sprintf(&CO1[0], "%3i,%3i,%3i,%3i,%3i,%5.2f,\0",N_tree_inside,N_cand_inside,CountMatches,Countomissions,CountCommissions,MatchP); strcat(CO,CO1); i = sprintf(&CO1[0], "%5.2f,%5.2f,%5.2f,%5.2f,\0",dZave,Zbias,dXave,dYave); strcat(CO,CO1); i = sprintf(&CO1[0], "%5.2f,%5.2f,%8.5f,%8.5f,%8.6f,%6.4f,%8.2f,%8.2f\n\0",RMSEXY, RMSEZ, aa,bb, seb, r2, VolMatch, VolMissed); strcat(CO,CO1); i = fputs(&CO[0],fp); */ fclose(fp); free(TreeMatchStruct); free(ClusMatchStruct); return double(mc); } static int kutsuttavafunktio(double *para) { int i; char buffer[100]; double a; a = *para; i = sprintf(buffer,"Kutsuit parametrilla ja globaali on %f\n", a); MessageBox ( NULL, buffer,"MyFuncMatlabCall",MB_OK); return 1; } #define NRANSI #include "nrutil.h" #define SWAP(a,b) itemp=(a);(a)=(b);(b)=itemp; #define M 7 #define NSTACK 50 double __stdcall MyFunc_indexx(long n, float *arr, long *indx) { long i,indxt,ir=n,itemp,j,k,l=1; int jstack=0,*istack; float a; istack=ivector(1,NSTACK); for (j=1;j<=n;j++) indx[j]=j; for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { indxt=indx[j]; a=arr[indxt]; for (i=j-1;i>=1;i--) { if (arr[indx[i]] <= a) break; indx[i+1]=indx[i]; } indx[i+1]=indxt; } if (jstack == 0) break; ir=istack[jstack--]; l=istack[jstack--]; } else { k=(l+ir) >> 1; SWAP(indx[k],indx[l+1]); if (arr[indx[l+1]] > arr[indx[ir]]) { SWAP(indx[l+1],indx[ir]) } if (arr[indx[l]] > arr[indx[ir]]) { SWAP(indx[l],indx[ir]) } if (arr[indx[l+1]] > arr[indx[l]]) { SWAP(indx[l+1],indx[l]) } i=l+1; j=ir; indxt=indx[l]; a=arr[indxt]; for (;;) { do i++; while (arr[indx[i]] < a); do j--; while (arr[indx[j]] > a); if (j < i) break; SWAP(indx[i],indx[j]) } indx[l]=indx[j]; indx[j]=indxt; jstack += 2; if (jstack > NSTACK) nrerror("NSTACK too small in indexx."); if (ir-i+1 >= j-l) { istack[jstack]=ir; istack[jstack-1]=i; ir=j-1; } else { istack[jstack]=j-1; istack[jstack-1]=l; l=i; } } } free_ivector(istack,1,NSTACK); return 1.0; } #undef M #undef NSTACK #undef SWAP #undef NRANSI /* Below three functions from nrutil.c needed by indexx above */ #include #include #include #define NR_END 1 #define FREE_ARG char* void nrerror(char error_text[]) /* Numerical Recipes standard error handler */ { fprintf(stderr,"Numerical Recipes run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); exit(1); } int *ivector(long nl, long nh) /* allocate an int vector with subscript range v[nl..nh] */ { int *v; v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); if (!v) nrerror("allocation failure in ivector()"); return v-nl+NR_END; } void free_ivector(int *v, long nl, long nh) /* free an int vector allocated with ivector() */ { free((FREE_ARG) (v+nl-NR_END)); } /* (C) Copr. 1986-92 Numerical Recipes Software "!15L1. */