00001 #include "sadie.h"
00002
00003
00004 typedef struct {
00005 unsigned short hstart;
00006 unsigned short hend;
00007 unsigned short *vstart;
00008 unsigned short *vend;
00009 } LOCAL_INDEX;
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 void ROI_GRADIENT (
00033 IMAGE *in,
00034 LOCAL_INDEX *localindex,
00035 PIXEL *vert,
00036
00037
00038 PIXEL *horz,
00039
00040
00041 short size,
00042 IMAGE **out1,
00043
00044 IMAGE **out2
00045
00046
00047 ) { register int i, j, k, l, m, n=size/2;
00048 char msg[SLEN];
00049 double dv, dh, pinc, psum;
00050
00051 int progmeter;
00052 double progress;
00053 int progcount, progtotal, progincr;
00054 int cancel=0;
00055
00056 progmeter = ProgMeter_Create("Gradient");
00057 progress = 0.0;
00058
00059
00060 if (TIMES) TIMING(T_GRADIENT);
00061 if (NAMES) {
00062 MESSAGE('I',"");
00063 MESSAGE('I',"GRADIENT");
00064 MESSAGE('I',"");
00065 sprintf(msg," Input image: %s",in->text);
00066 MESSAGE('I',msg);
00067 sprintf(msg," Kernel size: lines and pixels: %d",size);
00068 MESSAGE('I',msg);
00069 if (size <= SLEN/16) {
00070 MESSAGE('I'," Horizontal kernel:");
00071 for (i=0; i<size; i++) {
00072 for (j=k=0,msg[k++]=' '; j<size; j+=1,k=strlen(msg)) {
00073 sprintf(msg+k,"%12.4e",(double)horz[i*size+j]);
00074 }
00075 MESSAGE('I',msg);
00076 }
00077 MESSAGE('I'," Vertical kernel:");
00078 for (i=0; i<size; i++) {
00079 for (j=k=0,msg[k++]=' '; j<size; j+=1,k=strlen(msg)) {
00080 sprintf(msg+k,"%12.4e",(double)vert[i*size+j]);
00081 }
00082 MESSAGE('I',msg);
00083 }
00084 } else {
00085 MESSAGE('I'," Kernels too large to list.");
00086 }
00087 MESSAGE('I'," ...............");
00088 }
00089
00090
00091 if (!CHECKIMG(in)) {
00092 MESSAGE('E'," Can't identify image.");
00093 goto the_end;
00094 } else if (size < 2 || size%2 == 0) {
00095 MESSAGE('E'," Size of convolution mask must be an odd number greater than two.");
00096 goto the_end;
00097 } else if (size > in->nlin || size > in->npix) {
00098 MESSAGE('E'," Image size must be equal to or greater than convolution mask size.");
00099 goto the_end;
00100 }
00101
00102
00103 if (!CHECKIMG(*out1)) GETMEM(in->nbnd,in->nlin,in->npix,out1);
00104 if (!*out1) goto the_end;
00105 if (!CHECKIMG(*out2)) GETMEM(in->nbnd,in->nlin,in->npix,out2);
00106 if (!*out2) goto the_end;
00107
00108
00109 for(j = 0; j < in->nlin; j++) {
00110 for(k = 0; k < in->npix; k++) {
00111 (*out1)->data[0][j][k] = (PIXEL)0;
00112 (*out2)->data[0][j][k] = (PIXEL)0;
00113 }
00114 }
00115
00116
00117 if (LINES && !PROGRESS(psum=0.0)) goto the_end;
00118 pinc = 1.0/(double)in->nbnd/(double)(in->nlin-2*n);
00119 progtotal = in->nbnd * (in->nlin-2*n);
00120 progcount = 0;
00121 progincr = rint((double)progtotal/20.0);
00122
00123
00124 for (k = localindex->hstart + n + 1; k <= localindex->hend - n - 1; k++) {
00125 for (j = localindex->vstart[k] + n + 1; j <= localindex->vend[k] - n - 1; j++) {
00126 dv = dh = 0.0;
00127 for (l=0; l<size; l++) {
00128 for (m=0; m<size; m++) {
00129 dv += (double)vert[l*size+m]*(double)in->data[0][j-n+l][k-n+m];
00130 dh += (double)horz[l*size+m]*(double)in->data[0][j-n+l][k-n+m];
00131 }
00132 }
00133 (*out1)->data[0][j][k] = (PIXEL)sqrt((dv*dv+dh*dh));
00134 if (dh == 0.0) {
00135 if (dv > 0.0) {
00136 (*out2)->data[0][j][k] = (PIXEL)(PI/2.0);
00137 } else if (dv < 0.0) {
00138 (*out2)->data[0][j][k] = (PIXEL)(-PI/2.0);
00139 } else {
00140 (*out2)->data[0][j][k] = (PIXEL)0.0;
00141 }
00142 } else {
00143 (*out2)->data[0][j][k] = (PIXEL)atan2(dv,dh);
00144 }
00145 }
00146 if (LINES && !PROGRESS(psum+=pinc)) goto the_end;
00147 if (progmeter) {
00148 progcount++;
00149 if ((progcount%progincr) == 0) {
00150 progress = ((double)progcount/(double)progtotal)*100.0;
00151 cancel = ProgMeter_Update(progmeter,progress);
00152 if (cancel != 0) goto the_end;
00153 }
00154 }
00155 }
00156
00157
00158 the_end:
00159 if (TIMES) TIMING(T_EXIT);
00160
00161 if (progmeter) {
00162 ProgMeter_Destroy(progmeter);
00163 }
00164
00165 if (cancel != 0) {
00166 if (*out1) RELMEM(*out1);
00167 *out1 = NULL;
00168 if (*out2) RELMEM(*out2);
00169 *out2 = NULL;
00170 }
00171 }