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