Changeset 4185
- Timestamp:
- Jun 8, 2005, 8:18:19 PM (21 years ago)
- Location:
- trunk/psModules/src
- Files:
-
- 2 edited
-
pmImageCombine.c (modified) (3 diffs)
-
pmImageSubtract.c (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmImageCombine.c
r4030 r4185 3 3 * This file will ... 4 4 * 5 * @author Paul Price, IfA 5 6 * @author GLG, MHPCC 6 7 * 7 * @version $Revision: 1. 4$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-0 5-25 20:28:32$8 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-06-09 06:16:43 $ 9 10 * 10 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 11 *12 * XXX: None of this has been tested.13 *14 *15 *16 *17 *18 12 * 19 13 */ … … 269 263 /****************************************************************************** 270 264 XXX: Directly from Paul Price 271 XXX: Confirm that this gradient calculation is still what's required.265 XXX: Must add mask parameter, use it in gradient calculation. 272 266 *****************************************************************************/ 273 267 static psF32 CalcGradient(psImage *image, … … 316 310 gradLimit) 317 311 318 XXX: Use static variables.319 312 XXX: Optimization: we don't need to transform the entire mask image. 320 313 XXX: The inToOut and outToIn transforms are confusing. Verify that what 321 314 I think they mean syncs with PWP. 315 XXX: Must add mask parameter, use it in gradient calculation. 322 316 *****************************************************************************/ 323 317 psArray *pmRejectPixels(const psArray *images, ///< Array of input images -
trunk/psModules/src/pmImageSubtract.c
r4168 r4185 4 4 * 5 5 * @author GLG, MHPCC 6 * @author Paul Price, IfA 6 7 * 7 * @version $Revision: 1. 7$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-06-09 0 0:59:53$8 * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-06-09 06:18:19 $ 9 10 * 10 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 52 53 psFree(stamp->matrix); 53 54 psFree(stamp->vector); 54 psFree(stamp);55 55 } 56 56 … … 392 392 393 393 // Determine if this pixel is larger than the max, and unmasked. 394 if ((image->data.F32[y][x] > max) && 395 !(mask->data.U8[y][x] & mask->data.U8[y][x])) { 396 max = image->data.F32[y][x]; 397 bestx = x; 398 besty = y; 394 if (image->data.F32[y][x] > max) { 395 if ((mask == NULL) || !((mask->data.U8[y][x]) & maskVal)) { 396 max = image->data.F32[y][x]; 397 bestx = x; 398 besty = y; 399 } 399 400 } 400 401 } … … 405 406 // Otherwise, mark the stamp as PM_STAMP_NONE 406 407 // 407 if (image->data.F32[besty][bestx] > threshold) {408 if (image->data.F32[besty][bestx] >= threshold) { 408 409 stamp->x = bestx; 409 410 stamp->y = besty; … … 451 452 452 453 /******************************************************************************* 454 GeneralKernelConvolve(input, kernels, kernelID, col, row): This routine 455 convolves a single kernel basis function with a pixel in an image. 456 457 XXX: merge this code with the other convolution code in this file. 458 459 XXX: Must we ensure that the pixels being accessed are inside the image? Is 460 there any garantee that the kernelSize, footprint stuff will be set correctly 461 to ensure no seg faults. 462 ******************************************************************************/ 463 psF32 GeneralKernelConvolve(const psImage *input, 464 const psSubtractionKernels *kernels, 465 psS32 kernelID, 466 psS32 col, 467 psS32 row) 468 { 469 psS32 spatialOrder = kernels->p_spatialOrder; 470 psS32 kernelSize = kernels->p_size; 471 psS32 xOrder = (psS32) kernels->xOrder->data.F32[kernelID]; 472 psS32 yOrder = (psS32) kernels->yOrder->data.F32[kernelID]; 473 psF32 numColsHalf = 0.5 * (psF32) input->numCols; 474 psF32 numRowsHalf = 0.5 * (psF32) input->numRows; 475 psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position 476 psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position 477 psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY); 478 psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; 479 480 psImage *preCalc = (psImage *) kernels->preCalc->data[kernelID]; 481 psF32 conv = 0.0; 482 for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { 483 for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { 484 //printf("HERE: (%d, %d)\n", yy, xx); 485 //printf("HERE: (%d, %d)\n", yy+row, xx+col); 486 //printf("HERE: (%d, %d)\n", yy-kernelSize, xx-kernelSize); 487 //printf("KERNEL SIZE is %d\n", kernelSize); 488 conv += input->data.F32[yy+row][xx+col] * 489 preCalc->data.F32[yy+kernelSize][xx+kernelSize] * 490 polyVal; 491 } 492 } 493 psFree(polyValues); 494 495 return(conv); 496 } 497 498 499 500 501 /******************************************************************************* 453 502 XXX: How is the spatial order factor calculated? Is it simply a x^iy^j power 454 503 evaluated at the (x, y) center for the stamp? Why bother with the 2-D 455 504 polynomial in evaluating it? 505 506 XXX: Should we assert that the footprint is equal to the kernel size, or does 507 they have nothing in common. 456 508 ******************************************************************************/ 457 509 bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, … … 494 546 for (psS32 s = 0; s < stamps->n; s++) { 495 547 pmStamp *stamp = (pmStamp *) stamps->data[s]; 496 548 psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d\n", s); 497 549 if (stamp->status == PM_STAMP_RECALC) { 550 psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: status is PM_STAMP_RECALC.\n", s); 498 551 psImage *stampMatrix = stamp->matrix; 499 552 psVector *stampVector = stamp->vector; … … 501 554 if (stampMatrix == NULL) { 502 555 stampMatrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); 556 stamp->matrix = stampMatrix; 503 557 } else { 504 558 PS_ASSERT_IMAGE_TYPE(stampMatrix, PS_TYPE_F64, false); … … 509 563 if (stampVector == NULL) { 510 564 stampVector = psVectorAlloc(numSolveParams, PS_TYPE_F64); 565 stamp->vector = stampVector; 511 566 } else { 512 567 PS_ASSERT_VECTOR_TYPE(stampVector, PS_TYPE_F64, false); … … 514 569 } 515 570 PS_VECTOR_SET_F64(stampVector, 0.0); 571 psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: allocate matrix and vector.\n", s); 516 572 517 573 // … … 524 580 ((psF64) (stamp->y - numHalfRows)) / ((psF64) numHalfRows)); 525 581 582 psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s); 526 583 // 527 584 // Iterate over all pixels surrounding this stamp. … … 534 591 // based on whether the kernels are ISIS or POIS. 535 592 593 psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x); 536 594 if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { 537 595 // … … 547 605 // First convolution. This will set the value for the stampVector. 548 606 // 607 // 608 // XXX: verify the [y-v2][x-u2] subscript. This generated errors in 609 // testing, depending on kernel size and footprint. 610 // 549 611 psF32 conv1 = polyValues->data.F64[j1][i1] * reference->data.F32[y - v1][x - u1]; 550 612 … … 561 623 // 562 624 for (psS32 k2 = k1; k2 < numKernels; k2++) { 563 psS32 u2 = kernels->u->data.F32[k2]; // Offset in x 564 psS32 v2 = kernels->v->data.F32[k2]; // Offset in y 565 psS32 i2 = kernels->xOrder->data.F32[k2]; // Polynomial order in x 566 psS32 j2 = kernels->yOrder->data.F32[k2]; // Polynomial order in y 567 625 psS32 u2 = (psS32) kernels->u->data.F32[k2]; // Offset in x 626 psS32 v2 = (psS32) kernels->v->data.F32[k2]; // Offset in y 627 psS32 i2 = (psS32) kernels->xOrder->data.F32[k2]; // Polynomial order in x 628 psS32 j2 = (psS32) kernels->yOrder->data.F32[k2]; // Polynomial order in y 629 // 630 // XXX: verify the [y-v2][x-u2] subscript. This generated errors in 631 // testing, depending on kernel size and footprint. 632 // 633 //printf("footprint is %d\n", footprint); 634 //printf("Stamp (%d, %d).\n", stamp->y, stamp->x); 635 //printf("HERE (y, x) is (%d, %d). (v2, u2) is (%d, %d).\n", y, x, v2, u2); 636 //printf("HERE 00 (%d %d) (%d %d)\n", j2, i2, y-v2, x-u2); 568 637 // 569 638 // Second convolution 570 639 // 571 psF32 conv2 = polyValues->data.F64[j2][i2] * reference->data.F32[y - v2][x - u2]; 640 psF32 conv2 = polyValues->data.F64[j2][i2] * 641 reference->data.F32[y-v2][x-u2]; 572 642 573 643 // … … 600 670 601 671 } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { 602 printf("XXX: put some warning message here (ISIS kernels not implemented).\n");603 return(false);672 // printf("XXX: put some warning message here (ISIS kernels not implemented).\n"); 673 // return(false); 604 674 // XXX: HEY: code this 605 /*606 675 for (psS32 k1 = 0; k1 < numKernels; k1++) { 607 psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1 );676 psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1, x, y); 608 677 609 678 for (psS32 k2 = k1; k2 < numKernels; k2++) { 610 psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2 );679 psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2, x, y); 611 680 612 681 stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2; 613 682 614 683 } 615 vector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2;684 stampVector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2; 616 685 stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2; 617 686 } 618 */619 687 } else { 620 688 printf("XXX: put some warning message here (bad kernel->type).\n"); … … 662 730 663 731 /******************************************************************************* 732 XXX: Assert correct vector matrix sizes. 664 733 ******************************************************************************/ 665 734 psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL … … 668 737 { 669 738 PS_ASSERT_PTR_NON_NULL(stamps, NULL); 739 PS_ASSERT_IMAGE_NON_NULL(((pmStamp *) stamps->data[0])->matrix, NULL); 740 PS_ASSERT_VECTOR_NON_NULL(((pmStamp *) stamps->data[0])->vector, NULL); 670 741 psS32 size = ((pmStamp *) stamps->data[0])->vector->n; 671 742 … … 804 875 // Iterate over the kernel basis functions 805 876 for (psS32 k = 0; k < nBF; k++) { 806 // XXX: What's the story with this?807 #if 0808 877 psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; 809 878 psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; 810 879 psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; 811 #else 812 813 psF32 polyVal = 1.0; 814 #endif 815 816 psImage *preCalc = (psImage *) kernels->preCalc->data[nBF]; 880 881 psImage *preCalc = (psImage *) kernels->preCalc->data[k]; 817 882 for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) { 818 883 for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) { … … 820 885 conv += solution->data.F64[k] * 821 886 input->data.F32[yy+row][xx+col] * 822 preCalc->data.F32[yy -kernelSize][xx-kernelSize] *887 preCalc->data.F32[yy+kernelSize][xx+kernelSize] * 823 888 polyVal; 824 889 } … … 854 919 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); 855 920 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); 856 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); 921 if (kernels->preCalc != NULL) { 922 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); 923 } else { 924 if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { 925 psError(PS_ERR_BAD_PARAMETER_NULL, true, 926 "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); 927 return(NULL); 928 } 929 } 857 930 psS32 nBF = kernels->u->n; 858 931 PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); … … 943 1016 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); 944 1017 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); 945 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); 1018 if (kernels->preCalc != NULL) { 1019 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); 1020 } else { 1021 if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { 1022 psError(PS_ERR_BAD_PARAMETER_NULL, true, 1023 "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); 1024 return(NULL); 1025 } 1026 } 946 1027 psS32 nBF = kernels->u->n; 947 1028 PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL); … … 1026 1107 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false); 1027 1108 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false); 1028 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); 1109 if (kernels->preCalc != NULL) { 1110 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false); 1111 } else { 1112 if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { 1113 psError(PS_ERR_BAD_PARAMETER_NULL, true, 1114 "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); 1115 return(false); 1116 } 1117 } 1118 1029 1119 psS32 nBF = kernels->u->n; 1030 1120 PS_ASSERT_VECTOR_SIZE(solution, nBF+1, false); … … 1091 1181 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL); 1092 1182 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL); 1093 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); 1183 if (kernels->preCalc != NULL) { 1184 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL); 1185 } else { 1186 if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) { 1187 psError(PS_ERR_BAD_PARAMETER_NULL, true, 1188 "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n"); 1189 return(NULL); 1190 } 1191 } 1094 1192 PS_ASSERT_INT_EQUAL(1+kernels->u->n, solution->n, NULL); 1095 1193 … … 1097 1195 psS32 spatialOrder = kernels->p_spatialOrder; 1098 1196 psS32 kernelSize = kernels->p_size; 1099 psS32 xCenter; // The pixel location for the center of the kernel in out img.1100 psS32 yCenter; // The pixel location for the center of the kernel in out img.1101 psS32 numRows;1102 psS32 numCols;1103 psS32 imgSize;1104 1197 1105 1198 if (out != NULL) { 1106 xCenter = (psS32) ( ((x+1.0)/2.0) * (psF32) (out->numCols)); 1107 yCenter = (psS32) ( ((y+1.0)/2.0) * (psF32) (out->numRows)); 1108 if ( ((xCenter - kernelSize) < 0) || 1109 ((xCenter + kernelSize) > out->numCols) || 1110 ((yCenter - kernelSize) < 0) || 1111 ((yCenter + kernelSize) > out->numRows)) { 1112 1199 if ((out->numCols < (1+2*kernelSize)) || (out->numRows < (1+2*kernelSize))) { 1113 1200 printf("XXX: generate WARNING: out image is not large enough.\n"); 1114 1201 return(out); 1115 1202 } 1116 1203 } else { 1117 // 1118 // We calculate the minimize size image that can contain the kernel 1119 // centered at the specified (x, y) position. 1120 // 1121 numRows = (psS32) ceilf( ((psF32) kernelSize) / ((y + 1.0)/2.0)); 1122 numRows = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((y + 1.0)/2.0)))); 1123 numCols = (psS32) ceilf( ((psF32) kernelSize) / ((x + 1.0)/2.0)); 1124 numCols = PS_MAX(numRows, (((psF32) kernelSize) / (1.0 - ((x + 1.0)/2.0)))); 1125 imgSize = PS_MAX(numRows, numCols); 1126 1127 out = psImageAlloc(imgSize, imgSize, PS_TYPE_F32); 1128 xCenter = (psS32) ( ((x+1.0)/2.0) * (psF32) (out->numCols)); 1129 yCenter = (psS32) ( ((y+1.0)/2.0) * (psF32) (out->numRows)); 1204 out = psImageAlloc(1+2*kernelSize, 1+2*kernelSize, PS_TYPE_F32); 1130 1205 } 1131 1206 PS_IMAGE_SET_F32(out, 0.0); 1207 1208 // 1209 // Generate the spatial-order polynomial. The [i][j]-th element of 1210 // the psImage polyValues will hold (x^i * y^j) for the stamp. 1211 // 1212 psImage *polyValues = GenSpatialOrder(spatialOrder, x, y); 1132 1213 1133 1214 // XXX: switch (i, j) so they correspond to (x, y). 1134 1215 if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) { 1135 1216 for (psS32 k = 0 ; k < nBF ; k++) { 1136 // XXX: verify these loop bounds with the size of the preCalc images. 1217 psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; 1218 psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; 1219 psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; 1220 1221 // XXX: Verify that this is correct. 1137 1222 for (psS32 i = -kernelSize ; i <= kernelSize ; i++) { 1138 1223 for (psS32 j = -kernelSize ; j <= kernelSize ; j++) { 1139 psImage *preCalc = (psImage *) kernels->preCalc->data[nBF]; 1140 out->data.F32[yCenter + i][xCenter + j]+= solution->data.F64[nBF] * 1141 preCalc->data.F32[i+kernelSize][j+kernelSize]; 1224 psImage *preCalc = (psImage *) kernels->preCalc->data[k]; 1225 out->data.F32[i+kernelSize][j+kernelSize]+= 1226 solution->data.F64[k] * 1227 preCalc->data.F32[i+kernelSize][j+kernelSize] * 1228 polyVal; 1142 1229 } 1143 1230 } 1144 1231 } 1145 1232 } else if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) { 1146 //1147 // Generate the spatial-order polynomial. The [i][j]-th element of1148 // the psImage polyValues will hold (x^i * y^j) for the stamp.1149 //1150 psImage *polyValues = GenSpatialOrder(spatialOrder, x, y);1151 1152 1233 for (psS32 k = 0 ; k < nBF ; k++) { 1153 // XXX: Why don't we have compilation warnings on type here? 1154 psS32 u = kernels->u->data.F32[nBF]; 1155 psS32 v = kernels->v->data.F32[nBF]; 1156 psS32 xOrder = kernels->xOrder->data.F32[nBF]; 1157 psS32 yOrder = kernels->yOrder->data.F32[nBF]; 1234 // XXX: Why don't we have compilation warnings on type here (if 1235 // we remove the (psS32) cast)? 1236 psS32 u = (psS32) kernels->u->data.F32[k]; 1237 psS32 v = (psS32) kernels->v->data.F32[k]; 1238 psS32 xOrder = (psS32) kernels->xOrder->data.F32[k]; 1239 psS32 yOrder = (psS32) kernels->yOrder->data.F32[k]; 1240 psF32 polyVal = polyValues->data.F32[yOrder][xOrder]; 1158 1241 // XXX: Verify that this is correct. 1159 out->data.F32[yCenter - (psS32) v][xCenter - (psS32) u]+= solution->data.F64[nBF] * 1160 polyValues->data.F64[yOrder][xOrder]; 1161 } 1162 psFree(polyValues); 1163 } 1242 1243 out->data.F32[kernelSize - v][kernelSize - u]+= 1244 solution->data.F64[k] * 1245 polyValues->data.F64[yOrder][xOrder] * 1246 polyVal; 1247 polyVal = polyVal; 1248 } 1249 } 1250 psFree(polyValues); 1164 1251 1165 1252 return(out);
Note:
See TracChangeset
for help on using the changeset viewer.
