00001 #include "sadie.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 typedef struct {
00017 unsigned short hstart;
00018 unsigned short hend;
00019 unsigned short *vstart;
00020 unsigned short *vend;
00021 } LOCAL_INDEX;
00022
00023
00024 static PIXEL firstdiff_x[3][3] = { { -1.0, 0.0, 1.0 },
00025 { -2.0, 0.0, 2.0 },
00026 { -1.0, 0.0, 1.0 } };
00027
00028 static PIXEL firstdiff_y[3][3] = { { 1.0, 2.0, 1.0 },
00029 { 0.0, 0.0, 0.0 },
00030 {-1.0, -2.0, -1.0 } };
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 void ROI_TREERINGMAG (
00041 IMAGE *in,
00042 LOCAL_INDEX *localindex,
00043 double sigma,
00044 int size,
00045 IMAGE **out
00046
00047
00048 ) { register int i, j, k, n;
00049 char msg[SLEN];
00050 IMAGE *gaussian, *smoothed;
00051 IMAGE *FoG_mag, *FoG_dir;
00052 PIXEL *psf=0, factor;
00053 PIXEL theta;
00054 int sectorval;
00055
00056
00057 int progmeter;
00058 double progress;
00059 double localpercent;
00060 double localprogress;
00061 int progcount, progtotal, progincr;
00062 int cancel=0;
00063
00064
00065 progmeter = ProgMeter_Create("Compute Tree Ring Magnitudes");
00066 progress = 0.0;
00067
00068
00069
00070
00071 if (NAMES) {
00072 MESSAGE('I',"");
00073 MESSAGE('I',"TREERINGMAG");
00074 MESSAGE('I',"");
00075 sprintf(msg," Input image: %s",in->text);
00076 MESSAGE('I',msg);
00077 sprintf(msg," Gaussian smoothing filter size: %d",size);
00078 MESSAGE('I',msg);
00079 sprintf(msg," Standard Deviation: %f",sigma);
00080 MESSAGE('I',msg);
00081 MESSAGE('I'," ...............");
00082 }
00083
00084
00085 if (!CHECKIMG(in)) {
00086 MESSAGE('E'," Can't identify image.");
00087 goto the_end;
00088 } else if (in->nbnd > 1) {
00089 MESSAGE('E'," Image must be single-band.");
00090 goto the_end;
00091 }
00092
00093
00094
00095 if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,in->npix,out);
00096 if (!*out) goto the_end;
00097
00098
00099 if (size > 0) {
00100
00101 CREATEGAUSS(size, sigma, &gaussian);
00102 if (CHECKIMG(gaussian)) {
00103 sprintf(gaussian->text, "%dx%d Gaussian -- sigma=%f", size,size,sigma);
00104 } else {
00105 MESSAGE('E'," Can't identify gaussian filter.");
00106 goto the_end;
00107 }
00108 if (progmeter) {
00109 progress = 10.0;
00110 cancel = ProgMeter_Update(progmeter,progress);
00111 if (cancel != 0) goto the_end;
00112 }
00113
00114
00115 if (!(psf = (PIXEL *) malloc(size*size*sizeof(PIXEL)))) {
00116 MESSAGE('E'," Memory allocation failed.");
00117 goto the_end;
00118 }
00119 for (n=0, j=0, factor=0.0; j<size; j++) {
00120 for (k=0; k<size; k++) {
00121 psf[n] = gaussian->data[0][j][k];
00122 factor += psf[n++];
00123 }
00124 }
00125 for (n=0, j=0; j<size; j++) {
00126 for (k=0; k<size; k++) {
00127 psf[n++] /= factor;
00128 }
00129 }
00130 ROI_SCONVL(in, localindex, psf, size, size, &smoothed);
00131 } else {
00132 IMGCOPY(in,&smoothed);
00133 }
00134
00135 if (progmeter) {
00136 progress = 30.0;
00137 cancel = ProgMeter_Update(progmeter,progress);
00138 if (cancel != 0) goto the_end;
00139 }
00140
00141
00142 if (CHECKIMG(smoothed)) {
00143 sprintf(smoothed->text, "Smoothed image");
00144 } else {
00145 MESSAGE('E'," Can't identify smoothed image.");
00146 goto the_end;
00147 }
00148
00149
00150 ROI_GRADIENT(smoothed, localindex, (PIXEL *)firstdiff_y, (PIXEL *)firstdiff_x, 3, &FoG_mag, &FoG_dir);
00151
00152 if (progmeter) {
00153 progress = 70.0;
00154 cancel = ProgMeter_Update(progmeter,progress);
00155 if (cancel != 0) goto the_end;
00156 }
00157
00158
00159
00160 if (progmeter) {
00161 localpercent = 0.3;
00162 localprogress = 0.0;
00163 progcount = 0;
00164 progtotal = in->npix * in->nlin;
00165 progincr = rint((double)progtotal/8.0);
00166 }
00167 for (j=0; j<in->nlin; j++) {
00168 for (k=0; k<in->npix; k++) {
00169
00170 if (progmeter) {
00171 progcount++;
00172 if ((progcount%progincr) == 0) {
00173 localprogress = ((double)progcount/(double)progtotal)*100.0;
00174 cancel = ProgMeter_Update(progmeter,(localprogress*localpercent)+progress);
00175 if (cancel != 0) goto the_end;
00176 }
00177 }
00178
00179 if (((k>0)&&(FoG_mag->data[0][j][k-1] > FoG_mag->data[0][j][k]))
00180 || ((k<in->npix-1)&&(FoG_mag->data[0][j][k+1] > FoG_mag->data[0][j][k]))) {
00181
00182 (*out)->data[0][j][k] = 0.0;
00183 } else if ((k<in->npix-1)&&(k>0)
00184 &&(smoothed->data[0][j][k+1] < smoothed->data[0][j][k-1])) {
00185
00186 (*out)->data[0][j][k] = 0.0;
00187 } else {
00188 (*out)->data[0][j][k] = FoG_mag->data[0][j][k];
00189 }
00190 }
00191 }
00192
00193
00194 the_end:
00195 if (psf) free(psf);
00196 if (CHECKIMG(gaussian)) RELIMG(&gaussian);
00197 if (CHECKIMG(smoothed)) RELIMG(&smoothed);
00198 if (CHECKIMG(FoG_mag)) RELIMG(&FoG_mag);
00199 if (CHECKIMG(FoG_dir)) RELIMG(&FoG_dir);
00200
00201
00202 if (progmeter) {
00203 ProgMeter_Destroy(progmeter);
00204 }
00205
00206 if (cancel != 0) {
00207 if (*out) RELMEM(*out);
00208 *out = NULL;
00209 }
00210 }