00001 #include "sadie.h"
00002 #include "proto.h"
00003 #include <math.h>
00004
00005
00006
00007
00008
00009
00010 extern int **x_array_buffer;
00011 extern int point;
00012 extern int count_point;
00013
00014
00015 typedef struct {
00016 short x;
00017 short y;
00018 double graddir;
00019 double avgdir;
00020 } RINGPIXEL_DATA;
00021
00022 typedef struct {
00023 short counter;
00024 PIXEL *array;
00025 double width;
00026 double hwidth;
00027 } WIDTHLINE_DATA;
00028
00029 typedef struct {
00030 short size;
00031 short counter;
00032 RINGPIXEL_DATA *pixels;
00033 int numwsamp;
00034 WIDTHLINE_DATA *wsample;
00035 double avgwidth;
00036 double avghwidth;
00037 double index;
00038 double hindex;
00039 PIXEL *avgprofile;
00040 } TREERING_DATA;
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 void RINGWIDTHS_NOCHECK (
00054 IMAGE *in,
00055 IMAGE *rings,
00056 IMAGE *graddir,
00057 IMAGE **out
00058
00059 ) { register short i, j, k, l, n;
00060 char msg[SLEN];
00061 IMAGE *lblimg;
00062 TREERING_DATA *ringdata=NULL;
00063 int curlbl;
00064 double theta, curtheta, deltatheta;
00065 double tempdata, tempdata2;
00066 double targetslope;
00067 double curslope;
00068 short startx, starty, endx, endy;
00069 short deltay, deltax;
00070 short wcount;
00071 short hwcount;
00072 int width;
00073 double percent;
00074 int sampindex;
00075 PIXEL *avgprofile=NULL;
00076 int avgprofsize;
00077 int nbrhdsize = 35;
00078 PIXEL nbrhd[nbrhdsize];
00079 int index;
00080 PIXEL med;
00081
00082 short nlev;
00083 short inc = 1;
00084 double mean;
00085 double dev;
00086 double max;
00087 double min;
00088 double range;
00089 double sum, sumsq;
00090 int zeros;
00091 PIXEL *hist=NULL;
00092
00093 double hmean;
00094 double hdev;
00095 double hmax;
00096 double hmin;
00097 double hrange;
00098 double hsum, hsumsq;
00099 int hzeros;
00100
00101 double avgringwidth;
00102 double avghringwidth;
00103
00104
00105
00106
00107 if (!CHECKIMG(rings)) {
00108 MESSAGE('E'," Can't identify rings image.");
00109 goto the_end;
00110 } else if (!CHECKIMG(graddir)) {
00111 MESSAGE('E'," Can't identify gradient direction image.");
00112 goto the_end;
00113 }
00114
00115
00116 if (NAMES) {
00117 MESSAGE('I',"");
00118 MESSAGE('I',"RINGWIDTHS_NOCHECK");
00119 MESSAGE('I',"");
00120 sprintf(msg," Rings image: %s",rings->text);
00121 MESSAGE('I',msg);
00122 MESSAGE('I',"");
00123 sprintf(msg," Gradient Direction image: %s",graddir->text);
00124 MESSAGE('I',msg);
00125 MESSAGE('I'," ...............");
00126 }
00127
00128
00129 if (rings->nbnd > 1) {
00130 MESSAGE('E'," Only using first band of rings image.");
00131 }
00132 if (graddir->nbnd > 1) {
00133 MESSAGE('E'," Only using first band of gradient direction image.");
00134 }
00135 if ((in->nlin != graddir->nlin) || (in->npix != graddir->npix)) {
00136 MESSAGE('E'," Input images must be same size.");
00137 goto the_end;
00138 }
00139 if ((rings->nlin != graddir->nlin) || (rings->npix != graddir->npix)) {
00140 MESSAGE('E'," Input images must be same size.");
00141 goto the_end;
00142 }
00143
00144 RANGE(graddir);
00145 if (graddir->gmin < -2.0*PI || 2.0*PI < graddir->gmax) {
00146 MESSAGE('E'," This is not a gradient direction image!");
00147 MESSAGE('E'," Direction must be between -360 and +360 degrees.");
00148 goto the_end;
00149 }
00150
00151
00152
00153 if (!CHECKIMG(*out)) GETMEM(1,rings->nlin,rings->npix,out);
00154 if (!*out) goto the_end;
00155
00156
00157 for (j=0; j<rings->nlin; j++) {
00158 for (k=0; k<rings->npix; k++) {
00159 (*out)->data[0][j][k] = (PIXEL)0.0;
00160 }
00161 }
00162
00163
00164
00165
00166 CMPLBL8(rings,&lblimg);
00167 if (!CHECKIMG(lblimg)) {
00168 MESSAGE('E'," Problem labeling rings image.");
00169 goto the_end;
00170 }
00171 RANGE(lblimg);
00172
00173
00174
00175
00176 if (!(ringdata=(TREERING_DATA *)malloc((int)rnd(lblimg->gmax)*sizeof(TREERING_DATA)))) {
00177 MESSAGE('E'," Could not allocate memory for treering data.");
00178 goto the_end;
00179 }
00180
00181
00182 for (i = (int)rnd(lblimg->gmin); i < (int)rnd(lblimg->gmax); i++) {
00183 ringdata[i].size = 0;
00184 ringdata[i].counter = 0;
00185 ringdata[i].pixels = NULL;
00186 ringdata[i].numwsamp = 0;
00187 ringdata[i].wsample = NULL;
00188 ringdata[i].avgwidth = 0.0;
00189 ringdata[i].avghwidth = 0.0;
00190 ringdata[i].index = 0.0;
00191 ringdata[i].hindex = 0.0;
00192 ringdata[i].avgprofile = NULL;
00193 }
00194
00195
00196 for (j=0; j<lblimg->nlin; j++) {
00197 for (k=0; k<lblimg->npix; k++) {
00198 if (lblimg->data[0][j][k] > (PIXEL)0.0) {
00199 curlbl = (int)rnd(lblimg->data[0][j][k]) - 1;
00200 ringdata[curlbl].size += 1;
00201 }
00202 }
00203 }
00204
00205
00206 for (i = (int)rnd(lblimg->gmin); i < (int)rnd(lblimg->gmax); i++) {
00207 ringdata[i].pixels = (RINGPIXEL_DATA *) malloc((long)(ringdata[i].size)*sizeof(RINGPIXEL_DATA));
00208 if (ringdata[i].pixels == (RINGPIXEL_DATA *) NULL) {
00209 MESSAGE('E'," Could not allocate memory for treering pixel data.");
00210 goto the_end;
00211 }
00212
00213 ringdata[i].wsample = (WIDTHLINE_DATA *) malloc((long)(ringdata[i].size)*sizeof(WIDTHLINE_DATA));
00214 if (ringdata[i].wsample == (WIDTHLINE_DATA *) NULL) {
00215 MESSAGE('E'," Could not allocate memory for width sample profile data.");
00216 goto the_end;
00217 }
00218 for (j=0; j<ringdata[i].size; j++) {
00219 ringdata[i].wsample[j].counter = 0;
00220 ringdata[i].wsample[j].array = NULL;
00221 ringdata[i].wsample[j].width = 0.0;
00222 ringdata[i].wsample[j].hwidth = 0.0;
00223 }
00224 }
00225
00226
00227
00228 for (j=0; j<lblimg->nlin; j++) {
00229 for (k=0; k<lblimg->npix; k++) {
00230 if (lblimg->data[0][j][k] > (PIXEL)0.0) {
00231 curlbl = (int)rnd(lblimg->data[0][j][k]) - 1;
00232 ringdata[curlbl].pixels[ringdata[curlbl].counter].x = k;
00233 ringdata[curlbl].pixels[ringdata[curlbl].counter].y = j;
00234
00235 theta = (double)(graddir->data[0][j][k]);
00236
00237
00238
00239
00240
00241
00242 ringdata[curlbl].pixels[ringdata[curlbl].counter].graddir = theta;
00243
00244 ringdata[curlbl].counter += 1;
00245 }
00246 }
00247 }
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 for (i = (int)rnd(lblimg->gmin); i < (int)rnd(lblimg->gmax); i++) {
00278 if (ringdata[i].size > 10) {
00279 tempdata = 0.0;
00280 for (l=0; l<=10; l++) {
00281 ringdata[i].pixels[l].avgdir = 0.0;
00282 }
00283 for (l=((nbrhdsize+1)/2); l<ringdata[i].size-(((nbrhdsize+1)/2)+1); l++) {
00284 for (j=0; j<nbrhdsize; j++) {
00285 nbrhd[j] = ringdata[i].pixels[l+((nbrhdsize+1)/2)-j].graddir;
00286 }
00287
00288
00289
00290 for (j=0; j<(nbrhdsize+1)/2; j++) {
00291 med = nbrhd[index=j];
00292 for (k=j+1; k<nbrhdsize; k++) {
00293 if (med > nbrhd[k]) {
00294 med = nbrhd[index=k];
00295 }
00296 }
00297 nbrhd[index] = nbrhd[j];
00298 }
00299 ringdata[i].pixels[l].avgdir = med;
00300
00301
00302 }
00303 for (; l<ringdata[i].size; l++) {
00304 ringdata[i].pixels[l].avgdir = 0.0;
00305 }
00306 } else {
00307 for (l=0; l<ringdata[i].size; l++) {
00308 ringdata[i].pixels[l].avgdir = 0.0;
00309 }
00310 }
00311 }
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324 printf("\n");
00325 for (i = (int)rnd(lblimg->gmin); i < (int)rnd(lblimg->gmax)-1; i++) {
00326 wcount = 0;
00327 hwcount = 0;
00328 tempdata = 0.0;
00329 tempdata2 = 0.0;
00330
00331 for (l=0; l<ringdata[i].size; l++) {
00332 if (l%1 == 0) {
00333 theta = (double) ringdata[i].pixels[l].avgdir;
00334 targetslope = tan(theta);
00335
00336
00337
00338 if (l < ringdata[i+1].size) {
00339 j = l;
00340 } else {
00341 j = ringdata[i+1].size - 1;
00342 }
00343
00344
00345
00346 if (ringdata[i+1].pixels[j].y < ringdata[i].pixels[l].y) {
00347 while (ringdata[i+1].pixels[j].y < ringdata[i].pixels[l].y) {
00348 j++;
00349 }
00350 } else {
00351 while (ringdata[i+1].pixels[j].y > ringdata[i].pixels[l].y) {
00352 j--;
00353 }
00354 }
00355 ringdata[i].wsample[hwcount].hwidth = (double)(ringdata[i+1].pixels[j].x - ringdata[i].pixels[l].x);
00356 tempdata2 += ringdata[i].wsample[hwcount].hwidth;
00357 hwcount++;
00358
00359
00360
00361 deltax = ringdata[i+1].pixels[j].x - ringdata[i].pixels[l].x;
00362 deltay = ringdata[i].pixels[l].y - ringdata[i+1].pixels[j].y;
00363 curslope = (double)(deltay) / (double)(deltax);
00364 if (curslope > targetslope) {
00365 while ((curslope > targetslope) && (j < ringdata[i+1].size-1)) {
00366 deltax = ringdata[i+1].pixels[j].x - ringdata[i].pixels[l].x;
00367 deltay = ringdata[i].pixels[l].y - ringdata[i+1].pixels[j].y;
00368 curslope = (double)(deltay) / (double)(deltax);
00369 (double)(deltax);
00370
00371
00372
00373
00374
00375 j++;
00376 }
00377 } else {
00378 while ((curslope < targetslope) && (j > 0)) {
00379 deltax = ringdata[i+1].pixels[j].x - ringdata[i].pixels[l].x;
00380 deltay = ringdata[i].pixels[l].y - ringdata[i+1].pixels[j].y;
00381 curslope = (double)(deltay) / (double)(deltax);
00382
00383
00384
00385
00386
00387 j--;
00388 }
00389 }
00390
00391
00392
00393 curtheta = (double)ringdata[i+1].pixels[j].avgdir;
00394 deltatheta = curtheta - theta;
00395 if (deltatheta < 0.0) {
00396 deltatheta = 0.0 - deltatheta;
00397 }
00398
00399
00400
00401
00402
00403
00404 ringdata[i].wsample[ringdata[i].numwsamp].width = sqrt((double)(deltax*deltax)+(double)(deltay*deltay));
00405 tempdata += ringdata[i].wsample[ringdata[i].numwsamp].width;
00406 wcount++;
00407
00408
00409
00410
00411
00412 startx = ringdata[i].pixels[l].x;
00413 starty = ringdata[i].pixels[l].y;
00414 endx = ringdata[i+1].pixels[j].x;
00415 endy = ringdata[i+1].pixels[j].y;
00416
00417
00418
00419
00420 count_point = TRUE;
00421 point = 0;
00422 brshnm(startx,starty,endx, endy);
00423
00424
00425 x_array_buffer = (int **) malloc(point*sizeof(int * ));
00426 if (x_array_buffer == (int **) NULL){
00427 printf("No memory for x_array_buffer\n");
00428 } else{
00429
00430 for (n = 0; n < point; n++){
00431 x_array_buffer[n] = (int *) calloc(2, sizeof(int));
00432 if (x_array_buffer[n] == (int) NULL){
00433 printf ("No memory for element %d\n",x_array_buffer[n]);
00434 }
00435 }
00436
00437
00438 count_point = FALSE;
00439 point = 0;
00440 brshnm(startx,starty,endx, endy);
00441
00442
00443
00444 ringdata[i].wsample[ringdata[i].numwsamp].counter = point;
00445 ringdata[i].wsample[ringdata[i].numwsamp].array = (PIXEL *) malloc(point * sizeof(PIXEL));
00446 for (n=0; n<point; n++) {
00447 ringdata[i].wsample[ringdata[i].numwsamp].array[n] = in->data[0][x_array_buffer[n][1]][x_array_buffer[n][0]];
00448 }
00449 ringdata[i].numwsamp += 1;
00450
00451
00452 if (l%10 == 0) {
00453
00454 for (n=0; n < point; n++) {
00455 (*out)->data[0][x_array_buffer[n][1]][x_array_buffer[n][0]] = (PIXEL)1.0;
00456 }
00457 }
00458
00459 if(x_array_buffer) {
00460
00461 for (n=0; n < point; n++){
00462 if(x_array_buffer[n]) {
00463 free(x_array_buffer[n]);
00464 x_array_buffer[n] = NULL;
00465 }
00466 }
00467
00468 free(x_array_buffer);
00469 x_array_buffer = NULL;
00470 }
00471
00472 }
00473
00474
00475
00476
00477 }
00478 }
00479
00480
00481 if (wcount > 0) {
00482 ringdata[i].avgwidth = tempdata / (double)wcount;
00483 } else {
00484 ringdata[i].avgwidth = 0.0;
00485 }
00486
00487 if (hwcount > 0) {
00488 ringdata[i].avghwidth = tempdata2 / (double)hwcount;
00489 } else {
00490 ringdata[i].avghwidth = 0.0;
00491 }
00492
00493
00494
00495
00496 if (ringdata[i].avgwidth > 1.0) {
00497 width = (int)rnd(ringdata[i].avgwidth);
00498 ringdata[i].avgprofile = (PIXEL *) malloc(width * sizeof(PIXEL));
00499
00500 for (n=0; n<width; n++) {
00501 ringdata[i].avgprofile[n] = (PIXEL)0.0;
00502 }
00503
00504
00505 min = max = ringdata[i].wsample[0].width;
00506 zeros = 0;
00507 hmin = hmax = ringdata[i].wsample[0].hwidth;
00508 hzeros = 0;
00509 for (sum=sumsq=0.0,hsum=hsumsq=0.0,j=0; j<ringdata[i].numwsamp; j++) {
00510
00511 sum += ringdata[i].wsample[j].width;
00512 sumsq += (ringdata[i].wsample[j].width)*(ringdata[i].wsample[j].width);
00513 min = min(min,ringdata[i].wsample[j].width);
00514 max = max(max,ringdata[i].wsample[j].width);
00515 if (ringdata[i].wsample[j].width == 0.0) {
00516 zeros++;
00517 }
00518
00519 hsum += ringdata[i].wsample[j].hwidth;
00520 hsumsq += (ringdata[i].wsample[j].hwidth)*(ringdata[i].wsample[j].hwidth);
00521 hmin = min(hmin,ringdata[i].wsample[j].hwidth);
00522 hmax = max(hmax,ringdata[i].wsample[j].hwidth);
00523 if (ringdata[i].wsample[j].hwidth == 0.0) {
00524 hzeros++;
00525 }
00526
00527
00528
00529
00530
00531 for (n=0; n<width; n++) {
00532 percent = (double)(n+1) / (double)width;
00533 sampindex = (int)rnd((double)ringdata[i].wsample[j].counter * percent) - 1;
00534 if (sampindex < 0) sampindex = 0;
00535 if (sampindex >= ringdata[i].wsample[j].counter) sampindex = ringdata[i].wsample[j].counter-1;
00536 ringdata[i].avgprofile[n] += ringdata[i].wsample[j].array[sampindex];
00537 }
00538 }
00539 mean = sum / (double)ringdata[i].numwsamp;
00540 dev = sqrt((sumsq / (double)ringdata[i].numwsamp) - (mean*mean));
00541 range = max - min;
00542
00543 hmean = hsum / (double)ringdata[i].numwsamp;
00544 hdev = sqrt((hsumsq / (double)ringdata[i].numwsamp) - (hmean*hmean));
00545 hrange = hmax - hmin;
00546
00547 MESSAGE('I',"");
00548 sprintf(msg," Statistics for orthogonal width measurements of Ring %d:",i+1);
00549 MESSAGE('I',msg);
00550 MESSAGE('I'," Minimum Maximum Mean Deviation Zero/Non-Zero");
00551 sprintf(msg," %12.4e%12.4e%12.4e%12.4e%10ld/%ld",min,max,mean,dev,(long)zeros,(long)ringdata[i].numwsamp-(long)zeros);
00552 MESSAGE('I',msg);
00553 sprintf(msg," Statistics for horizontal width measurements of Ring %d:",i+1);
00554 MESSAGE('I',msg);
00555 MESSAGE('I'," Minimum Maximum Mean Deviation Zero/Non-Zero");
00556 sprintf(msg," %12.4e%12.4e%12.4e%12.4e%10ld/%ld",hmin,hmax,hmean,hdev,(long)hzeros,(long)hwcount-(long)hzeros);
00557 MESSAGE('I',msg);
00558 MESSAGE('I',"");
00559
00560 for (n=0; n<width; n++) {
00561 ringdata[i].avgprofile[n] /= (PIXEL)(ringdata[i].numwsamp);
00562 }
00563 } else {
00564 ringdata[i].avgprofile = (PIXEL *) malloc(sizeof(PIXEL));
00565 ringdata[i].avgwidth = 1.0;
00566 ringdata[i].avghwidth = 1.0;
00567 ringdata[i].avgprofile[0] = (PIXEL) 0.0;
00568 }
00569
00570
00571 }
00572
00573
00574 tempdata = 0.0;
00575 tempdata2 = 0.0;
00576 for (i = 0; i < (int)rnd(lblimg->gmax)-1; i++) {
00577 tempdata += (double)ringdata[i].avgwidth;
00578 tempdata2 += (double)ringdata[i].avghwidth;
00579 }
00580 avgringwidth = (double) tempdata / ((double)(lblimg->gmax) - 1.0);
00581 avghringwidth = (double) tempdata2 / ((double)(lblimg->gmax) - 1.0);
00582 avgprofsize = (int)rnd(tempdata);
00583 avgprofile = (PIXEL *) malloc(avgprofsize * sizeof(PIXEL));
00584
00585 for (i=0,k=0; i < (int)rnd(lblimg->gmax)-1; i++) {
00586
00587 ringdata[i].index = ringdata[i].avgwidth / avgringwidth;
00588 ringdata[i].hindex = ringdata[i].avghwidth / avghringwidth;
00589
00590 for (j=0; j<(int)rnd(ringdata[i].avgwidth); j++) {
00591 avgprofile[k++] = ringdata[i].avgprofile[j];
00592 }
00593 }
00594
00595
00596 MESSAGE('I',"");
00597 MESSAGE('I',"");
00598 MESSAGE('I'," Ring Orth. width Orth. index Horiz. width Horiz. index");
00599 MESSAGE('I',"------ ------------- ------------- -------------- --------------");
00600 for (i=0; i < (int)rnd(lblimg->gmax)-1; i++) {
00601 sprintf(msg,"%3d %10.3f pixels %4.3f %10.3f pixels %4.3f",i+1,ringdata[i].avgwidth,ringdata[i].index,ringdata[i].avghwidth,ringdata[i].hindex);
00602 MESSAGE('I',msg);
00603 }
00604 MESSAGE('I',"");
00605 MESSAGE('I',"");
00606
00607
00608 Sadie_Plot_PlotArrayCmd(avgprofile, avgprofsize);
00609
00610
00611
00612
00613 the_end:
00614
00615
00616
00617 if (ringdata) {
00618 for (i = (int)rnd(lblimg->gmin); i < (int)rnd(lblimg->gmax); i++) {
00619 if(ringdata[i].pixels) {
00620 free(ringdata[i].pixels);
00621 ringdata[i].pixels = NULL;
00622 }
00623 if(ringdata[i].wsample) {
00624 for (j=0; j<ringdata[i].size; j++) {
00625 if (ringdata[i].wsample[j].array) {
00626 free(ringdata[i].wsample[j].array);
00627 ringdata[i].wsample[j].array = NULL;
00628 }
00629 }
00630
00631 free(ringdata[i].wsample);
00632 ringdata[i].wsample = NULL;
00633 }
00634 if(ringdata[i].avgprofile) {
00635 free(ringdata[i].avgprofile);
00636 ringdata[i].avgprofile = NULL;
00637 }
00638 }
00639
00640 free(ringdata);
00641 ringdata = NULL;
00642 }
00643
00644 if (avgprofile) {
00645 free(avgprofile);
00646 avgprofile = NULL;
00647 }
00648
00649 if (hist) {
00650 free(hist);
00651 hist = NULL;
00652 }
00653
00654 if (CHECKIMG(lblimg)) RELIMG(&lblimg);
00655
00656 }